A new scripting library for modeling flow and transport in fractured rock with channel networks

Accepted Manuscript A new scripting library for modeling flow and transport in fractured rock with channel networks Benoît Dessirier, Chin-Fu Tsang, Auli Niemi PII:

S0098-3004(17)30456-9

DOI:

10.1016/j.cageo.2017.11.013

Reference:

CAGEO 4054

To appear in:

Computers and Geosciences

Received Date: 20 April 2017 Revised Date:

6 November 2017

Accepted Date: 11 November 2017

Please cite this article as: Dessirier, Benoî., Tsang, C.-F., Niemi, A., A new scripting library for modeling flow and transport in fractured rock with channel networks, Computers and Geosciences (2017), doi: 10.1016/j.cageo.2017.11.013. 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

3 4

Benoît Dessiriera,*, Chin-Fu Tsanga,b, Auli Niemia

6

*

corresponding author

7

a

Uppsala University, Villavägen 16, Uppsala, Sweden

8

b

Lawrence Berkeley National Laboratory, Berkeley, CA, USA

9

[email protected]

10

[email protected]

11

[email protected]

M AN U

5

RI PT

2

Manuscript: A new scripting library for modeling flow and transport in fractured rock with channel networks

SC

1

12

Abstract:

14

Deep crystalline bedrock formations are targeted to host spent nuclear fuel owing to their overall low

15

permeability. They are however highly heterogeneous and only a few preferential paths pertaining to a

16

small set of dominant rock fractures usually carry most of the flow or mass fluxes, a behavior known as

17

channeling that needs to be accounted for in the performance assessment of repositories. Channel

18

network models have been developed and used to investigate the effect of channeling. They are usually

19

simpler than discrete fracture networks based on rock fracture mappings and rely on idealized full or

20

sparsely populated lattices of channels. This study reexamines the fundamental parameter structure

21

required to describe a channel network in terms of groundwater flow and solute transport, leading to an

22

extended description suitable for unstructured arbitrary networks of channels. An implementation of

23

this formalism in a Python scripting library is presented and released along with this article. A new

24

algebraic multigrid preconditioner delivers a significant speedup in the flow solution step compared to

AC C

EP

TE D

13

1

ACCEPTED MANUSCRIPT

previous channel network codes. 3D visualization is readily available for verification and interpretation

26

of the results by exporting the results to an open and free dedicated software. The new code is applied

27

to three example cases to verify its results on full uncorrelated lattices of channels, sparsely populated

28

percolation lattices and to exemplify the use of unstructured networks to accommodate knowledge on

29

local rock fractures.

30

Keywords: channel network, groundwater flow, solute transport, graph theory

SC

31

RI PT

25

1. Introduction

Groundwater flow and solute transport in deep crystalline bedrock are central to the safety assessment

33

of deep underground applications such as geological waste disposal (Pusch, 2009; Tsang et al., 2015).

34

These deep formations have generally undergone some degree of fracturing during their geological

35

history and exhibit a highly heterogeneous and anisotropic permeability field related to the embedded

36

network of rock fractures.

37

The mean transmissivity varies widely between fractures at the scale of the network (Dverstorp et al. ,

38

1992). Consequently, the large-scale flow through a network of fractures is usually focused on a limited

39

number of channels formed by a small subset of fractures with high transmissivities (Neretnieks et al.,

40

1982; Tsang and Neretnieks, 1989; Tsang and Tsang, 1989), sometimes called the backbone of the

41

network (e.g. Aldrich et al., 2017). Besides this large-scale channeling effect, in each flowing fracture,

42

the local aperture variability, in particular under confining stress, has also been shown to lead to

43

preferential flow pathways, resulting in flow channeling, this time at the scale of a single fracture (de

44

Dreuzy et al., 2012; Ishibashi et al., 2015; Kang et al., 2016; Larsson et al., 2012; Tsang and Tsang, 1989).

45

Flow channeling has been proposed as a direct cause for anomalous transport (Nowamooz et al., 2013;

46

Kang et al., 2016; Tsang and Neretnieks, 1998) and makes it difficult as well to precisely estimate the

47

flow-wetted surface, i.e. the surface of rock in contact with actively flowing groundwater, which directly

AC C

EP

TE D

M AN U

32

2

ACCEPTED MANUSCRIPT

regulates solute retardation through diffusion and sorption inside the rock matrix (Gylling et al., 1999;

49

Neretnieks, 1987), two other processes leading to anomalous transport (Carrera et al., 1998).

50

Conceptual models for groundwater flow and transport in fractured rock are built in practice by

51

upscaling the hydraulic rock properties considering different data/information supports (Selroos et al.,

52

2002). Stochastic continua assume permeability values distributed on a volumetric grid thus assuming

53

full connectivity within the grid blocks (Tsang et al., 1996). Discrete fracture networks (DFN) use thin

54

plates as the basic flowing units (e.g. de Dreuzy et al., 2012; Dverstorp et al., 1992; Hyman et al., 2015).

55

Channel networks consider linear, possibly tortuous, preferential paths to be the representative flowing

56

units (Gylling et al., 1999; Moreno and Neretnieks, 1993). This latter approach relies traditionally on a

57

lattice structure with binary or lognormally distributed lattice bond properties to calculate a steady state

58

flow solution (Black et al. 2016; Figueiredo et al., 2016; Margolin et al. 1998; Moreno and Neretnieks,

59

1993). The solute transport is then calculated often by a particle-tracking technique and assuming a one-

60

dimensional advection equation in the channels with retardation by diffusion and sorption in the

61

adjacent rock matrix (Moreno and Neretnieks, 1993). Different conceptualizations of the rock matrix, in

62

terms of thickness reachable by transport and layering, can lead to analytical solutions in the Laplace

63

domain (Mahmoudzadeh et al., 2013; Moreno and Crawford, 2009) and sometimes directly in the time

64

domain (Moreno and Neretnieks, 1993), which keeps the computational cost low and enables rapid

65

calculations even on a regular desktop computer.

66

The channel network approach has been successfully applied to describe tracer test experiments at

67

kilometer scale, such as the long-term pumping test (LPT2) at the Äspö Hard Rock Laboratory in

68

southeastern Sweden (Gylling et al., 1998), and at the tunnel scale (10-100 meters), such as the tracer

69

retention understanding experiment (TRUE) also at Äspö (Gylling et al., 1999). More recent

70

developments, have been focused on the quantification of the impact of solute diffusion into stagnant

71

water zones around the flowing channels (Mahmoudzadeh et al., 2013; Shahkarami et al., 2016) and the

AC C

EP

TE D

M AN U

SC

RI PT

48

3

ACCEPTED MANUSCRIPT

treatment of whole decay chains of tracers (Mahmoudzadeh et al., 2014, 2016; Shahkarami et al., 2015).

73

Black et al. (2016) also used percolation networks based on sparsely populated lattices of channels to

74

propose an explanation to the so-called positive skin effect around tunnels, i.e. a reduced tunnel inflow

75

compared with radially convergent (two-dimensional) flow which is usually interpreted as an additional

76

head drop near the tunnel wall. They showed that sparse networks of long channels just above the

77

percolation threshold exhibit behavior consistent with the observed positive skin effect (Black et al.,

78

2016).

