Groundwater flow into underground openings in fractured crystalline rocks: an interpretation based on long channels
Groundwater flow into underground openings in fractured crystalline rocks: an interpretation based on long channels
John H. Black 0 1
Nicholas D. Woodman 0
John A. Barker 0
0 Department of Engineering and Environment, University of Southampton , University Road, Southampton SO17 1BJ , UK
1 In Site Hydro Ltd , East Bridgford, Nottingham NG13 8NW, UK 2
Rethinking an old tracer experiment in fractured crystalline rock suggests a concept of groundwater flow in sparse networks of long channels that is supported by results from an innovative lattice network model. The model, HyperConv, can vary the mean length of 'strings' of connected bonds, and the gaps between them, using two independent probability functions. It is found that networks of long channels are able to percolate at lower values of (bond) density than networks of short channels. A general relationship between mean channel length, mean gap length and probability of percolation has been developed which incorporates the wellestablished result for 'classical' lattice network models as a special case. Using parameters appropriate to a 4m diameter drift located 360 m below surface at Stripa Mine Underground Research Laboratory in Sweden, HyperConv is able to reproduce values of apparent positive skin, as observed in the socalled Macropermeability Experiment, but only when mean channel length exceeds 10 m. This implies that such channel systems must cross many fracture intersections without bifurcating. A general relationship in terms of flow dimension is suggested. Some initial investigations using HyperConv show that the commonly observed feature, 'compartmentalization', only occurs when channel density is just above the percolation threshold. Such compartments have been observed at Kamaishi Experimental Mine (Japan) implying a sparse flow network. It is suggested that compartments and skin are observable in the field, indicate sparse channel systems, and could form part of site characterization for deep nuclear waste repositories.
Conceptual models; Crystalline rocks; Hydraulic properties; Numerical modeling; Radioactive waste repositories