79

Although very versatile and physically consistent with the observed channeling effect, the channel

80

network approach comes at the cost of some level of abstraction because channel network properties

81

(lattice spacing, conductance, aperture, flow wetted surface) can only indirectly be inferred by the

82

number of inflow points on tunnel walls or transmissivity distributions in borehole segments (Gylling et

83

al., 1999; Neretnieks, 1987). First order estimates of the parameter distributions for the channel

84

network model are usually obtained by considering additional assumptions, e.g. a cubic lattice

85

arrangement of the channels (Moreno and Neretnieks, 1993). Channel conductances and lattice spacing

86

have been linked for example to transmissivity distributions in borehole intervals (Moreno and

87

Neretnieks, 1993). The estimation rule for flow wetted surface from Gylling et al. (1999) relies however

88

on channels being uniformly oriented in space. Black et al. (2016) indicate that the two basic parameters

89

underlying their network generator are related to the mean fracture trace length and fracture intensity

90

respectively but do not elaborate on potential functional relationships.

91

While it has helped to define simple generators for stochastic channel network realizations and

92

demonstrated the versatility of the channel network approach, the lattice topology has also limited the

93

flexibility and ability to include deterministic features without a cumbersome mapping procedure. In the

94

case of spatially uncorrelated conductances, the lattice spacing acts de-facto as the mean correlation

AC C

EP

TE D

M AN U

SC

RI PT

72

4

ACCEPTED MANUSCRIPT

length for the network, whereas it becomes the model resolution while mapping deterministic features

96

such as tunnels and fracture zones.

97

To generalize some of the modeling techniques routinely applied to lattice networks of channels, the

98

present study starts with a general reflection on the data structure adequate to describe an

99

unstructured channel network, i.e. a network without any specific assumption on the position of the

100

nodes (channel intersections) or the connectivity between nodes (channels). This will show that the

101

existing methods for treating lattice networks can be extended to arbitrary network topologies and

102

geometries. This paper then introduces a new open-source scripting library, Pychan3d, based on the

103

new extended formalism and suitable to represent unstructured channel networks and shows significant

104

gains in computing time obtained by using algebraic multigrid preconditioners rather than more

105

traditional ones based incomplete LU factorization. Two examples of applications are given as

106

verification cases against the existing lattice-based channel network codes of Gylling et al. (1999) and

107

Black et al. (2016). A third and final example illustrates some new capabilities for building a custom

108

channel network based on an a priori fracture geometry which could originate from a conceptual

109

geological model or from a realization of a discrete fracture network model.

SC

M AN U

TE D

2. Methods

EP

110

RI PT

95

This section presents an analysis of the basic requirement to define a channel network with regard to

112

the steady-state flow problem and the solute transport problem respectively. In each case the

113

corresponding implementation in the new software library Pychan3d is also briefly described.

114

2.1 Steady-state flow network

115

In terms of flow, a single channel between two channel intersections can be represented by its

116

conductance C, a single parameter that lumps the influences of its average conductivity K and its

AC C

111

5

ACCEPTED MANUSCRIPT

117

geometry (cross-section A and length L), through the following relationship (Moreno and Neretnieks,

118

1993):

120

×

(1).

RI PT

=

119

The flow rate through the channel Q can then be calculated by Darcy’s law: =

121

× ∆ℎ

(2)

with C the previously defined conductance and ∆h the hydraulic head difference between the

123

extremities of the channel. To completely define the steady state flow problem, one needs to know the

124

topology of the network, i.e. the connectivity map between channels, the conductance of each channel

125

and the boundary conditions.

126

From the mathematical stand point, a channel network can thus be represented by a flow graph where

127

the nodes of the graph are the channel intersections and the edges of the graph are channels

128

characterized by their conductances (Li et al., 2014). Solving the flow problem consists of solving a

129

sparse system:

M AN U

TE D

× ℎ =

(3)

EP

130

SC

122

where h is the vector of unknown heads at the channel intersections, the symmetric positive definite

132

matrix M is the Laplacian matrix of the conductance graph and has the generic term:

133 134

AC C

131

(

( ) = ∑

)

= − +

(4a)

,

×

!

(4b)

135

with Cij the conductance of the channel between nodes i and j, δDirichlet, i is an indicator function equal to

136

1 at nodes where a fixed head boundary condition is in place and equal to 0 elsewhere, Cbnd is an

6

ACCEPTED MANUSCRIPT

137

arbitrary constant with a value significantly larger than all conductances in the system, and b is a sparse

138

vector of boundary conditions with general term: =

,

×

!

×ℎ

!,

+

" #$% ,

×&

!,

(5)

RI PT

139

where δNeuman, i is an indicator function for nodes with a prescribed flow boundary condition, hbnd, i and

141

qbnd, i are the applied head or flow values at the boundaries. Simultaneously prescribing head and flux in

142

this case yields a Robin boundary condition (3rd kind) in which case Cbnd is a weight factor that

143

corresponds to the conductance of the insulating layer.

144

It appears clearly from this graph representation that the lattice topology used by previous channel

145

networks (e.g. Moreno and Neretnieks, 1993) is a means to generate and parametrize stochastic

146

networks but is not a necessary assumption to build the graph representation of channel networks.

147

Once a flow solution has been obtained and a groundwater head value has been associated to each

148

channel intersection, the graph defined by the channel intersections as nodes and the channels oriented

149

by the flow direction as edges is a directed acyclic graph (DAG) over which the values of groundwater

150

heads at the nodes provide a partial topological order.

151

The new code Pychan3d is developed as a Python package with dependencies on the Numpy and Scipy

152

modules to allow easy portability and rapid code development (Figure 1). A scripting approach is

153

adopted, whereby the user builds a script, i.e. a list of consecutive operations of data or model

154

management. This approach is in line with a current trend in computational geosciences (e.g. Bakker et

155

al. 2016; Croucher, 2011). The scripting approach has the advantage to preserve an integral trace of the

156

modeling procedure as well as shortening the development cycle by avoiding the huge effort that is the

157

development and maintenance of a decent quality graphical user interface. Unstructured flow networks

158

are implemented as a main Python class. To initialize a network, one needs an array of nodes (size Nn)

159

and a sparse matrix of dimension (NnxNn) with Nc non-zero entries each corresponding to a channel. The

AC C

EP

TE D

M AN U

SC

140

7

ACCEPTED MANUSCRIPT

entry on row i and column j is the conductance between nodes i and j. The nodes are represented by

161

their three-dimensional coordinates, which can be used for visualization purposes via an optional

162

dependency on the vtk package, but the node coordinates have no bearing on the flow solution. Two

163

additional sparse arrays (of dimension Nn) corresponding to the Dirichlet and Neuman boundary

164

conditions at each node are required to completely define the flow problem. Two subclasses are

165

provided for convenience and for testing that correspond respectively to the full lattice network

166

parametrization (Gylling et al., 1999; Moreno and Neretnieks, 1993) and the sparse lattice network

167

parametrization (Black et al., 2007; Margolin et al., 1998). In each case, the lattice network is generated

168

and then internally translated into the generic form previously described.

EP

169

TE D

M AN U

SC

RI PT

160

Figure 1: Modeling workflow and dependencies of Pychan3d.

171

Irrespective of the network type, the flow solution can be obtained by several types of solvers (Figure 1,

172

Table 1). The first one is available through Scipy as a wrapper of the serial SuperLU library, which

173

provides a direct solver. It is unconditionally stable and suitable for small network (Nn and Nc < 105),

174

which is interesting for testing and for sparse networks. The second solver is based on an additional

175

package Pyamg which offers Algebraic Multi-Grid (AMG) solvers and preconditioners (Bell et al., 2013). A

176

classic AMG preconditioner is built using the system matrix M and then used to solve the flow problem

AC C

170

8

ACCEPTED MANUSCRIPT

with a preconditioned conjugate gradient solver from the Scipy library. On a single core, this new

178

preconditioner is shown to greatly reduce the time required to obtain the flow solution in comparison

179

with the more commonly used preconditioner based on an incomplete LU (ILU) factorization used by

180

Gylling et al. (1999) (see upcoming section 3.1 Full Lattice Network), a result consistent with the

181

observations of Dreuzy et al. (2013). Pychan3d has also been successfully linked to the PETSc library

182

(Balay et al., 2016) through Python bindings available via two additional packages, namely mpi4py and

183

petsc4py (Dalcin et al, 2011). PETSc includes parallel algorithms for building block ILU or AMG

184

preconditioners and parallel iterative solvers. Adequate graph partitioning is available thanks to the

185

Pymetis bindings to the program METIS (Karypis and Kumar, 1998).

186

One set of utilities is provided to carve different domain shapes such as spheres, cylinders or arbitrary

187

convex shapes to represent tunnels or outer domain shapes and help assign corresponding boundary

188

conditions (Table 1). Another utility is also available to identify percolating clusters of channels between

189

two sets of nodes prior to calculating the flow solution, which allows to trim disconnected channels as

190

well (Table 1). Networks and their corresponding flow solution or percolating clusters can easily be

191

exported to vkt file formats for quick visualization and control of the output with the external open

192

visualization software Paraview (Ahrens et al, 2005; Table 1).

193

Table 1: comparison of code functionalities for CHAN3D (Gylling et al., 1999), Hyperconv (Black et al., 2006) and Pychan3d.

AC C

EP

TE D

M AN U

SC

RI PT

177

Functionality

CHAN3D

HyperConv

Pychan3d

(Gylling et al.,

(Black et al.,

This article

1999)

2016)

Supported network types -

Full Lattice networks

Yes

Yes

Yes

-

Sparse Lattice networks

No

Yes

Yes

9

ACCEPTED MANUSCRIPT

Unstructured networks

No

No

Yes

Percolation-check/trimming

Not relevant

Yes

Yes

Independent clusters analysis

Not relevant

Yes

Yes

Carving utilities (cylinder, sphere, convex)

Presumed No

Yes

Yes

XYZ-plane boundary conditions

Yes

Yes

Distributed boundary conditions

Yes,

Presumed.

Steady-state flow solution:

-

-

Direct solver

No

Serial iterative preconditioned

Unknown

Yes,

TE D

-

M AN U

needed if #>50.

Yes Yes

SC

recompilation

RI PT

-

Yes (Serial, small networks: Nc < 105)

Yes,

Yes,

conjugate gradient solver

ILU precond.

ICC precond.

AMG precond.

Parallel iterative solvers and

No

No

Yes, if PETSc installed, optional

AC C

EP

preconditioners

graph partitioning with METIS

Transport (Particle-tracking): -

Mixing rules at intersections

Perfect mixing

N/A

Perfect mixing

Channel types

Plug flow

N/A

Plug flow

Sorption/Diffusion

Sorption/Diffusion

in infinite matrix

in infinite matrix

10

ACCEPTED MANUSCRIPT

Sorption/Diffusion

in finite matrix

in finite matrix

Sorption/Diffusion

(under test)

in layered matrix

Sorption/Diffusion

RI PT

Sorption/Diffusion

in layered matrix (under test)

-

Flow solution

Partial, MATLAB

SC

3D Visualization: Yes,

Full, Paraview

M AN U

dedicated program.

-

Particle-tracking

Partial, MATLAB

194

No transport

Full, Paraview

2.2 The transport problem

196

The solute transport problem over a flow network can be broken down in two parts: a mixing rule at

197

each node and a transport model in each channel. The simplest mixing rule is that of perfect mixing at all

198

nodes, which means that a hypothetical particle of solute gets to leave a node through a certain exit

199

channel with a probability equal to the outflow rate to this channel divided by the total throughflow at

200

that node (Berkowitz et al., 1994). Stream tube routing, in contrast, matches stream tubes between

201

entry and exit channels deterministically and mixing occurs only by diffusion between the stream tubes

202

at the intersection or later by diffusion along the exit channels if they result from the merger of stream

203

tubes from different entry channels (Berkowitz et al., 1994). Stream tube routing and perfect mixing are

204

the two ends of a spectrum of intermediate mixing behaviors. Using four-member orthogonal

205

intersections, Berkowitz et al. (1994) showed that the applicable mixing rule should theoretically be

206

decided according to the effective Peclet number at each intersection. Furthermore, based on two-

AC C

EP

TE D

195

11

ACCEPTED MANUSCRIPT

dimensional square lattice networks, Kang et al. (2015) showed that the choice of the mixing rule can

208

have a significant impact on transverse dispersion in networks with low levels of heterogeneity, but

209

becomes insignificant for high levels of heterogeneity. Using two-dimensional networks with power-law

210

length distributions, Park et al. (2001) showed that the impact of the mixing rule is negligible in

211

networks with a low average coordination number (number of intersecting members at each node),

212

which is typically the case for sparse channel networks.

213

The transport model depends on the transport processes at work in the channel which is likely a

214

function of the channel geometry, but it may also depend on the layout and properties of the rock mass

215

around the channel and the type of particle being transported. In practical terms the transport model

216

determines a travel time probability distribution P(t) for the particle through this channel. The simplest

217

behavior would be to consider an inert particle that is only subject to advection within the channel in

218

which case the travel time distribution is:

M AN U

SC

RI PT

207

'(() = ()/ )

(6)

TE D

219

with δ the dirac delta function, V the channel volume and Q the channel flow rate. Another

221

configuration lending itself to an analytic solution is to consider that the particle is subject to advection

222

along the channel opening, diffusion into the pores of the adjacent rock matrix (perpendicularly to the

223

channel and with a matrix considered infinitely thick) and sorption on the channel wall and surfaces

224

inside the rock matrix (Moreno and Neretnieks, 1993). The corresponding travel time distribution is:

226

AC C

225

EP

220

'(() =

./0 12

×

'(() = 0 for ( ≤ - × )/

345

67(

89×:/2);

./0 1 ) 2