Many geoscientists are interested in how groundwater flows
through sparsely fractured crystalline rock, none more so than
those considering constructing nuclear waste repositories in
such formations. During the 1970s and 1980s, several
countries initiated programmes aimed at disposing of radioactive
waste at depth, some of which included the construction of
underground research laboratories (URLs; Witherspoon and
Bodvarsson 2006). One of the major intended activities in
these URLs was to study groundwater flow and radionuclide
transport within poorly transmissive fractures at roughly
repository depth (say ∼500 m below surface). Early tracer
experiments at Stripa Mine, Sweden (Abelin et al. 1985, 1991a),
attempted to inject tracers into the high head gradients induced
by the presence of the mine and then recover them in the
underground opening. This tended to yield low rates of tracer
recovery so subsequent tracer experiments in URLs were
usually located some distance laterally into the rock in artificially
created ‘point sink’ configurations (Sawada et al. 2000; Elert
and Svensson 2001; Andersson et al. 2002; Einsiedl and
Maloszewski 2005). In addition, these borehole to borehole
tests were always sited in major fractures and fracture zones
with concomitantly higher values of transmissivity compared
to surrounding rock: say ∼10−7 relative to ∼10−10 m2/s.
Obviously, these experimental choices were made in order to
ensure positive results since migration rates within low
transmissivity fractures for even conservative tracers over
small distances require many years.
However, one could reasonably question how much results
from fracture and fault zones reveal about flow and transport
in what is often referred to as ‘background rock’ (i.e. the
average rock between fault and fracture zones and referred
to as ‘fracture domains’ by Follin and Hartley 2014). The
origin of such zones ought to produce a morphology and
mineral content that differs from ordinary fractures in
background rock. Indeed, Andersson et al. (2007) concluded
after the extensive series of TRUE tracer tests at spö Hard
Rock Laboratory, Sweden, that uncertainty concerning the
tracer flowpath was greater than uncertainty about the water
−rock interaction. Also that background fractures were about
ten times more retentive than the main fractures involved in
the tracer experiments.
Notwithstanding that the objective of the tracer tests has
been to determine radionuclide transport properties, it is the
premise of this paper that an early, possibly considered
inconclusive, tracer test yielded a valuable insight into how
groundwater flows through background rock. Based on observations
from the aforementioned tracer test, a particular concept of
long channels and sparse intersections is proposed. The
concept is encapsulated within an innovative lattice network code
that is shown to reproduce some of the features noted or
measured in other URL experiments.
The 3D (Migration) Experiment
The ‘ThreeDimensional (Migration) Experiment’ was
undertaken at the Stripa Mine between the end of 1984 and early
1987. The experiment was based on the assumption that a drift
at atmospheric pressure, hundreds of metres below the water
table, would provide a suitable sink for surrounding
groundwater. Then, if different tracers were injected into the local
groundwater system at various locations, they could all be recovered
eventually as inflow to the drift (Abelin et al. 1991a, b).
An experimental drift (roughly 75 m long and 4 m
diameter) was located at least 100 m from existing mine openings
and connected to the old mine galleries by a 140mlong
access drift. Two ‘side drifts’ about 10 m long, were also
excavated about 10 m before the end of the drift, thus forming a
cruciformshaped experimental gallery (Fig. 1). A major
innovation of the experiment was the development of a
polythenesheet catchwater system so that inflow could be
identified as emanating from individual rectangles marked
on the upper walls and ceiling of the experimental drift.
Three hundred and seventyfive polythene sheets, each a
rectangle of 2 m2, were attached to the upper rock surfaces using
resin glue which enabled each area of the grid to yield a
measurement of inflow (Fig. 1a). Inflows varied greatly from sheet
to sheet; less than half of them gave any measurable inflow
and 11 sheets accounted for 50% of inflow. About half of the
gridbased inflow emanated from one of the side drifts. The
results clearly showed that natural inflow is limited to a few
areas with large dry areas between.
For the next stage of the experiment, three 70mlong tracer
injection boreholes were drilled vertically upwards into the
ceiling of the main drift (Fig. 1b). Nine different tracers were
injected in short (2.5 m), discrete, packerisolated zones at
distances between 10 and 56 m directly above the crown of
the drift. They were injected for periods between 60 and
95 weeks followed by 25 weeks of continued monitoring.
Only seven tracers were recovered in the experimental drift,
of which one, bromide, incurred analysis difficulties that
negated attempts to calculate cumulative recovery.
The distribution of tracer inflow was as noteworthy as the
distribution of water inflow. In contrast to the water inflow
that was highest into the NE arm of the crucifixshaped
excavation, all tracers that were recovered at all arrived in the
middle of the experimental drift (Fig. 1b). They ignored the
wettest part of the drift located in one of the ‘arms’ (Fig. 1a).
Also surprising, was the recovery of all four tracers from
borehole C given its proximity to the highinflow ‘arm’.
Perhaps equally surprising was the absence of recovery of
the tracers from the two more distant injection points in
borehole B bearing in mind that they both appeared to show a
pronounced head response to the experimental drift. A final
conundrum was posed by the subsequent discovery, in
significant quantities, of the tracer Eosin Y from borehole A when
excavating a connecting drift about 150 m away from the
injection point.
In summary, the 3D Migration Experiment revealed the
extremely uneven distribution of inflow to underground
excavations together with flow paths that ignore the concepts
underlying groundwater flow in a porous medium. Abelin et al.
(1991b) concluded that the results indicated that much flow
takes place in channels that interact infrequently with
neighbouring channels. Even though Abelin et al. (1991a, b)
conclude that flow occurs via channels within fractures, they
do not attempt to gain further insight by examining the pattern
of mixing exhibited by their experiment. The pattern is
enlightening because most catchwater sheets that yielded tracer,
yielded mixtures from several injection points, and not always
from the nearest injection points. Also, not all possible
mixtures were represented. It seems reasonable to assume that
mixing must occur at certain points where actively flowing
channels meet others. Based on the recovery of the six tracers,
it is possible to reproduce the observed mixtures based on only
eleven ‘mixing nodes’ (Fig. 2). The incidence of mixtures
containing tracers from injection points in boreholes A and
C, many metres apart, suggests that the mixing nodes are
relatively close to the catchwater sheets, possibly within an
excavation damaged zone (EDZ). This in turn suggests that
the channelized flowpaths from the injection points to the
Fig. 1 The 3D (Migration)
Experiment: a natural inflows of
groundwater; b the measured
distribution of tracer inflows from
nine injection points in the three
boreholes A, B and C (after Figs. 5
and 13 in Abelin et al. 1991a)
The sparse channel network concept
The concept of channelbased flow is not new. However, since
Rouleau and Gale (1985) report that fractures at Stripa Mine
are seldom more than 0.5 m apart, the concept of long
channels proposed here means that channels in background rock
cross many fracture intersections without branching.
The concept of sparse long channels with few
intersections agrees with the frequent observation in tunnels and
drifts that inflows are sporadic and wet patches are
separated by large regions of dry tunnel/drift surface (Tsang
and Neretnieks 1998; Tsang et al. 2015). An implication
of isolated points of drift inflow is that groundwater must
converge to those points and that the convergence must
cause head loss proportional to the ‘dimension’ of the
convergence (Barker 1988). Figure 3a shows theoretical
head distributions towards a drift based on a 4 m diameter
drift in a groundwater environment with a 200 m head
Fig. 2 Pseudo threedimensional
representation of the movement
of tracers in the 3D Migration
Experiment invoking imaginary
‘mixing’ and ‘dispersal’ nodes
(constructed from information
contained in Fig. 13 of Abelin
et al. 1991a)
Fig. 3 Head versus distance from the centre of a drift. a Calculated values and log linear regression lines for ideal systems with zero skin. b Results from
packer isolated sections in the radial boreholes of the Macropermeability Experiment at Stripa (based on Wilson et al. 1981)
boundary located 50 m from the drift centreline. Figure 3a
uses a logarithmic xaxis so that, based on porous medium
theory, cylindrically convergent flow should plot as a
straight line. If the convergent flow is of a higher
dimension than cylindrical (i.e. >2D), the line of loglinear
regression projected to the drift wall will intersect it at
positive values of head. Most hydrogeologists would interpret
this configuration as a positive skin, a region of locally
reduced hydraulic conductivity relative to ‘average’ rock
further from the drift (Fig. 3a). The general concept and
mathematics of skin and its application to a drift deep
below the water table is summarized in the Appendix.
Recognizing and estimating skin in fractured crystalline
rocks is difficult because of the inherent variability of
the transmissivity of adjacent fractures. Nevertheless, if
enough points around a drift (actually short straddle
packer intervals in nearby boreholes) are measured, then a
distance–drawdown relationship is usually apparent.
Figure 3b shows such a set of measurements from the
Macropermeability Experiment at Stripa Mine (Wilson
et al. 1981). Similar results have been reported for other
experiments at Stripa, such as the Simulated Drift
Experiment (Black et al. 2006), and for other URLs such
as the KD90 drift at Kamaishi in Japan (Chijimatsu et al.
1996).
The idea proposed here is that observations of isolated
patches of inflow on drift surfaces, positive skin effect and
channelling effects in tracer experiments in ‘background rock’
are all linked to flow in long channels with infrequent
junctions. It was decided to validate the concept by reproducing
the ‘apparent’ skin effect in an appropriate model. The
intention was to construct a numerical model of flow to
limited inflow points to a drift that would, when using parameter
values constrained by the Stripa experiments, yield similar
values of apparent skin effect to those that were observed.
A novel lattice network model
Choice of model type
Studies of flow in fractured crystalline rocks have employed a
wide variety of different conceptual and numerical models
(Selroos et al. 2002). Apart from porous medium and
discrete fracture network models, lattice network models
have sometimes been used to understand flow and solute
transport behaviour in three dimensions. Moreno and
Neretnieks (1993) and Margolin et al. (1998) both used a
cubic lattice model characterized by the simple probability
of any single bond being present. This approach, here termed
the ‘classical lattice model’, has the characteristic that long
channels (continuous lines of bonds) are associated with short
gaps and vice versa.
Since the experiments at Stripa described earlier suggest
long channels separated by large gaps, the classical lattice
model is considered inappropriate. A novel lattice network
concept has therefore been developed that allows the mean
lengths of the channels and gaps to be independent of each
other. This ‘twoparameter bond lattice model’ has been
implemented primarily as a 3D steadystate flow code named
HyperConv. This name derives from the term
‘hyperconvergent flow’ introduced by Black et al. 2006 to represent flow
converging to a point sink within cylindrical boundary
conditions. It was felt that important characteristics of the code
should include the ability to calculate flow and head quickly
in order to run multiple realizations. HyperConv does not
attempt to reflect rock or fracture structure except in as much
as it is threedimensional and channels are organized on an
orthogonal grid.
Hundreds of realizations are generated for two fixed
parameters (e.g. mean channel and gap lengths). For all realizations,
flow through each bond is determined for one of two types of
boundary conditions and the key measures of bulk hydraulic
conductivity and skin effect are evaluated. The measures are
then averaged over all percolating realizations of a particular
ensemble (i.e. of 100 realizations).
The network is formed by creating bonds between the nodes
of a 60 × 60 × 60 orthogonal lattice. A channel is conceived as
a group of bonds in a continuous row. Two probabilities are
used to generate channels along lines of nodes in each
coordinate direction:
1. PN (abbreviation of PNew) is the probability of starting a
new channel if the previous bond (within a particular line
of nodes) is absent.
2. PA (abbreviation of PAgain) is the probability of continuing
a channel to a subsequent node.
This method of network generation gives rise to channel
lengths and gap lengths (in the x, y and zdirections) that both
have a geometric distribution (Table 1). The mean and
variance are given in Table 2. The lengths are all scaled to the
length of a bond (single subchannel) which is given unit
length.
It should be noted that any two of the parameter set {L,
G, PON, PN, PA} are sufficient to characterize the new
lattice model and that the ‘classical lattice model’ is a
special case where PON = PN = PA or, equivalently, 1/L +
1/G = 1.
Applying boundary conditions and assigning parameters
The basic channel network is produced in the form of a cube.
Once this has been generated, different boundary conditions
can be applied, namely facetoface, and cylindrical (Fig. 4).
Facetoface boundary conditions are used to evaluate
percolation thresholds and ‘compartmentalization’, whilst
cylindrical boundaries are used to evaluate skin (in the context of a
tunnel or drift). For facetoface flow, four of the six faces are
set as noflow boundaries confining flow between the
remaining two opposing faces at fixed head. Cylindrical boundary
conditions are achieved by removing surplus nodes in the
region of four parallel edges to form a cylinder and setting
an inner cylindrical boundary at (zero) fixed head (Fig. 4).
The modelling was intended to simulate observations at
Stripa URL so values were assigned to boundary conditions
with the Stripa measurements in mind. The main choice
concerns the ‘equivalent’ distance between adjacent nodes within
the lattice and hence the length of a ‘subchannel’. A value of
1.5 m was chosen so that each 2 m × 1 m catchwater sheet of
the 3D Migration Experiment (Fig. 1) could receive inflow at
maximum model density (PON = 1). This in turn meant that the
external boundary of the cylindrical version of the network
was 44.25 m from the centre of the network. A value of
200 m head was chosen for the external boundary as this
appeared to be the undisturbed background value near to the
3D Migration Experiment (Abelin et al. 1991a). The inner
boundary was a 4m diameter central drift connected to the
atmosphere at 0 m head.
Inflow to the Macropermeability Experiment, in conjunction
with the neardrift hydraulic gradient, yielded an estimated value
of bulk hydraulic conductivity of about 2 × 10−10 m/s. Within
the HyperConv code, channel conductances (Cs) are selected
randomly from a lognormal distribution. All bonds
(subchannels) comprising a channel have the same conductance (Cs).
The bond conductance values are based on a constant base
(geometric mean) value of conductance, Csg. So,
where u is a randomly generated unit normal variable, and σ is
the standard deviation of log10 (Cs).
A value of Csg = 4.5 × 10−10 m3/s was used to yield a target
hydraulic conductivity of 2 × 10−10 m/s for a fully occupied
grid. By setting the standard deviation of the log normal
L = 1/(1 − PA)
G = 1/PN
Distribution of lengths of:
Prfl ¼ ng ¼ PnA−1ð1−PAÞ ¼ ðL−1Þn−1=Ln
Prfg ¼ ng ¼ PNð1−PNÞn−1 ¼ ðG−1Þn−1=Gn
Relationships between probability generators, mean lengths and the bond probability PON
L ¼ 1−1PA
PA ¼ 1− L1
distribution (σ) of log10(Cs) to zero, a constant conductance
network can be produced. In this way it is possible to compare
and therefore distinguish the significance of network structure
on flow from the effect of varying channel conductance.
Solving the flow equation
The last step is to solve the flow equation for the thinned
channel network. The network has mass balance applied at each
node with a linear, Darcian, relationship between flow and head
difference between each bond pair. The resulting linear
equations are solved for the nodal heads using different routines
depending on bond density. Where the network is sparse and
symmetric, a preconditioned conjugate gradient approach is
favoured. For poorly connected systems (less than 30% of the
nodes connected) a general purpose sparse matrix solver with a
partial Cholesky preconditioner is used on a system containing
only the connected nodes. For fuller systems, the full lattice is
solved using a specialized blockstructured preconditioner (see
Black et al. 2006 for more detail).
Finally, the channel network is assessed to identify which
nodes have significant flow in them and which contain
negligible flow. All the nodes containing significant flow (i.e.
values greater than 10−16 m3/s) are termed ‘flowing’ nodes.
The rest are termed ‘connected’ nodes.
The evaluation of skin is based on finding the bestfit
loglinear regression line through the heads calculated at the nodes
Fig. 4 Schematic of the sequence of operations involving the thinning of the HyperConv network during the application of different boundary
conditions and the solution of the flow equation
of the network in the manner identified in Fig. 3b. The
calculation is of ‘skin effect’ (hskin), and has units of metres of water
head. It comes in two versions, hsckin and hfskin that refer to the
calculation based either on all connected nodes (superscript c)
or only on the flowing nodes (superscript f). The values for all
connected nodes, hsckin, is used in the results below because
that is the more likely type of value measured in the field,
albeit at less frequency at greater radial distance.
HyperConv was designed to construct a relatively complex
network, solve the flow equation and derive a value for skin
effect quickly enough that it could be run many times. The
outcomes of many realizations can then be averaged to yield a
single value appropriate to the single set of parameters used to
generate the ensemble (of realizations). A key question
concerned how big the ensemble needed to be. Based on three
runs (σ = 0, 1 and 2), each of about 4,000 realizations, and
using the fact that the variance of a mean is inversely
proportional to the sample size, it is estimated that for 100
realizations the relative error is typically between 2 and 5%.
Therefore an ensemble size of 100 was chosen as a reasonable
compromise between accuracy and computer run time (Black
et al. 2007).
A fundamental property of any network is whether or not it
percolates (i.e. conductive elements connect across the
network). Nonpercolating realizations are generally of little
interest for flow or transport. This study used a series of
ensembles of 100 realizations for five values of mean channel length,
but with different sizes of intervening gaps. Starting with large
gaps (i.e. low density), only a very few realizations within an
ensemble would percolate. As gap size was decreased
increasing numbers of realizations would percolate until all
percolated. The results show that the average number of percolating
realizations is dependent on mean channel length as well as
network density (Fig. 5). More importantly, systems of long
channels continue to percolate at lower values of network
density than short channels. Directly linked to this finding is
that, of the systems examined using HyperConv, classical
lattice networks require the highest values of density to
percolate.
Based on this outcome, the percolation behaviour of a
much larger millionnode cubic lattice designed to evaluate
only connectivity (i.e. without the ability to solve for flow)
was evaluated. The range of behavior of a ‘twoparameter
bond lattice’ type model is shown in Fig. 6 and spans all
Fig. 5 The dependence of percolation probability on channel length and
network density. L = average number of bonds (subchannels) per
channel. The value of probability of percolation of 50% is commonly
referred to as the average percolation threshold or critical density
(Harter 2005)
combinations of average channel and gap length (i.e. all PA
and PN values). The boundary between percolating and
nonpercolating combinations of channel and gap length can be
approximated by the equation
Fig. 6 The combinations of channel and gap lengths that percolate and
do not percolate and their relationship to bond density (PON)
Figure 6 shows the two regions where percolation does and
does not occur, though the boundary is only sharp for infinite
networks. It also indicates that any model with a bond density,
PON, greater than about 0.3 should percolate and that, as PON
falls below this value, percolation only occurs for increasingly
long channels and gaps. Most of the important results from the
HyperConv simulations occupy this low density region. It can
be seen that the possible combinations of average channel and
gap lengths for a classical lattice network are limited to the
single line running diagonally across Fig. 6, corresponding to
1/L + 1/G = 1.
The intention of the modelling was to investigate the effect
of channel length and sparseness on the derivation of values
of skin effect. Margolin et al. (1998) have previously shown
that sparse flow systems can be produced by classical lattice
networks based on variable bond conductance. In order to
avoid confusion with effects caused by variable bond
conductance, all the results reported here (except for some
ensembles investigating variable conductance) concern
systems comprising subchannels that have the same value of
conductance.
The variation with density in Fig. 7 is surprising. A general
trend from extreme values at very low bond density to zero
skin when the lattice is full (i.e. high bond density) might have
been expected. Instead, while the ensemble results of some
shorter channels exhibit a reasonably gradual evolution from
large negative values of skin at low bond density towards zero
at high density, some longer channels exhibit a peak of
positive skin at intermediate values of density.
Fig. 7 The dependence of skin
effect on channel length and
network density. Each point
represents the average value of an
ensemble of 100 realizations or
the number of realizations that
percolated: reported as the
number beside data points. Where
numbers are absent, all 100
realizations of a given ensemble
percolate
Considering the results in more detail, it can be seen that all
networks yielded negative skin effect over at least some range
of density values; that is they indicated a region of apparently
increased permeability in the neartunnel region. However,
only some more dense ensembles of longer channels yielded
a positive skin effect (i.e. apparent reduced permeability near
the tunnel). The division between networks with a region of
positive skin effect and those without, appears to occur around
networks of 6.7 bonds average channel length (10 m). It
should probably be borne in mind that the modelled system
includes a 4m diameter drift as its inner boundary and that
this is of similar size to the 10m scale dividing the two types
of behaviour. Another point of note from the results is that
classical lattice networks show no tendency to yield positive
skin effects in the conditions applied here.
In line with the percolation results shown in Fig. 5, the
expected trend of percolation percentage decreasing with
bond density is maintained with a marked decline for PON less
than 0.04. With decreasing numbers of realizations
percolating, the variability of their average increases as evidenced by
the overlapping data curves at large negative values of skin
effect (Fig. 7).
The relatively sudden transition from negative to
positive skin effect exhibited by networks of longer channels
is a region of particular interest. It was decided to
investigate this transition by examining in more detail the
distance/drawdown behaviour of three ensembles of
10bonds average length per channel (PON ≈ 0.02, 0.04 and
0.07 as identified in Fig. 7).
Head was calculated as a function of distance from the
inner boundary for all connected nodes within every
percolating realization of each of the three ensembles. Values were
placed in annular bins of equal thickness and an average
Fig. 8 a Head versus distance
from the centre of a drift of 4m
diameter for percolating
ensembles based on 15m average
channel length and three values of
density. b Section across the axis
of the inner boundary of the
HyperConv model showing the
configuration of potential inflow
points. Note that the lattice
extends out in all directions
beyond the inner region depicted
derived (Fig. 8a). In all cases, the value of conductance of the
bonds is kept constant (σ = 0) so that all outcomes are a
function of geometry. It should be pointed out that these values of
density represent average ‘gap’ distances between channels
ranging between 125 and 500 bonds (i.e. beyond the
boundaries of the model).
The variation of head with the logarithm of distance from a
drift centreline is plotted in Fig. 8a against a background of
1D, 2D and 3D flow in a porous medium with zero skin. The
straightline sections extending from the outer boundary for
all three cases indicates reasonable adherence to cylindrically
convergent flow essentially imposed by the outer boundary.
However, as the three types of flow system approach the drift
wall, the least dense system diverges from cylindrical
behavior at greatest distance and yields an excellent fit to the model
of Barker (1988) with a flow dimension of 1.6. The most
dense system of the three (i.e. PON ≈ 0.07) adheres to
cylindrically convergent flow until very close to the drift where
almost half of the overall head loss occurs within 2 m of the
drift wall. If the procedure summarized in the Appendix is
applied, the data yields an apparent positive skin effect only
for the densest of the three realizations.
It is useful to consider some practical dimensions of the
lattice network at the point of inflow. Eight bonds (potential
inflow points) emanate from each lattice plane perpendicular
to a 4m diameter drift (Fig 8b), so the probability of one or
more of these bonds being ‘on’ is (1 – PON)8 and from this it
follows that the average distance between perpendicular
planes with inflow points to the tunnel is
For the results in Fig. 8a) this gives distances of: 2.1, 3.7
and 6.8 lattice intervals (3.3, 5.6 and 10.3 m), respectively,
increasing with decreasing density. It is notable that for the
two higher density results in Fig. 8a, these intervals are
reasonably consistent with the radial distance at which flow
ceases to be cylindrical, which is analogous to the flow pattern
for a partially penetrating well.
Although the objective of constructing a sparse long channel
model was to simulate skin effect, it is valuable to examine the
relationship between average channel length, channel density
and hydraulic conductivity. The HyperConv model is capable
of calculating overall network conductivity using different
boundary conditions: facetoface and radial flow to a line sink
(Fig. 4).
The results shown in Fig. 9 are much as expected: as bond
density decreases, hydraulic conductivity decreases. In line
with skin effect, the impact of density is most pronounced
where average channel lengths are shortest. It is noteworthy
that the percentage of percolating realizations drops below
100% at roughly the same value of hydraulic conductivity
(about 2.5 orders of magnitude less than the full lattice value)
for all lengths of channel.
Effect of spatially variable channel conductance on skin
effect and hydraulic conductivity
In the investigations reported in the preceding sections, the
conductance of the bonds of the networks in all the realizations
has the same value by setting the standard deviation of the log
normal distribution of log10 (conductance) to zero (σ = 0).
This means that all the effects observed in the model results
thus far are the product of network geometry alone. Obviously,
natural systems cannot be expected to be composed solely of
constant transmissivity channels, so variance in conductance
was introduced by setting σ greater than zero on the set of
ensembles based around an average channel length of 10
bonds (i.e. average length = 15 m and the same as the
investigation of transitional behaviour above) (Fig. 10).
Fig. 9 The dependence of hydraulic conductivity on channel length and
network density. Each point represents the average conductivity of the
percolating realizations from an ensemble of 100 realizations
The skin effect values (Fig. 10a) are displaced upwards but
by a relatively small amount given that the range of
conductances for one standard deviation differ by a factor of 30. This
suggests that the results for equal conductance channels
(σ = 0) are of reasonably general applicability.
The outcome for hydraulic conductivity can be understood
by reference to the analytical result obtained by Rubin (2003)
for an uncorrelated, spatially stochastic, lognormal
permeability field:
where Kg is the geometric mean of K, m is the geometrical
(network) dimension and σy2 is the variance of ln(K).
At high density, the network dimension, m, can be expected
to be at a maximum of 3 which yields a value for Keff of
roughly 4 × 10−10 m/s via Eq (4) given Kg = 2 × 10−10 m/s.
In other words, flow occurs preferentially via the more
conductive channels. Thus, it must be assumed that, for any given
density, the number of channels containing the bulk of flow is
fewer when conductance varies than when it is held constant.
This may also be the cause of the increase in skin effect at
maximum density because flow continues to converge to the
more conductive channels intersecting the inner boundary.
Based on Eq. (4), when the network dimension (m) = 2,
Keff = Kg which agrees well with Fig. 10b) for PON ≈ 1.
The transition from positive to negative values of skin
effect (at PON ≈ 0.04) also marks the boundary between all
Fig. 10 Results for networks averaging 10 bonds per channel but with
variable conductance: a skin effect versus density, and b hydraulic
conductivity. Numbers beside data points indicate percolating
percentage and apply to both figures (a). A ‘spanning’ channel is one
that extends the entire distance between the inner and outer boundaries of
the model
realizations percolating and decreasing percentages
percolating (Fig. 10). At network densities less than PON ≈ 0.04,
results are best understood in the context of a single individual
channel extending the entire distance between inner and outer
boundaries. The parameters appropriate for a single line of
‘flowing’ nodes (i.e. 29 bonds) are:
−62.27 m (loglinear regression to 29
equally spaced head steps of 7.1 m)
It can be seen in Fig. 10b that once node density decreases
to 0.01 and percolation frequency is at 21%, the value of
conductivity of the equal conductance networks decreases to
the value for a single (‘spanning’) channel intercepting the
inner boundary. For the spatially variable networks where
σ = 0.75, each single spanning channel can vary within the
limits of the log normal distribution and an effective hydraulic
conductivity given by the average conductance can be
expected. For a lognormal distribution the average is Kg exp(σy2/2)
giving Kmin = 2.5 × 10−13 m/s for σ = 0.75 which agrees
reasonably well with Fig. 10b).
Bounding the model behaviour
A ‘semigeneric’ approach
It is possible to derive semigeneric sets of results by linking
the handcalculated minimum results (for a single spanning
line of bonds) to those calculated by the lattice network model
using estimated results (Fig. 11). Some of the lattice model
results in the region where only small percentages of
ensembles percolate (and therefore results are variable) have been
manually smoothed (compare Figs. 7 and 9 with Fig. 11a, b).
The resulting diagrams, Fig. 11a, b are not generic since the
axes are not dimensionless. They are termed ‘semigeneric’
because it is considered that realizations based on different
parameters from the set calculated will exhibit similar relative
behaviour: in particular, the peak values of positive skin just
above the percolation threshold for longer channels.
Fig. 11 The relationship between
a skin effect and b hydraulic
conductivity, and bond (network)
density shown in semigeneric
form to indicate generalized
behaviour. Key applies to both
parts of the figure
The semigeneric form of the skin effect results
(Fig. 11a) indicates a division into two groups, those short
channels that always yield a negative skin effect and the
longer ones that exhibit an interval of positive skin effect.
It is surmised that the value of this division is linked to
the particular size of the inner boundary (4 m diameter)
relative to the average length of the channels within each
set of results (10 m in this case).
The envelope of possible values of hydraulic
conductivity (Fig. 11b) is bounded by the density of a single line
of bonds between the inner and outer boundaries (i.e.
1.88 × 10−4) and the density of a full lattice (i.e. 1). The
upper edge of the envelope is formed by the
density–conductivity relationship of multiples of single lines of
connected channels. One of the curiosities of the results set is
that a single line of channels does not yield the minimum
value of hydraulic conductivity. Instead, the minimum
values of hydraulic conductivity are produced by a small
number of very tortuous flow paths. It is not clear whether
short channel networks (i.e. less than three bonds average
length and including classical lattice networks) exhibit
such behavior because the probability of percolation
declines so rapidly with decreasing bond density that
numerical results were unobtainable in this study.
The skin effect results for head versus distance (Fig. 8)
and the explanation of the results for spatially variable
hydraulic conductivity Fig. 10b) indicate the possible
usefulness of a dimensionbased interpretation. A
conceptual model involving three different dimensions is
suggested.
The channels of the percolating lattice form a network
which has a network dimension, Dnet (‘m’ in Eq. 4) which
decreases from 3 at the highest densities (PON approaching
unity) to 1 when the density is below about 1%. As density
decreases, the number of inflow points to the tunnel decreases,
which increases the flow convergence close to the tunnel wall.
This convergence yields positive values of skin and is here
termed the ‘boundarycondition dimension’, Dbc, close to the
tunnel. The skin factor closely relates to the flow dimension,
Dflow (Fig. 10a).
Relating these three dimensions to two parameters of the
lattice model, PON and L, requires some speculative
assumptions. If it assumed that the network dimension is linearly
related to the number of intersections per channel, B, which,
for the twoparameter lattice network model is:
Then, combining Eq. (5) with relations in Table 2, and
bearing in mind that a 1D system will have zero intersections,
Dnet can be defined in terms of PON and L:
where Fnet is a constant.
Similarly, it is proposed that the boundary condition
dimension is linearly related to the average interval between inflow
planes along the tunnel (Δinflow as defined by Eq. 3). Since, if
each plane of nodes perpendicular to the axis of the inner
boundary has an inflow point to the tunnel and results in 2D
flow, then a value of Δinflow greater than 1 should yield a flow
dimension greater than 2; hence, Eq. (7) imposes cylindrical
flow when the inflow interval is equal to the lattice grid
spacing.
1−ð1−PONÞ
where Fbc is a constant.
A reasonable supposition is that the flow dimension is
limited by both the network dimension and the
boundarycondition dimension, so:
Since the intention is to relate apparent skin effect to flow
dimension, Dflow, the results from HyperConv of skin effect
for varying density, PON, and average channel length, L, are
used for appropriate scaling. It seems reasonable that, based
on Fig. 7, the maximum value of skin effect for L = 20 at
PON = 0.03 should equate to a flow dimension of 3. This
yields values for Fnet of 1.5 and for Fbc of 0.3 and gives the
variation of flow dimension with density and channel length
shown in Fig. 12. A better correspondence could not be
expected given the simplicity of the model, the use of averaged
parameters and functions, and the fact that the results in Fig. 7
are averages over percolating realizations alone.
Although Eqs. (6)–(8) do not include the size of the inner
boundary directly, the outcome involves Eq. (7) which is
based on eight possible inflow points connecting the lattice
to the inner boundary. This means that the channel lengths are
a function of tunnel diameter and the interpretation may have
general application.
The implications of a longchannel concept
The idea of long channels between bifurcations does not, at
first, appear to be a particularly radical concept; however, it
does have many implications that are not immediately
obvious. Firstly, the HyperConv modelling clearly shows that
longchannel systems are able to percolate at lower values of
network density than networks consisting of shorter channels.
Assuming that once water finds the shortest pathway between
regions of higher and lower pressure, it will not seek extra
pathways unless there is a change in boundary conditions or
the original pathways alter due to external processes (e.g.
seismicity, mineral precipitation, solution, etc.). This should leave
most groundwater in background crystalline rock flowing
within a channel network at a density just above the
percolation threshold for an average channel length available within
Fig. 12 Variation of flow dimension with density, PON, and channel
length based on Eqs. (6)–(8)
that particular rock. This should mean that the network
containing flowing groundwater is sparse.
Channel sparseness has its own implications. HyperConv
simulations show how, close to the percolation threshold, flow
within a sparse long channel network converges to pass
through a few key conductors that can be regarded as
‘chokes’. The same choke effect is reported in previous
variable conductance lattice network studies (Margolin et al.
1998) and in a fracture network study by Robinson (1984).
Within such a flow system containing chokes, the bulk of head
loss within the system will occur within the chokes, leaving
the network between the chokes with little head loss as a
consequence. These patches of similar head adjacent to
patches of quite different head are termed ‘compartments’.
These ideas of sparseness, chokes and compartments are
discussed in the following.
The main reason that channel networks such as those
modelled here using the HyperConv lattice network model,
emulate the skin effect behaviour observed at Stripa is their
sparseness. It is difficult to appreciate the impact of sparseness
without a visual cue. Figure 13 shows a single representative
2D slice through two HyperConvlike 3D lattice networks.
Figure 13a is of a classical network at critical density (i.e.
50% of realizations percolate), whilst Fig. 13b is of a long
channel network, also at critical density. The channels in
Fig. 13b are 5× longer than those in Fig. 13a and the bond
density is also 5× less. Neither slice has a direct connection
across the network from facetoface so relies for percolation
on connections perpendicular to the 2D slice.
There is ample field evidence of sparseness in the form of
isolated patches of inflow to drifts at Stripa URL—e.g. the
Buffer Mass Test area of Pusch et al. (1985) and the site
characterization and validation drift of Harding and Black
(1992)—into the Kymmen (Moreno and Neretnieks 1989)
and Bolmen tunnels in Sweden (Neretnieks 2006), and into
tunnels in the Aar and Gotthard massifs of Switzerland
(Masset and Loew 2010). At Pocos de Caldas uranium mine
in Brazil, patches of oxidized pyrite associated with
groundwater flow are observed together with channels of 100 m or
more in length (Chapman et al. 1992).
Chokes and compartments
The model of Margolin et al. (1998) was quite similar to
HyperConv, though being of the classical lattice network type
(see Fig. 6) it lacked the independent spatial correlation of
bonds and gaps that HyperConv is able to simulate. Even so,
Margolin et al. (1998) found that by either reducing the
density of channels or increasing the variation of channel
conductance they could reduce the number of percolating flow paths
in a region to a small number. They also reported that chokes
have a greater effect on overall network conductivity as
network density decreases towards the percolation threshold (see
Fig. 8 of Margolin et al. 1998). Hunt (2005), in his review of
percolation theory applied to hydrogeology, stated that,
Bcritical resistances on the optimal path control the entire field
of potential drops.^ In other words, a sparse channel network
should be associated with regions (i.e. clusters of channels) of
roughly equal head (termed ‘compartments’) separated by
head drops that reflect the head losses through the dominant
chokes in the flowing system.
Fig. 13 2D slices from two HyperConvtype lattice networks at critical
bond density: a a classical network (average channel length (L) = 1.331
bonds, average gap length (G) = 4.02 bonds), in which lines A and B
identify two of the many continuous lines of gaps across the entire
network; b a long channel network (L = 10, G = 200 bonds)
Almost all of the results reported here are based on constant
conductance (σ = 0) networks, where any chokelike
behaviour must result from the convergence of several channels
towards a single channel. Whilst that behaviour appears to
be exhibited by flow in long channels converging to a central
opening (Fig. 8a), it does not demonstrate chokes and
compartments within a more general flow regime.
It was decided therefore to search for chokes using the
facetoface boundary conditions shown in Fig. 4 and networks of
long channels close to critical density. The approach required
the development of an addon code (POOLS) that searches for
clusters of nodes, between which no adjacent pairs have a
head drop greater than a predetermined value (Δh). This
process is repeated a number of times using increasing values of
Δh until there is a small irreducible number of clusters,
sometimes a single cluster.
An initial investigation of this type was performed on three
realizations based on the network generation parameters used
to investigate skin around tunnels. That is, a mean channel
length of 10 (subchannels/bonds) and average subchannel
densities (PON) of 0.02, 0.04 and 0.08 where Pc is just less
than 0.02. The densest network did not yield any specific
chokes and the number of clusters decreased gradually with
increasing Δh. The two networks closer to the percolation
threshold displayed behaviour associated with chokes, and
heads and flows within the denser one (i.e. PON = 0.04) are
shown in Fig. 14. The system contains five clusters of
wideranging sizes, namely 5,301, 418, 289, 71 and 21 bonds. One
of the clusters (289) is entirely separate from the rest of the
network and contains flowpath 1 linking top and base directly.
The largest cluster virtually coincides with the high head
channels (coloured red) located in the top half of Fig. 14a and
connects via several chokes to the remaining three clusters.
The contrast in density between Fig. 14a, b shows how few
channels provide the dominant flow pathways. Given that
Fig. 14 is of a cubic system of 90 m edge, it shows how
unlikely it is for boreholes in reallife tracer tests to intercept
a flowing groundwater channel network.
In summary, this limited investigation suggests that only
networks close to the percolation threshold are likely to
Fig. 14 Heads and flows in a
single realization of a HyperConv
network where L = 10, G = 240
bonds and PON = 0.04: a head
variation [note that the chokes and
flowpath 1 have been visually
enhanced for diagram clarity]; b
Location of flowing channels
(thickness of line is proportional
to the logarithm of flow rate)
exhibit chokeandcompartment behavior. Hydrogeological
compartments (of head) have only recently begun to be
identified in fractured crystalline rocks since it requires
timeconsuming groundwaterbased investigations (Sawada et al.
2000). Notwithstanding this difficulty, two projects around the
world seeking sites for deep repositories for nuclear waste
have begun to describe fractured crystalline rocks in terms
of compartments (e.g. Kamaishi research mine in Japan
(Sawada et al. 2000) and Bátaapáti possible repository site
in Hungary (Benedek and Dankó 2009). Other repository
investigations have recorded abrupt head changes in deep
borehole profiles but possibly only the UK Sellafield investigation
had sufficient density to warrant a compartmentbased
interpretation (Black and Barker 2016).
The investigation at Kamaishi research mine in granitic
rock in Japan was probably the bestreported example of head
compartments (Sawada et al. 2000). Seven subhorizontal
boreholes (totalling 560 m) about 350 m below surface
intersected more than 3,000 fractures of which about 600 were
considered open by borehole TV logging. Some of the
fractures were extrapolated between boreholes based on their
location, similarity of orientation and other properties such as
alteration and mineralization (Fig. 15a). Hydrogeological
testing for a tracer test revealed a complex arrangement of
interborehole connections that did not seem to reflect the preceding
structural analysis; however, crosshole testing did agree with
preexisting knowledge about the disposition of head and a
system of compartments was hypothesized (Fig. 15b).
Sawada et al. (2000) conclude, Bthe cause of the
compartmentalization of flow at this site is not well understood.^ They
also concluded that, Bconductive fractures are heterogeneous
internally, and a conductive fracture may not produce flow at
all points of borehole penetration.^ The key aspect of the
compartment definition at Kamaishi is the form of the
boundary of the high head compartment with its complicated
intrusion into its neighbours. It is obviously impossible to employ
essentially planar features to define compartment boundaries
of such complexity; however, it is exactly how one would
expect to observe the complex interactions between the
clusters of channels of Fig. 14a).
Fig. 15 Plan views of the hydrogeological and tracer testing area at
Kamaishi Mine: a boreholes, open fractures and extrapolated fractures;
b head compartments as derived from crosshole testing and
environmental head measurements. The compartments are coloured
Physical appearance of longchannel networks
There were three phases of tracer testing performed at Stripa
Mine: smallscale double fracture testing (Abelin et al. 1985),
smallscale single fracture testing (Abelin et al. 1994) and the
mediumscale 3D experiment described in the preceding
(Abelin et al. 1991a, b). Channelling was seen to underly the
Fig. 16 Schematic of how the sparse channel network appears at Stripa:
a artist’s impression by the author, b a ‘perspectivized’ version of
Figure 19 of Abelin et al. (1994) with dimensions added; c a section of
borehole wall, a ‘perspectivized’ version of Fig. 9 of Abelin et al. (1994)
from dark equals high head to lightest equals low head. Numbers in
compartments are values of head in metres above the access drift floor
(figure adapted from Figs. 3, 7 and 8 of Sawada et al. 2000)
results from the outset and a detailed picture emerged from the
holetohole measurements in the single fracture experiment
(reproduced here as Fig. 16b, c).
The intertwined braided channels in the artist’s impression in
Fig. 16b are important to the overall long channel concept
because they help to explain certain phenomena and are similar to
the samplebased simulations of Hirano et al. (2010). The
arrangement satisfies the conclusions about deadend channels
containing water that is readily accessible to tracer in actively
flowing channelized groundwater (Neretnieks 2006). It also
accords with the observation of ‘channel swapping’ in the Site
Characterization and Validation Experiment at Stripa (Harding
and Black 1992), since it places ‘semiseparate’ channels of
likely similar flow resistance in close proximity. This means that if a
channel becomes clogged with either precipitate (Milodowski
et al. 1998) or exsolved gas (Long et al. 1992), then there is a
nearby alternative flowpath readily available. Over long periods
of geological time, precipitates of very different geological age
could end up in close proximity within the same fracture such as
was found in the Borrowdale Volcanic Group rocks beneath
West Cumbria, UK (Milodowski et al. 1998).
It should be noted that channelized flow has been proposed
many times before, though mainly at the scale within a single
fracture (Neretnieks et al. 1982; Hakami and Larsson 1996;
Brown et al. 1998 and Vilarrasa et al. 2011). Some authors
have proposed multifracture channel models, though most,
like Cacas et al. (1990), envisage channel bifurcation at every
fracture intersection. Only the model of Tsang et al. (1991)
explicitly conceptualizes channels crossing fracture
intersections without bifurcating. The concept visualized in Fig. 16
only differs from that of Tsang et al. (1991) in the
extensiveness of the nonbifurcating channels such that the distance
between major channel junctions—i.e. ignoring channel
choices within clusters (Fig. 16b)—is shown to be 10–20 m.
A review of a decadesold tracer experiment in Stripa URL
coupled with the widely noted patchiness of inflows to tunnels
prompted the concept of sparse networks of long channels—
that is, networks of channels that do not bifurcate at every
fracturetofracture intersection but form long sinuous
channels with only occasional junctions.
The concept has been investigated using an innovative 3D
twoparameter bond lattice network model named HyperConv
where the mean length of ‘strings’ of connected bonds
(channels) and the gaps between them can be controlled by
separate probability functions. It is found that networks of
long channels are able to percolate at lower values of (bond)
density than networks of short channels. A general
relationship between mean channel length, mean gap length and
probability of percolation has been developed which incorporates
the wellestablished result for classical lattice network models
as a special case.
Positive skin effects were common in the underground
experiments at Stripa and were interpreted as regions of reduced
hydraulic conductivity. HyperConv simulations were
performed in this study using cylindrically convergent boundary
conditions and yielded apparent positive skin under certain
circumstances. More specifically, only realizations of
channels longer than the diameter of the inner boundary (e.g. the
tunnel) exhibited positive skin and reached maximum values
just above critical density. Based on these numerical model
results and calculations of a single channel configuration, a
semigeneric pattern of behaviour is postulated that indicates
that networks of short channels (relative to the size of the
underground opening) are unlikely to yield positive values
of apparent skin. A few ensembles based on variable
conductance bonds (all other realizations used constant conductance)
had only a minor effect on this unexpected behaviour.
A relatively simple analytical relationship in terms of
overall channel density and the frequency of channels intersecting
the inner boundary is developed that largely replicates the
behaviour seen in the numerical experiments. Several
researchers have suggested in the past that flow in sparse
networks will be dominated by a few key connecting channels (or
fractures) termed ‘chokes’. A few numerical experiments were
undertaken where it was found that chokes only become
significant close to the percolation threshold. Values of channel
density greater than about 4× critical density seem to rule out
the creation of recognisable chokes. The occurrence of chokes
is linked to the creation of head ‘compartments’ which,
although the HyperConv simulations represent conditions at
the Stripa URL, appear very similar to actual head
compartments reported at Kamaishi experimental mine in Japan.
It is concluded that, based on this relatively limited
numerical study, the longchannel concept is capable of reproducing
many groundwaterrelated effects observed in underground
research laboratories in fractured crystalline rocks. Networks
of long channels are intimately linked to sparseness which is
itself linked to ‘apparent skin’, chokes and head
compartments. Where these effects are observed would seem to
indicate that the groundwater is flowing in a longchannel network
close to critical density. These links between sparse
longchannel networks and observable features in the field present
the possibility of recognising and characterising such
networks in future site characterization projects for nuclear
repositories in fractured crystalline rock.
Acknowledgements The model ‘HyperConv’ was developed by Peter
Robinson of Quintessa Ltd. This study was supported in part by the
Swedish Nuclear Fuel and Waste Management Company, SKB.
Appendix: the interpretation of (hydraulic) skin
around tunnels based on head measurements
‘Skin’ is a concept used to describe poor performance of
groundwater and oil production boreholes relative to what
would be expected under ‘ideal’ conditions. It is a catchall term
for effects due to permeability variations close to the borehole
wall. Within nuclear repository projects, the concept has been
adapted to apply to tunnels and drifts (Rhén et al. 1997), even
though tunnels are horizontal rather than vertical and
timescales, lengthscales and rocks are significantly different.
In groundwater hydraulics, skin is usually quantified in
terms of a dimensionless ‘skin factor’, sf, additional to any
well function describing head loss involved in flow through
an aquifer towards a sink. Normally, this factor would be
added to the well function evaluated at the well, but to adopt
the convention of the tunnel being at zero head, the skin factor
is added to the head in the formation outside the well:
where the mathematical terms are as follows:
Using the well function representing cylindrically
symmetrical flow in an infinite system, Eq. (9) can be rearranged to
give:
where an effective radius is introduced given by:
reff ¼ rte−s f
Thus, a positive skin factor (equating to a skin of reduced
permeability) corresponds to an effective radius of the tunnel
which is less than the actual radius, and vice versa (see
Fig. 17e, f).
A somewhat more realworld solution than Eq. (9) is
obtained by assuming a horizontal water table (at a constant
head) combined with constant pressure conditions at the
tunnel wall (like an airfilled tunnel) to yield a solution of more
general application. Equations for the head distribution were
derived by El Tani (2003). The full solution is too complex to
reproduce here but a key result is the flux to the tunnel:
ς ¼ rzt −
t
where zt = depth of the tunnel centre below the water table.
When zt/rt is large, the El Tani (2003) solution differs from
the Thiem equation (with zero head at radius 2rt) by less than
1% when zt/rt > 7.2.
As can be seen in Fig. 17a, the system is symmetrical about
the vertical plane intersecting the axis of the tunnel but not
about the horizontal plane, since all the inflow is sourced from
the water table above. As the ratio of depth below the water
table relative to tunnel diameter increases, the head
distribution becomes more radially symmetrical about the tunnel axis.
Cylindrically symmetrical flow at great depth is the
assumption underlying the wellknown approximate solution of
Goodman et al. (1965). The asymmetry of the solution is seen
as divergent heads dependent on direction from the tunnel
centreline (Fig. 17c).
The maximum amount of divergence occurs between
vertically upward and vertically downward directions and can be
generalized by expressing depth and distance in terms of
tunnel diameters (Fig. 17d). For example, a 4mdiameter tunnel
400 m below the water table (i.e. 100 tunnel diameters) would
have upward and downward measurements of head that
differed by 4 m (i.e. 1 tunnel diameter) at 6.1 tunnel diameters
distance from the drift centreline (i.e. 24 m; from Fig. 17d).
It can be concluded that for either a waterfilled tunnel or a
tunnel at atmospheric pressure, groundwater flow can be
approximated as cylindrically convergent flow to a line source
when at distances of more than a few times the tunnel
diameter. This can be represented by the Thiem equation with skin:
r
ln r þ s f
t
Hence, a plot of head versus the log of distance from the
tunnel centreline should yield a straight line passing through
zero at the tunnel wall (i.e. Fig. 17c). The deviation of real
measurements from this ideal gives rise to interpretations of
skin for tunnels.
Most commonly, head variation around tunnels is
measured in boreholes drilled ‘out’ from the tunnel itself so that
measurements are readily expressed in terms of distance from
the tunnel (Fig. 17b). In nuclear repository projects, boreholes
are instrumented with packers to hydraulically isolate
individual sections so that resultant head measurements apply to
specific borehole intervals (e.g. as in Fig. 17b, e, f).
All of the formulae given earlier are based on the
assumption of flow in a homogeneous isotropic porous medium. In
real fractured crystalline rock, measurements usually vary
much more widely than that caused by the asymmetry of the
theoretical solution. This variation is usually interpreted as
resulting from the basic heterogeneity of fractured rock and
the regression line is assumed to be indicative of the average
behaviour of the groundwater system interacting with the
undergroundopening. In many cases, the regression line does
not intercept the tunnel wall at zero head (Fig. 17e, f).
Open Access This article is distributed under the terms of the Creative
C o m m o n s A t t r i b u t i o n 4 . 0 I n t e r n a t i o n a l L i c e n s e ( h t t p : / /
creativecommons.org/licenses/by/4.0/), which permits unrestricted use,
distribution, and reproduction in any medium, provided you give
appropriate credit to the original author(s) and the source, provide a link
to the Creative Commons license, and indicate if changes were made.
Abelin H , Neretnieks I , Tunbrant S , Moreno L ( 1985 ) Final report of the migration in a single fracture: experimental results and evaluation . OECD/NEA, Stripa Project TR 8503, Swedish Nuclear Fuel and Waste Management Co ., Stockholm, 103 , pp
Abelin H , Birgersson L , Gidlund J , Moreno L , Neretnieks I ( 1991a ) A largescale flow and tracer experiment in granite: 1. experimental design and flow distribution . Water Resour Res 27 : 3107  3117 . doi:10.1029/91WR01405
Abelin H , Birgersson L , Moreno L , Widén H , Ågren T , Neretnieks I ( 1991b ) A largescale flow and tracer experiment in granite: 2. results and interpretation . Water Resour Res 27 : 3119  3135 . doi:10.1029/91WR01404
Abelin H , Birgersson L , Widén H , Ågren T , Moreno L , Neretnieks I ( 1994 ) Channeling experiments in crystalline fractured rocks . J Contaminant Hydrol 15 : 129  158 . doi:10.1016/ 0169 7722(94 )90022 1
Andersson P , Byegård J , Winberg A ( 2002 ) TRUE Block Scale Project Final Report 2: tracer tests in the block scale . Technical report TR0214, Swedish Nuclear Fuel and Waste Management Co ., Stockholm
Andersson P , Byegård J , Billaux D , Cvetkovic V , Dershowitz W , Doe T , Hermanson J , Poteri A , Tullborg EL , Winberg A ( 2007 ) TRUE Block Scale Continuation Project final report . TR0214, Technical report, Swedish Nuclear Fuel and Waste Management Co ., Stockholm
Barker JA ( 1988 ) A generalised radial flow model for hydraulic tests in fractured rock . Water Resour Res 24 : 1796  1804 . doi:10.1029 /WR024i010p01796
Benedek K , Dankó G ( 2009 ) Stochastic hydrogeological modelling of fractured rocks: a generic case study in the Mórágy Granite Formation (South Hungary) . Geol Carpath 60 /4: 271  281 . doi:10.2478/v10096 009  0019 y
Black JH , Barker JA ( 2016 ) The puzzle of high heads beneath the West Cumbrian coast , UK: a possible solution . Hydrogeol J . 24 : 439  457 . doi: 10.1007/s10040 015  1340 4
Black JH , Robinson PC , Barker JA ( 2006 ) A preliminary investigation of the concept of hyperconvergence . Research report R0630, Swedish Nuclear Fuel and Waste Management Co ., Stockholm, 66 pp
Black JH , Barker JA , Woodman N ( 2007 ) An investigation of 'sparse channel networks': characteristic behaviours and their causes . Research report R0735, Swedish Nuclear Fuel and Waste Management Co ., Stockholm, 123 pp
Brown S , Caprihan A , Hardy R ( 1998 ) Experimental observation of fluid flow channels in a single fracture . J Geophys Res 103 (B3): 5125  5132 . doi:10.1029/97JB0354
Cacas MC , 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: calibration and validationthe flow model . Water Resour Res 26 /3: 479  489 . doi:10.1029 /WR026i003p00479
Chapman NA , McKinley IG , Penna Franca E , Shea MJ , Smellie JAT ( 1992 ) The Pocos de Caldas project: an introduction and summary of its implications for radioactive waste disposal . J Geochem Explor 45 : 1  24 . doi:10.1016/ 0375 6742(92) 90120  W
Chijimatsu M , Fujita T , Sugita Y , Ishikawa H ( 1996 ) Coupled thermohydromechanical experiment at Kamaishi Mine: initial data around THM experiment area . PNC report 039503 , JAEA, Ibaraki, Japan
Einsiedl F , Maloszewski P ( 2005 ) Tracer tests in fractured rocks with a new fluorescent dye: pyrene 1 , 3 , 6 , 8  tetra sulphonic acid (PTS). Hydrol Sci J 50 ( 3 ): 543  554 . doi:10.1623/hysj.50.3.543.65026
El Tani M ( 2003 ) Circular tunnel in a semiinfinite aquifer . Tunn Undergr Space Technol 18 : 49  55 . doi:10.1016/S08867798(02)00102 5
Elert M , Svensson H ( 2001 ) Evaluation of modelling of the TRUE1 radially converging tests with conservative tracers . The spö Task Force on Modelling of Groundwater Flow and Transport of Solutes, Tasks 4E and 4F . Technical report TR0112, Swedish Nuclear Fuel and Waste Management Co ., Stockholm
Follin S , Hartley L ( 2014 ) Approaches to confirmatory testing of a groundwater flow model for sparsely fractured crystalline rock, exemplified by data from the proposed highlevel nuclear waste repository site at Forsmark , Sweden. Hydrogeol J 22 : 333  349 . doi:10.1007/s10040 013  1079 8
Goodman RE , Moye DG , van Schalkwyk A , Javandel I ( 1965 ) Ground water inflows during tunnel driving . Eng Geol 2 ( 1 ): 39  56
Hakami E , Larsson E ( 1996 ) Aperture measurements and flow experiments on a single natural fracture . Int J Rock Mech Min Sci Geomech Abstr 33 ( 4 ): 395  404 . doi:10.1016/ 0148 9062(95 )00070 4
Harding W , Black J ( 1992 ) Site characterization and validation: inflow to the validation drift . OECD/NEA technical report 9218 of the Stripa Project, Swedish Nuclear Fuel and Waste Management Co ., Stockholm, 84 pp
Harter T ( 2005 ) Finitesize scaling analysis of percolation in threedimensional correlated binary Markov chain random fields . Phys Rev E 72:026120 . doi:10.1103/PhysRevE.72.026120
Hirano N , Ishibashi T , Watanabe N , Okamoto A , Tsuchiya N ( 2010 ) New concept discrete fracture network model simulator, GeoFlow, and three dimensional channeling flow in fracture network . Proceedings World Geothermal Congress 2010 Bali, Indonesia, 25  29 April 2010
Hunt A ( 2005 ) Percolation theory and the future of hydrogeology . Hydrogeol J 13 : 202  205 . doi:10.1007/s10040 004  0405 6
Long JCS , Olsson O , Martel S , Black JH ( 1992 ) Effects of excavation on water inflow to a drift . In: Myer LR, Tsang C F, Cook NGW , Goodman RE (eds) Fractured and jointed rock masses . Proc. of Conf. on Fractured and Jointed Rock Masses , Lake Tahoe , California. Balkema, Rotterdam, The Netherlands, pp 543  549
Margolin G , Berkowitz B , Scher H ( 1998 ) Structure, flow, and generalized conductivity scaling . Water Resour Res 34 : 2103  2121 . doi:10.1029/98WR01648
Masset O , Loew S ( 2010 ) Hydraulic conductivity distribution in crystalline rocks, derived from inflows to tunnels and galleries in the Central Alps , Switzerland. Hydrogeol J 18 ( 4 ): 863  892 . doi:10.1007/s10040 009  10569 1
Milodowski AE , Gillespie MR , Naden J , Fortey NJ , Shepherd TJ , Pearce JM , Metcalfe R ( 1998 ) The petrology and paragenesis of fracture mineralization in the Sellafield area, west Cumbria . Proc Yorkshire Geol Soc 52 ( 2 ): 215  242
Moreno L , Neretnieks I ( 1989 ) Channeling in fractured zones and its potential impact on the transport of radionuclides . Symposium on Scientific Basis for Nuclear Waste Management XII , Berlin, 10  13 October 1988, pp 779  786
Moreno L , Neretnieks I ( 1993 ) Fluid flow and solute transport in a network of channels . J Contaminant Hydrol 14 : 163  192 . doi:10.1016 / 0169 7722(93) 90023  L
Neretnieks I ( 2006 ) Channeling with diffusion into stagnant water and into a matrix in series . Water Resour Res 42 , W11418, 15 pp. doi:10.1029/2005WR004448
Neretnieks I , Eriksen T , Tähtinen P ( 1982 ) Tracer movement in a single fissure in granitic rock: some experimental results and their interpretation . Water Resour Res 18 : 849  863 . doi:10.1029/WR018i004 p00849
Pusch R , Nilsson J , Ramqvist G ( 1985 ) Final report of the Buffer Mass Test, vol I: scope preparative fieldwork and test arrangement . Tech. Report no . 85 11 , Report of the Stripa Project, Swedish Nuclear Fuel and Waste Management Co ., Stockholm
Rhén I , Gustafson G , Wikberg P ( 1997 ) Äspö HRL : geoscientific evaluation 1997/4results from preinvestigations and detailed site characterizationcomparison of predictions and observationshydrogeology, groundwater chemistry and transport of solutes . Technical report no . 97 05 , Swedish Nuclear Fuel and Waste management Co ., Stockholm, 334 pp
Robinson P ( 1984 ) Connectivity, flow and transport in network models of fractured media . PhD Thesis , Oxford
Rouleau A , Gale JE ( 1985 ) Statistical characterization of the fracture system in the Stripa Granite , Sweden. Int J Rock Min Sci Geomech Abstr 22 /6: 353  367 . doi:10.1016/ 0148 9062(85)90001 4
Rubin Y ( 2003 ) Applied stochastic hydrogeology . Oxford University Press, Oxford, 390 pp
Sawada A , Uchida M , Shimo M , Yamamoto H , Takahara H , Doe TW ( 2000 ) Nonsorbing tracer migration experiments in fractured rock at the Kamaishi Mine, northeast Japan . Eng Geol 56 ( 12 ): 75  96 . doi:10.1016/S00137952(99)00135 0
Selroos JO , Walker DD , Strom A , Gylling B , Follin S ( 2002 ) Comparison of alternative modelling approaches for groundwater flow in fractured rock . J Hydrol 257 : 174  188 . doi:10.1016/S0022 1694(01)00551 0
Tsang CF , Neretnieks I ( 1998 ) Flow channeling in heterogeneous fractured rocks . Rev Geophys 36 ( 2 ): 275  298 . doi:10.1029/97RG03319
Tsang CF , Tsang YW , Hale FV ( 1991 ) Tracer transport in fractures: analysis of field data based on a variable aperture channel model . Water Resour Res 27 ( 12 ): 3095  3106 . doi:10.1029/91WR02270
Tsang CF , Neretnieks I , Tsang Y ( 2015 ) Hydrologic issues associated with nuclear waste repositories . Water Resour Res 51 : 6923  6972 . doi:10.1002/2015WR017641
Vilarrasa V , Koyama T , Neretnieks I , Jing L ( 2011 ) Shear induced flow channels in a single rock fracture and their effect on solute transport . Transp Porous Media 87 : 503  523 . doi:10.1007/s11242 010  9698 1
Wilson CR , Long JCS , Galbraith RM , Karasaki K , Endo HK , DuBois AO , McPherson MJ , Ramqvist G ( 1981 ) Geohydrological data from the Macropermeability Experiment at Stripa , Sweden. Report no . LBL12520, Lawrence Berkeley Natl. Lab., Berkeley, CA, 194 pp
Witherspoon PA , Bodvarsson GS ( 2006 ) Geological challenges in radioactive waste isolation, fourth worldwide review , Rep. LBNL59808 , Lawrence Berkeley Natl. Lab., Berkeley, CA