× exp (−(

×

(7a)

345 ? ) @( 89×:/2)

for ( > - × )/

(7b)

227

with R the surface retardation factor, V the volume of the channel, Q the channel flow rate (obtained

228

from the flow solution), FWS the flow-wetted surface (which can be decomposed as twice the product 12

ACCEPTED MANUSCRIPT

229

of the channel length L by the mean channel width W) and MPG the material property group that

230

describes the joint retardation due to diffusion and sorption in the rock matrix and that can be

231

decomposed as

232

matrix sorption coefficient and E the bulk density of the rock (Moreno and Crawford, 2009; Moreno

233

and Neretnieks, 1993). More elaborate configurations, that offer analytical solutions in the Laplace

234

space, can address more complexity such as a layered matrix with a finite extent (Mahmoudzadeh et al.,

235

2013, 2014, 2016; Moreno and Crawford, 2009), diffusion into stagnant water zones along the channel

236

(Mahmoudzadeh et al., 2013; Shahkarami et al., 2016).

237

As shown by the aforementioned relations for the travel time probabilities, solving the solute transport

238

problem requires more input or assumptions regarding the geometry of the channels and the properties

239

of the rock. At a minimum, the volume V of each channel is required to calculate the advective travel

240

time. If matrix diffusion and sorption are considered, one must also know the flow wetted surface FWS

241

and the material property group MPG. Transport is implemented in the new Python library by a particle

242

tracking technique and by the injection of a predefined number of particles. Each particle can be

243

assigned an entry node in the network between a list of specified injection nodes either uniformly ---

244

which is sometimes called resident injection --- or following a flux-weighted rule --- flux injection

245

(Frampton and Cvetkovic, 2009; Kang et al., 2017; Dagan, 2016). Resident injection corresponds to a

246

source that would release particles at a given fixed rate whereas flux injection is typical of a large source

247

operating at constant head (Dagan, 2016). Because it samples more densely areas of the source with

248

slowly flowing channels, resident injection produces longer tails for the travel time distribution than flux

249

injection (Frampton and Cvetkovic, 2009; Kang et al., 2017). This impact of the injection mode on

250

particle dispersion in the network increases with the degree of heterogeneity on the system (Frampton

251

and Cvetkovic, 2009). A piecewise linear concentration-time relationship for tracer injection can also be

252

used to assign an injection time to each particle, consistently with an injection profile from the field for

∗ !

× E with C the effective diffusivity of the rock matrix,

∗ !

the bulk

AC C

EP

TE D

M AN U

SC

RI PT

'B = 6C ×

13

ACCEPTED MANUSCRIPT

example. The perfect mixing rule is implemented in Pychan3d at each node by a discrete random

254

variable generator giving the next visited node based on the calculated channel in- and outflow rates

255

and boundary conditions. The travel times along each channel segment are generated by randomly

256

sampling the selected travel time distribution based on the channel geometry and adjacent rock

257

parameters. As a result, each simulated particle is represented by an array of nodes through which it has

258

passed and an array of the corresponding times of passage. Routines are also included to export the

259

particles to VTK files for rapid visualization with Paraview (Ahrens et al., 2005) and diagnostics of the

260

particle injections (Table 1 and Supplementary material).

SC

M AN U

261

RI PT

253

3. Verification cases 3.1 Full lattice network (example 1)

263

As a first example, we used the newly implemented code to reproduce the flow solution over a lattice

264

network as parametrized in the predecessor CHAN3D (Gylling et al., 1999). The test domain is a cubic

265

lattice with a spacing of 1.5m and 60 nodes in each direction. Dirichlet boundary conditions (fixed head

266

values) are applied to two opposite faces of the lattice to obtain a unit average hydraulic gradient and

267

the four other faces are set to no-flow. This results to a network with 216,000 nodes, 637,200 channels

268

and 7,200 Dirichlet boundary conditions. The conductances were assumed to be lognormally distributed

269

and uncorrelated in space, so

270

example, we used the following values:

271

The flow solution was calculated for σ=0, which corresponds to uniform flow, and 10 stochastic

272

realizations of the channel network were built and solved for flow for σ=1 and σ=2 (Figure 2). Each

273

realization was solved with the original CHAN3D code, and with the two main solvers available in

274

Pychan3d (Table 2). The direct solver provides the accurate solution. The accuracy of the other solvers

275

was judged acceptable as the maximum error for the head solution relative the direct solution was kept

AC C

EP

TE D

262

=

F

× 10"×H , where N is a normal random variable. For this F

= 3.0 × 108KF

14

1

/L and values of σ ranging from 0 to 2.

ACCEPTED MANUSCRIPT

under 10-5 (Table 2). A comparison of the average solving times for similar levels of final residuals

277

obtained from the different solvers over each set of realizations are also presented in Table 2. The

278

results show that the direct solver from the Scipy library is not time efficient for a network of this size.

279

Separate tests showed that the solving time was kept under a minute for networks with less than 105

280

channels on the test desktop computer (Table 2). The iterative solver with the classical AMG

281

preconditioner, however, can handle large networks and deliver a significant speedup over CHAN3D,

282

both running in serial mode (Table 2).

TE D

M AN U

SC

RI PT

276

283

Figure 2: Flow solution over Lattice networks for uniform channel conductances (left) and two realizations of conductances with

285

a lognormal distribution - σ=1 (center) and σ=2 (right). Note that a higher value of σ increases the contrast between the

286

channels with high and low flow rates (increased channeling).

287

Table 2: Mean total time for assembly and solve with different codes and solvers (for 10 realizations, tested on the same desktop

288

computer with processor Intel Xeon quadcore 4x3.30 GHz and 16 GB RAM). All runs converge and yield a solution with a

289

maximum relative error for the hydraulic head under 10 compared with the solution from the direct solver.

Code

AC C

EP

284

-5

CHAN3D

Pychan3d-direct

Pychan3d-iterative

Solver library

SLATEC

Scipy (SuperLU)

Pyamg+Scipy

σ=0

7.7 s

10 min

1.9 s

15

ACCEPTED MANUSCRIPT

11.1 s

9 min

3.5 s

σ=2

44.8 s

9 min

23.6 s

SC

RI PT

σ=1

M AN U

290

Figure 3: Flow through percolating clusters of channels (color scale) in sparse lattice network realizations for L=10 and PON=0.02

292

(left), PON=0.04 (center) and PON=0.08 (right). The non-percolating channel members are represented by semi-transparent black

293

segments. Notice the different scales for the flow rates.

294

3.2 Sparse lattice networks (example 2)

295

This example attempts a verification against the results on sparse lattice networks reported by Black et

296

al. (2016). The modeling domain is a cubic percolation lattice of binary conductive/non-conductive

297

bonds, populated by the following two probabilistic rules: considering two consecutive channels along

298

one of the lattice directions, the second one has a probability PN to be conductive if the previous one

299

was non-conductive and a probability PA to be open if the previous one was already conductive (Black et

300

al., 2016). This gives a corresponding bond density, or probability for a channel to be conductive,

301

PON=PN/(1-PA+PN) and an average length L=1/(1-PA) of consecutive conductive channels. Two parameters,

302

for example (L, PON), suffice to define the network generator. We used here the same domain

303

dimensions as Black et al. (2016) and as the previous example, with a lattice spacing of 1.5m and 60

304

nodes in each direction in three dimensions. The channel network generation was tested for mean

305

channel length values L ranging from 1.33 to 20 lattice spacings and for densities PON ranging from 10-3

AC C

EP

TE D

291

16

ACCEPTED MANUSCRIPT

to 1 (Figure 3). Each conductive bond is assigned a conductance value of 3.0x10-10 m2/s as in Black et al.

307

(2016). Ensemble results were computed based on 500 realizations for each tested couple of values (L,

308

PON). Figure 4 shows the fraction of percolating realizations for the parameter values tested. For finite

309

size lattices, the percolation threshold can be approached by the average bond density at which 50% of

310

realizations percolate (Harter 2005). It is shown to agree with previous estimations of the percolation

311

threshold for single parameter percolation lattices (P=PON=PA=PN; Lorenz and Ziff, 1998; Figure 4 – star

312

symbol, P=0.2488). Figure 4 shows that percolation networks with long channels (L>1.33) have a lower

313

percolation threshold than the single parameter case, as previously reported by Black et al. (2016).

314

Although overall close to the values reported by Black et al. (2016), our estimates of the respective

315

percolation thresholds show some minor differences despite our best efforts to replicate their numerical

316

experiment (Figure 4). Part of these differences could well result from the interpolation between

317

different tested (L, PON) pairs or from the ensemble size, however further comparison between the two

318

codes would be in order to explain these differences.

319

Figure 3 illustrates the visualization capabilities for three network realizations with different bond

320

densities but the same mean channel length. One can see that the examples shown in the left and

321

central panel in Figure 3 are close to the estimated percolation threshold (Figure 4). They indeed exhibit

322

one and two sparse percolating clusters respectively. The realization shown in the right panel of Figure

323

3, however, is near the bond density where 100% of realizations percolate (Figure 4) and all flowing

324

channels have morphed into a single percolating cluster.

325

Figure 5 shows the equivalent conductivity of the lattice for the same ensembles of realizations. The

326

ensemble averages show good agreement with the values presented by Black et al. (2016). The

327

uncertainty on the ensemble permeability increases dramatically for very low-density values as the

328

number of percolating realizations reduces to a handful (Figure 5). For bond densities below 10-2, the

329

conductivity values form clusters around the values 5.65x10-14, 1.13x10-13 and 1.69x10-13 m/s which

AC C

EP

TE D

M AN U

SC

RI PT

306

17

ACCEPTED MANUSCRIPT

respectively correspond to one, two and three straight channels crossing the lattice directly from the

331

inflow face to the outflow face.

EP

332

TE D

M AN U

SC

RI PT

330

Figure 4: Percentage of percolating realizations for different parameter pairs (L, PON). Each dot is the average over 500

334

realizations. The error bars represent a 95% confidence interval estimated by resampling 1000 times from the corresponding

335

500 realizations (Efron and Tibshirani, 1994). The star shows previous estimates for the percolation threshold of a single

336

parameter percolation lattice (Lorenz and Ziff, 1998). The squares and dashed lines show values and interpolations reported by

337

Black et al. (2016) for two-parameter percolation lattices – each square representing the average of 100 realizations.

AC C

333

18

M AN U

SC

RI PT

ACCEPTED MANUSCRIPT

338 339

Figure 5: Mean value of the equivalent ensemble permeability (lines) of the percolating realizations for different parameter pairs

340

(L, PON). The left panel is reproduced from Black et al. (2016) . The right panel displays corresponding results from this study. In

341

the right panel, the crosses represent the single values obtained for each percolating realization, which shows the spread around

342

the ensemble means represented by the interpolated lines.

TE D

1

4. Connecting Discrete fracture networks to channel networks

343

(example 3)

EP

344

The graph representation underlying Pychan3d allows to add, delete or edit nodes and channels at any

346

location in the system. This last example presents a configuration that illustrates the gain in flexibility

347

with the new library and how it potentially helps to include information on rock fractures into the

348

channel network design without any mapping on a background lattice. Small Networks can be built from

349

scratch to test conceptual models consisting of a few fracture zones and boreholes. To demonstrate the

350

full range of flexibility of Pychan3d, we consider first a simple synthetic configuration, that can be

AC C

345

1

Published under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/)

19

ACCEPTED MANUSCRIPT

defined manually. It includes a cylindrical tunnel and three rock fractures in cascade (Figure 6). It is

352

assumed that the fracture intersections are very conductive, so that we represent them as a single

353

channel intersection where perfect mixing applies. For the sake of this example we assume further that

354

there are three parallel flowing channels in the fracture intersecting the tunnel (F1, Figure 6), four

355

parallel channels in the intermediate fracture (F2, Figure 6) and three channels in the fracture further

356

inside the rock mass (F3, Figure 6). These last three channels connect each to a different borehole that

357

penetrates this fracture, supplying a tracer source (shown with assorted colors: blue, red, green, Figure

358

6). We also consider an additional channel connecting to the intersection between fractures F2 and F3

359

that delivers a constant flow rate (inflowing Neumann boundary condition) and a channel connecting to

360

the intersection between fractures F1 and F2 that draws water away from the isolated network, also at a

361

constant flow rate (outflowing Neumann boundary condition). We assume that the hydraulic head is 0

362

at the tunnel wall and fixed to a single common value of 20m in the three boreholes. Conductances are

363

arbitrarily selected to achieve a significant contrast between parallel channels. For solute transport, we

364

assumed perfect mixing at the two fracture intersections, i.e. we consider for the sake of this example

365

that the fracture intersections form a very well connected hydraulic connector. This example

366

demonstrates the ability of the code to handle an arbitrary number of contributing channels at each

367

intersection, in this particular case eight channels, with seven channels proper and one Neuman

368

boundary condition at each intersection (Figure 6). In terms of transport model in each channel, we

369

assumed one-dimensional diffusion and possible sorption in the rock matrix as parametrized by Moreno

370

and Crawford (2009). For the channel lengths, we adopted the geometrical distance between

371

intersections. All channels have a width W of 0.2m and their aperture b is related to the calculated flow

372

rate through the cubic law (Witherspoon et al., 1980). The channel flow-wetted surface (FWS) is then

373

defined as FWS=2×L×W and the channel volume V as V=b×W×L. For matrix diffusion/sorption

374

parameters we show here the results for a surface retardation factor R=1 and a material property group

AC C

EP

TE D

M AN U

SC

RI PT

351

20

ACCEPTED MANUSCRIPT

MPG=2.0x10-8 m.s-1/2 which is representative of a non-sorbing tracer (matrix diffusion only) such as

376

Iodide in Äspö granite at approximately 400m depth (Moreno and Crawford, 2009). Cumulative

377

breakthrough curves for simultaneous pulse injections of 5000 particles in each of the three blue, red,

378

and green boreholes of Figure 6 are shown in Figure 7 by the matching colors (thick dashed curves). To

379

verify the accuracy of the new code, one can calculate the analytical transfer function of the network in

380

the Laplace space where the convolution of travel time distributions through successive channels

381

becomes the product of their Laplace transforms. The analytical cumulative recovery curve can then be

382

inverted back numerically to the time domain, using the De Hoog algorithm for example (De Hoog et al,

383

1982). These semi-analytical results are shown in Figure 7 as solid curves corresponding each to a

384

particle tracking dashed curve.

AC C

EP

TE D

M AN U

SC

RI PT

375

21

M AN U

SC

RI PT

ACCEPTED MANUSCRIPT

385

Figure 6: Unstructured network model example: Tunnel in gray, fracture 1 in red, fracture 2 in pink, fracture 3 in turquoise. The

387

fully-mixed intersections are represented by thin gray cylinders and the channels are depicted by their steady-state flow rates

388

(color scale). The three boreholes are also shown in blue, red and green.

AC C

EP

TE D

386

22

M AN U

SC

RI PT

ACCEPTED MANUSCRIPT

389

Figure 7: cumulative mass recovery after four different injections in the unstructured network. Blue, red and green correspond to

391

simultaneous instant injections of 5000 particles in each of the three boreholes and considering channel transport models

392

including diffusion and sorption into the rock matrix (See also animation appended as supplementary material). The black curves

393

show the outcome of an injection in the green well but neglecting rock matrix interactions (purely advective transport). The final

394

mass loss is due to the outflowing Neumann boundary condition. The particle tracking results (Pychan3d - dashed lines) are

395

plotted against a numerical inversion of the solution in the Laplace domain (semi-analytical - solid lines).

396

The comparison between semi-analytical solution and the particle-based solution delivered by Pychan3d

397

shows only minor differences (Figure 7). The partial final recovery is due to the outflow boundary

398

condition applied the intersection between fractures F1 and F2. The mass loss in this case corresponds

399

exactly to the outflow rate (Figure 7). An injection in the green borehole was also simulated for purely

400

advective flow in the channels and appears in Figure 7 as the black set of curves (solid and dashed). Here

401

also the agreement between the two methods is very good.

402

If one has a larger DFN, i.e. a network of polygonal fractures, Pychan3d includes one simple routine

403

presented by Cacas et al. (1990) to convert this discrete fracture network into a corresponding channel

AC C

EP

TE D

390

23

ACCEPTED MANUSCRIPT

network. This algorithm creates channels between the midpoints of all fracture intersection lines and

405

the center points of the fractures these connect. This simple method requires knowledge of the fracture

406

centers and average transmissivity as well as the midpoints of all fracture intersection lines. It can be

407

implemented very efficiently but potentially leads to some channels crossing other fracture intersection

408

lines without interaction which is not physically consistent. The design of an equivalent channel (or pipe)

409

network from a DFN has been and remains an active area of research (Dershowitz and Fidelibus, 1999;

410

Dershowitz et al., 1999; Xu et al., 2014). The effort has generally been to construct channel networks

411

that would yield flow results consistent with those of the complete DFN calculation but at a fraction of

412

the computational cost, and not on including channels as a physical phenomenon taking place in the

413

complex geometries of natural fracture networks. As such, the simple method of Cacas et al. (1990) is

414

provided here as a mere example of how to implement a transformation from DFN to channel network.

415

The broader vision for Pychan3d is to further develop methods for channel network design based on

416

DFN and test our understanding of fractured systems and the impact of channeling. For example, one

417

could test the impact of including fracture intersection lines as channels themselves connected within

418

each fracture planes by other channels the length, spacing and width of which could be based on

419

fracture aperture statistics (Larsson et al., 2012).

SC

M AN U

TE D

EP

5. Conclusions

AC C

420

RI PT

404

421

This study promotes a graph representation for arbitrary unstructured channel networks to model flow

422

and solute transport in deep fractured crystalline rock. The new representation is more general than

423

previous lattice-based channel networks but can still fully include them as specific cases.

424

A new Python library implementing the new graph representation for channel networks, Pychan3d, is

425

presented and applied to three different example cases. The results obtained with the new code are in

426

conformity with results from previous codes for lattice networks (example 1; Gylling et al., 1999) and 24

ACCEPTED MANUSCRIPT

consistent with published results on percolation lattice networks (example 2, Black et al., 2016). The

428

new implementation uses a more efficient AMG preconditioner to solve the steady state flow problem,

429

resulting in a speedup by a factor two to three on a single core (example 1). The new code has also been

430

successfully linked to the high-performance computing library PETSc to support parallel solvers and

431

preconditioners if very large networks need to be simulated. The new code supports output to VTK file

432

formats which allows the user to quickly visualize the results with dedicated software, e.g. Paraview

433

(examples 1, 2 and 3). The last example illustrates the enhanced capability to directly model flow and

434

transport in networks with arbitrary numbers of channel members at each intersection, which is not

435

possible with previous lattice-based models without resorting to some mapping procedure. The new

436

library grants full control to manually define or customize channel networks and it includes an example

437

of a simple automated method to transform a DFN into a channel network. The Python language and

438

the scripting approach provide both speed and flexibility while building channel networks. The AMG

439

preconditioners and/or parallel solvers allow for rapid solutions which is a prerequisite for efficient

440

comparison between alternative models, global sensitivity analysis or simulating stochastic ensembles.

441

And the flexible graph representation is instrumental in the effort to build channel networks that would

442

conform to information on known rock fractures or observed fracture statistics. This new library,

443

Pychan3d, is thus a suitable platform to develop and test new channel networks for the modeling of

444

flow and transport in the deep subsurface.

445

Acknowledgments

446

This work was supported by the Swedish Radiation Safety Authority (SSM). The authors are grateful to

447

Björn Gylling, the creator of CHAN3D, for agreeing to share the original code and granting us permission

448

to use it. Our thanks go also to the three anonymous reviewers who helped significantly improve the

449

original manuscript.

AC C

EP

TE D

M AN U

SC

RI PT

427

25

ACCEPTED MANUSCRIPT

450

References

451 Ahrens, J., Geveci, B., Law, C., 2005. ParaView: An End-User Tool for Large Data Visualization,

Visualization Handbook, Elsevier, 717pp.

RI PT

452

453 Aldrich, G., Hyman, J., Karra, S., Gable, C., Makedonska, N., Viswanathan, H., Woodring, J. and

Hamann, B., 2016. Analysis and Visualization of Discrete Fracture Networks Using a Flow

455

Topology Graph. IEEE Transactions on Visualization and Computer Graphics PP, no. 99: 1–1.

456

doi:10.1109/TVCG.2016.2582174.

SC

454

M AN U

457 Bakker, M., Post, V., Langevin, C. D., Hughes, J. D., White, J. T., Starn, J. J. and Fienen, M. N. 458

2016. Scripting MODFLOW model development using Python and FloPy. Groundwater 54,

459

no. 5 (2016): 733-739.

460 Balay S., Abhyankar, S., Adams, M.F., Brown, J., Brune, P., Buschelman, K., Dalcin, L., Eijkhout,

V., Gropp, W.D., Kaushik, D., Knepley, M.G., Curfman McInnes L., Rupp, K., Smith, B.F.,

462

Zampini, S. and Zhang, H., 2016. PETSc User’s Manual, Argonne National Laboratory, ANL-

463

95/11 - Revision 3.7, http://www.mcs.anl.gov/petsc

TE D

461

464 Bell W. N., Olson, L. N. and Schroder, J., 2013. PyAMG: Algebraic Multigrid Solvers in Python,

http://www.pyamg.org, Version 2.1

EP

465

AC C

466 Berkowitz, B., Naumann, C., Smith, L., 1994. Mass transfer at fracture intersections: An 467

evaluation of mixing models. Water Resources Research 30, 1765–1773.

468

doi:10.1029/94WR00432

469 Black, J. H., Robinson, P. C. and Barker, J. A., 2006. A Preliminary Investigation of the Concept 470

of" hyper-convergence" Using" sparse" Channel Networks. Swedish Nuclear fuel and Waste

471

Management Co. (SKB) SKB R-06-30.

26

ACCEPTED MANUSCRIPT

472 Black, J. H., Barker, J. A. and Woodman, N. D. 2007. An investigation of ‘sparse channel 473

networks'. Characteristic behaviours and their causes. No. SKB-R--07-35. Swedish Nuclear

474

Fuel and Waste Management Co. (SKB).

RI PT

475 Black, J. H., Woodman, N. D. and Barker, J. A., 2016. Groundwater Flow into Underground 476

Openings in Fractured Crystalline Rocks: An Interpretation Based on Long Channels.

477

Hydrogeology Journal, 1–19. doi:10.1007/s10040-016-1511-y.

SC

478 Cacas, M.C., Ledoux, E., de Marsily, G., Tillie, B., Barbreau, A., Durand, E., Feuga, B.,

Peaudecerf, P., 1990. Modeling fracture flow with a stochastic discrete fracture network:

480

calibration and validation: 1. The flow model. Water Resources Research 26, 479–489.

481

doi:10.1029/WR026i003p00479

M AN U

479

482 Carrera, J., Sánchez-Vila, X., Benet, I., Medina, A., Galarza, G., Guimerà, J., 1998. On matrix

diffusion: formulations, solution methods and qualitative effects. Hydrogeology Journal 6,

484

178–190. doi:10.1007/s100400050143

TE D

483

485 Croucher, A. “PyTOUGH: a Python scripting library for automating TOUGH2 simulations." In

Proceedings of New Zealand Geothermal Workshop 2011, 21 - 23 November 2011 Auckland,

487

New Zealand. 2011.

EP

486

488 Dagan, G., 2017. Solute plumes mean velocity in aquifer transport: Impact of injection and

detection modes. Advances in Water Resources, Tribute to Professor Garrison Sposito: An

490

Exceptional Hydrologist and Geochemist 106, 6–10.

491

doi:10.1016/j.advwatres.2016.09.014Dalcin, L.D., Paz, R.R., Kler, P.A. and Cosimo, A., 2011.

492

Parallel distributed computing using python. Advances in Water Resources 34(9), pp.1124-

493

1139.

AC C

489

27

ACCEPTED MANUSCRIPT

494 De Hoog, F. R., Knight, J. H. and Stokes, A. N., 1982. An Improved Method for Numerical 495

Inversion of Laplace Transforms. SIAM Journal on Applied Mathematics 3, no. 3 (1982): 357–

496

66.

RI PT

497 Dershowitz, W.S., Fidelibus, C., 1999. Derivation of equivalent pipe network analogues for three498

dimensional discrete fracture networks by the boundary element method. Water Resources

499

Research 35, 2685–2691. doi:10.1029/1999WR900118

SC

500 Dershowitz, B., Eiben, T., Follin, S., Andersson, J., 1999. SR 97 - Alternative models project.

Discrete fracture network modelling for performance assessment of Aberg (No. SKB-R--99-

502

43). Swedish Nuclear Fuel and Waste Management Co.

M AN U

501

503 Dreuzy, J.-R., Méheust, Y. and Pichot, G., 2012. Influence of Fracture Scale Heterogeneity on the 504

Flow Properties of Three dimensional Discrete Fracture Networks (DFN). Journal of

505

Geophysical Research: Solid Earth 117, no. B11. doi:10.1029/2012JB009461.

TE D

506 Dreuzy, J.-R., Pichot, G., Poirriez, B., Erhel, J., 2013. Synthetic benchmark for modeling flow in 507

3D fractured media. Computers & Geosciences, Benchmark problems, datasets and

508

methodologies for the computational geosciences 50, 59–71. doi:10.1016/j.cageo.2012.07.025

EP

509 Dverstorp, B., Andersson, J., and Nordqvist, W., 1992. Discrete Fracture Network Interpretation of

Field Tracer Migration in Sparsely Fractured Rock. Water Resources Research 28, no. 9:

511

2327–43. doi:10.1029/92WR01182.

AC C

510

512 Efron, B., and Tibshirani, R. J. 1994. An introduction to the bootstrap. CRC press. 513 Figueiredo, B., Tsang, C.-F., Niemi, A. and Lindgren, G., 2016. Review: The state-of-art of sparse 514

channel models and their applicability to performance assessment of radioactive waste

515

repositories in fractured crystalline formations. Hydrogeology Journal 24, no. 7: 1607-1622.

28

ACCEPTED MANUSCRIPT

516 Frampton, A., Cvetkovic, V., 2009. Significance of injection modes and heterogeneity on spatial

and temporal dispersion of advecting particles in two-dimensional discrete fracture networks.

518

Advances in Water Resources, Dispersion in Porous Media 32, 649–658.

519

doi:10.1016/j.advwatres.2008.07.010

RI PT

517

520 Gylling, B., Birgersson, L., Moreno, L., and Neretnieks, I., 1998. “Analysis of a Long-Term

Pumping and Tracer Test Using the Channel Network Model.” Journal of Contaminant

522

Hydrology 32, no. 3–4: 203–22. doi:10.1016/S0169-7722(97)00082-X.

SC

521

523 Gylling, B., Neretnieks, I. and Moreno, L., 1999. The Channel Network Model—A Tool for

Transport Simulations in Fractured Media - Gylling – Groundwater 37 no 3: 367-375.

M AN U

524

525 Harter, T., 2005. Finite-size scaling analysis of percolation in three-dimensional correlated binary 526

Markov chain random fields. Physical Review E 72, no. 2: 026120.

527 Hyman, J. D., Karra, S., Makedonska, N., Gable, C. W., Painter, S. L. and Viswanathan, H. S.,

2015. dfnWorks: A Discrete Fracture Network Framework for Modeling Subsurface Flow and

529

Transport. Computers & Geosciences 84: 10–19. doi:10.1016/j.cageo.2015.08.001.

TE D

528

530 Ishibashi, T., Watanabe, N., Hirano, N., Okamoto, A., Tsuchiya, N., 2015. Beyond laboratory

scale prediction for channeling flows through subsurface rock fractures with heterogeneous

532

aperture distributions revealed by laboratory evaluation. Journal of Geophysical Research:

533

Solid Earth 120, 106–124. doi:10.1002/2014JB011555

AC C

EP

531

534 Kang, P.K., Brown, S., Juanes, R., 2016. Emergence of anomalous transport in stressed rough 535

fractures. Earth and Planetary Science Letters 454, 46–54. doi:10.1016/j.epsl.2016.08.033

536 Kang, P.K., Dentz, M., Borgne, T.L., Juanes, R., 2015. Anomalous transport on regular fracture 537

networks: Impact of conductivity heterogeneity and mixing at fracture intersections. Physical

538

Review E 92, 022148. doi:10.1103/PhysRevE.92.022148

29

ACCEPTED MANUSCRIPT

539 Kang, P.K., Dentz, M., Le Borgne, T., Lee, S., Juanes, R., 2017. Anomalous transport in

disordered fracture networks: Spatial Markov model for dispersion with variable injection

541

modes. Advances in Water Resources, Tribute to Professor Garrison Sposito: An Exceptional

542

Hydrologist and Geochemist 106, 80–94. doi:10.1016/j.advwatres.2017.03.024

RI PT

540

543 Karypis, G., and Kumar, V., 1998. A Fast and High Quality Multilevel Scheme for Partitioning

Irregular Graphs. SIAM Journal on Scientific Computing 20 no 1: 359-392.

545

doi:10.1137/S1064827595287997.

SC

544

546 Larsson, M., Niemi, A. and Tsang, C.-F., 2012. A Study of Flow-Wetted Surface Area in a Single

Fracture as a Function of Its Hydraulic Conductivity Distribution. Water Resources Research

548

48, no. 1: W01508. doi:10.1029/2011WR010686.

M AN U

547

549 Li, S. C., Xu, Z. H., and Ma, G. W., 2014. A Graph-Theoretic Pipe Network Method for Water

Flow Simulation in Discrete Fracture Networks: GPNM. Tunnelling and Underground Space

551

Technology 42: 247–263. doi:10.1016/j.tust.2014.03.012.

TE D

550

552 Lorenz, C. D., and Ziff, R. M., 1998. Precise determination of the bond percolation thresholds and 553

finite-size scaling corrections for the sc, fcc, and bcc lattices. Physical Review E 57, no. 1: 230.

EP

554 Mahmoudzadeh, B., Liu, L. Moreno, L. and Neretnieks, I., 2013. Solute transport in fractured

rocks with stagnant water zone and rock matrix composed of different geological layers—

556

Model development and simulations. Water Resources Research 49, no. 3: 1709-1727.

AC C

555

557 Mahmoudzadeh, B., Liu, L. Moreno, L. and Neretnieks, I., 2014. Solute transport in a single 558

fracture involving an arbitrary length decay chain with rock matrix comprising different

559

geological layers. Journal of contaminant hydrology 164: 59-71.

30

ACCEPTED MANUSCRIPT

560 Mahmoudzadeh, B., Liu, L. Moreno, L. and Neretnieks I., 2016. Solute transport through fractured 561

rock: Radial diffusion into the rock matrix with several geological layers for an arbitrary length

562

decay chain. Journal of Hydrology 536: 133-146.

RI PT

563 Margolin, G., Berkowitz, B., and Scher, H., 1998. Structure, Flow, and Generalized Conductivity 564

Scaling in Fracture Networks. Water Resources Research 34, no. 9: 2103–21.

565

doi:10.1029/98WR01648.

SC

566 Moreno, L., and Crawford, J., 2009. Can We Use Tracer Tests to Obtain Data for Performance

Assessment of Repositories for Nuclear Waste? Hydrogeology Journal 17, no. 5: 1067–80.

568

doi:10.1007/s10040-008-0418-7.

M AN U

567

569 Moreno, L., and Neretnieks, I., 1993. Fluid Flow and Solute Transport in a Network of Channels. 570

Journal of Contaminant Hydrology 14 no. 3-4: 163-192.

571 Neretnieks, I., Eriksen, T. and Tähtinen, P., 1982. Tracer Movement in a Single Fissure in Granitic

Rock: Some Experimental Results and Their Interpretation. Water Resources Research 18, no.

573

4: 849–58. doi:10.1029/WR018i004p00849.

TE D

572

574 Neretnieks, I., 1987. Channeling effects in flow and transport in fractured rocks—Some recent

observations and models, In proceedings of GEOVAL-87 International Symposium, Swed.

576

Nucl. Power Insp. (SKI), Stockholm, April 7-9, 1987.

EP

575

AC C

577 Nowamooz, A., Radilla, G., Fourar, M., Berkowitz, B., 2013. Non-Fickian Transport in 578

Transparent Replicas of Rough-Walled Rock Fractures. Transport in Porous Media 98, 651–

579

682. doi:10.1007/s11242-013-0165-7

580 Park, Y.-J., de Dreuzy, J.-R., Lee, K.-K., Berkowitz, B., 2001. Transport and intersection mixing 581

in random fracture networks with power law length distributions. Water Resources Research

582

37, 2493–2501. doi:10.1029/2000WR000131

31

ACCEPTED MANUSCRIPT

583 Pusch, R., 2009. Geological Storage of Highly Radioactive Waste. Springer. 584 Selroos, J.-O., Walker, D. D., Ström, A., Gylling, B. and Follin, S., 2002. Comparison of

Alternative Modelling Approaches for Groundwater Flow in Fractured Rock. Journal of

586

Hydrology 257, no. 1–4:174–88. doi:10.1016/S0022-1694(01)00551-0.

RI PT

585

587 Shahkarami, P., Liu, L. Moreno, L. and Neretnieks, I., 2015. Radionuclide Migration through

Fractured Rock for Arbitrary-Length Decay Chain: Analytical Solution and Global Sensitivity

589

Analysis. Journal of Hydrology 520: 448–60. doi:10.1016/j.jhydrol.2014.10.060.

SC

588

590 Shahkarami, P., Liu L., Neretnieks, I. and Moreno, L., 2016. The Effect of Stagnant Water Zones

on Retarding Radionuclide Transport in Fractured Rocks: An Extension to the Channel

592

Network Model. Journal of Hydrology 540: 1122-1135.

M AN U

591

593 Tsang, Y. W., and Tsang, C. F., 1989. Flow Channeling in a Single Fracture as a Two

dimensional Strongly Heterogeneous Permeable Medium. Water Resources Research 25, no. 9:

595

2076–80. doi:10.1029/WR025i009p02076.

TE D

594

596 Tsang, C.-F., and Neretnieks, I., 1998. Flow Channeling in Heterogeneous Fractured Rocks. 597

Reviews of Geophysics 36, no. 2: 275–98. doi:10.1029/97RG03319.

EP

598 Tsang, C.-F., Neretnieks, I. and Tsang, Y., 2015. Hydrologic Issues Associated with Nuclear

Waste Repositories. Water Resources Research 51, no. 9: 6923–6972.

600

doi:10.1002/2015WR017641.

AC C

599

601 Tsang, Y. W., Tsang, C. F., Hale, F. V., and Dverstorp, B., 1996. Tracer Transport in a Stochastic 602

Continuum Model of Fractured Media. Water Resources Research 32, no. 10: 3077–92.

603

doi:10.1029/96WR01397.

32

ACCEPTED MANUSCRIPT

604 Xu, C., Fidelibus, C. and Dowd, P.A., 2014. Realistic pipe models for flow modelling in discrete 605

fracture networks. In Proceedings of the First International Discrete Fracture Network

606

Engineering Conference, Vancouver, Canada.

Supplementary material

TE D

M AN U

SC

Visualization of individual particle transport through a channel network (example 3 – section 4.)

EP

-

AC C

608 609

RI PT

607

33

ACCEPTED MANUSCRIPT

Highlights: New scripting library to model flow and transport in unstructured channel networks.

-

An algebraic multigrid preconditioner speeds up the flow solution significantly.

-

All results can be exported for quick 3D visualization in Paraview.

AC C

EP

TE D

M AN U

SC

RI PT

-