An alternative approach to understanding groundwater flow in sparse channel networks supported by evidence from ‘background’ fractured crystalline rocks

Size and shape of individual flow-features, and not their ‘organization’ in sets of predominant orientation, are the major influences on the ability of groundwater to percolate through sparse channel networks. Measurements in background fractured crystalline rocks proposed for nuclear waste repositories provide useful insight. Flow-features are observed as locations of increased transmissivity during packer or flow testing in boreholes. They are conceived here as channels on fracture surfaces. Findings are based on numerical modelling and a general formula by Barker (2018) for the percolation of two-dimensional (2D) objects in 3D space. Equidimensional shapes are found to be the least efficient at forming percolating networks. As discs are evolved into highly eccentric ellipses, percolation thresholds for number, area and intersection density decrease. At the same time, the percentage of features forming the active flow path declines from about 10% for discs to a few per cent for 50:1 ellipses. Compiling recent field measurements of area density of flow-features reveals low values within a limited range (0.01–0.8 m−1). When this range is combined with practical values of likely channel width, long narrow flow-features are the only reasonable components of a sparse percolating network. Conventional equidimensional discrete fracture networks are considered unlikely. Innovative field investigation and modelling methods based only on hydrogeological measurements are suggested. It is concluded that this consideration of shape supports the approach, broadly termed the ‘long channel’ concept. Barker J.A. (2018) Intersection statistics and percolation criteria for fractures of mixed shapes and sizes. Comput Geosci 112:47–53.


Your article is published under the Creative
Commons Attribution license which allows users to read, copy, distribute and make derivative works, as long as the author of the original work is cited. You may selfarchive this article on your own website, an institutional repository or funder's repository and make it publicly available immediately.

Channel networks and sparseness
Groundwater flow in sparse networks of channels defies easy characterization (i.e. derivation of parameters such as average channel length and width, frequency of intersections, etc.) largely because of the difficulty of meaningful sampling of such sparse systems. In channel network models such as that of Margolin et al. (1998), sparseness is possibly best indicated by the creation of 'chokes', which are links within a network to which many other links converge as networks are thinned. Hence, intersection density and its associated parameter channel density are probable indicators of sparseness.
In fractured rocks, channels are generally conceived as ribbon-like features on fracture surfaces. As a result, the density of channels is readily viewed as the 'surface area' of active channels per unit volume, or 'area density'. Again, in fractured rocks, it is widely accepted that only a small proportion of all observable fractures (containing channels) actually facilitate flow. Thus, the question arises: 'What is the minimum channel density required to ensure a continuous network of interconnected (flowing) channels across a region under consideration?' Some karstic flow systems can clearly be characterized as sparse networks, but their immense scale prohibits easy and spatially frequent experimental observation. Recently, however, Black et al. (2017) proposed a sparse channel network concept for flow in 'background' fractured crystalline rocks. Such low-conductivity networks offer an opportunity for understanding sparse channel networks because of their relatively small scale and the intensive investigations sponsored by the nuclear repository industry often performed from underground research laboratories (URLs).

Nuclear repositories in fractured crystalline rocks
Crystalline rocks are being considered in many countries as suitable candidates to host deep repositories for long-lived nuclear waste. Their main attributes from a safety viewpoint are their strength, which allows the excavation of large caverns, their tectonic stability and their generally low permeability. Their main drawback is their fractured nature since the joints and fractures they contain enable groundwater flows of variable magnitude in between blocks of intact rock of usually very low matrix permeability (say less than 10 −10 m s −1 ). Crystalline rocks also contain faults, deformation zones and major discontinuities, the most extensive of which are often significantly transmissive (say greater than 1 × 10 −8 m s −1 ) and can potentially transport radioactive leachate to the surface. Fortunately, they are usually few and far between, are identifiable by mapping, are treated deterministically (Rhén et al. 2007;Hartley et al. 2012) and are avoided in the design of repository layouts (Vaitennen et al. 2011;Hartley et al. 2012). The intervening rock, referred to by Hartley et al. (2012) as 'sparsely fractured rock' is here termed 'background crystalline rock' and forms the major geosphere barrier around the canisters of waste.

Characterizing and modelling background crystalline rock (BCR)
It has become common practice in the last couple of decades to describe, characterize and calculate likely groundwater behaviour in BCR probabilistically using a structural interpretation and model as a framework. Such models, known as 'discrete fracture network' (DFN) models assume that every fracture they contain is equidimensional (i.e. roughly circular) because this is required by the structural interpretation that is used to form their basic framework.
Recently, however, Black et al. (2017) proposed that head regimes around underground openings in BCR are best explained by strongly convergent flows caused by the sparseness of the networks of channels in which groundwater is able to flow. In other words, the channel systems involved are exhibiting the presence of 'chokes'. In addition, a tracer test performed at Stripa Mine URL in Sweden is best understood in terms of long sinuous channels that cross many fracture-tofracture intersections without bifurcating and have only very occasional interconnections. Such networks are termed 'long channel networks' (LCNs; Black et al. 2017). Such channels are therefore 'longer' than the average distance between adjacent fractures and a practical indication of sparseness is suggested.
In contrast, in DFN models, flow is assumed to take place within the (entire) plane of each fracture with exchange of mass between fractures taking place at (all) their intersections to form a conductive network in three dimensions (Follin et al. 2014). Clearly, the concept of 'long-channel networks' (LCNs) is inconsistent with the DFN concept of flow within networks of equidimensional fractures. Fortunately, this dichotomy can be enlightened using work recently reported by Barker (2018) who has developed a general formula to define the value of critical area density for percolation in threedimensional (3D) networks of plates of varying and mixed shapes and sizes.

Guide to the paper
The objective of this paper is to gain insights into the size, shape and density combinations that enable the sparsest networks to percolate. This objective is achieved by applying the dimensionless percolation density result of Barker (2018) to identify what are the most important parameters affecting percolation at the limit in sparse networks. The results of Barker (2018) are briefly introduced. The impact of size and shape on percolation is then presented in the form of shape parameters versus area density at the percolation threshold. Recognising that area density of flow-features is closely linked to sparseness in BCRs, the relatively rare results of area density reported to date are presented in section 'Measuring the key parameters of a groundwater flow network in BCR' and their origins in field measurements listed. The next step is to assume that groundwater flow within BCR occurs at or just above critical density so that field measurements of density can be regarded as values of critical density. Assembling a number of practical constraints on shape and size results in relatively small ranges of size-shape combinations that could enable groundwater flow in BCR. Finally, in the section 'Discussion', a putative test and modelling strategy is suggested and the implications for the DFN approach are discussed.
A general formula for percolation applicable to a wide range of geometries Barker (2018) considers a fracture model in the form of a mixture of elliptical plates of zero thickness, each of area A and perimeter P at number density, n (centres per unit volume). Consider the dimensionless density (ρ) defined by: where the angular brackets '〈 〉' represent averaging over all fractures and κ is a dimensionless constant. It is found through Monte Carlo simulations that by choosing κ = 0.774, at the percolation threshold density (n c ) this dimensionless density is always close to a single value: for aspect ratios up to 16. The same result was found to apply to any mixture of convex plate shapes and sizes, provided that for each plate, A and P are replaced by the area and perimeter of an ellipse with the same aspect ratio and product AP. Many previously published percolation thresholds (two examples of ellipse investigations plus one concerning rectangles are shown in Fig. 1a), when converted to the same dimensionless density (Eq. 1), were found to lie within the range given in Eq.
(2) (Fig. 1b). For identical plates with the same aspect ratio, α = a/b, Eq. (3) gives reasonably good correspondence with Monte Carlo critical densities up to at least α = 128: Effects of the shape, size and arrangement of network components on network behaviour Barker (2018) terms the components of his networks as plates.
Since fractured rock is the subject here, it might seem obvious to refer to network elements as fractures. However, since fractures in geological structural interpretation as well as DFN modelling are almost always envisaged as being equidimensional, this seems to prejudge the outcome of the investigation described here. The individual elements of the networks described in the following are termed 'flow-features' and are planar ellipses of variable size and aspect ratio.
In order to make the outcomes from Eqs. (1)-(3) relevant to practical applications, the results below are mainly rendered as values of area density in units of m −1 (L 2 /L 3 ).
where A Vc refers to critical area density.
Effect of flow-feature shape The development of a general formula for the dimensionless value of density of mixed shape plate-type objects in 3D space (Eq. 1) enables the impact of flow-feature shape on network behaviour to be readily evaluated. The initial investigation is of ellipse-shaped flow features located and oriented randomly in space ranging from discs to ellipses with aspect ratios up to 100:1 (Fig. 2). The critical area density (A Vc ) given in Fig. 2, is for identical ellipses of area A and perimeter P (see Appendix A for formulae for A and P), and is based on Eq. (1): where ρ c is given by Eq.
(3). The critical intersection density (A Vc ) is given by (for derivation see Appendix A): For a given area, discs have the smallest perimeter so a network of discs requires the highest density in order to percolate: the difference between discs and 100:1 ellipses is a factor of about 12. Figure 2 also includes the density of intersections at percolation threshold and it can be seen that the effect of shape is more profound than in the case of area density: the difference between discs and 100:1 ellipses is a factor of about 24. ; blue is Yi and Tawerghi (2009); light amber is Mourzenko et al. (2005) Actively flowing part of the network The decrease in intersection density with aspect ratio suggests that the proportion of ellipses actively involved in flow also declines with increasing eccentricity. This phenomenon was investigated using ensembles of 50 realisations of a periodic 'box model' containing up to 16,000 same-size ellipses and extrapolating to an infinite system. Two opposing faces were envisaged as upstream and downstream faces and ellipses within the box were then assigned as 'connected only to either upstream or downstream boundary', 'connected between both boundaries (active flow path)', 'connected to active flow path by a single connection (dead ends)' and 'isolated ellipses/clusters)'. Figure 3 demonstrates that as ellipses become more eccentric, a smaller proportion of them enable percolation ( Fig. 3a) and that this proportion is very small (Fig. 3b,c).
Effect of flow-feature size, size distribution and shape Figure 4 presents the critical area density as a function of an 'equivalent disc diameter', D e ¼ ffiffiffiffiffiffiffiffiffiffiffi 4A=π p , which gives a clearer impression of the size of the flow-features involved. The area density is inversely proportional to D e .
It is clear that both size and shape have significant impacts on critical area density; however, examination of the extremes of the x and y axes of Fig. 4 indicates some practical limits. For instance, values of area density greater than 10 m −1 are unlikely to be indicative of flow in BCR and it is difficult to envisage a disc-shaped fracture of 100 m diameter with fullface flow at 500 m depth in BCR.
Entire systems of same-size flow-features are impossible in nature so most interpretations adopt some sort of mixture. Consider a simple mixture of two sizes of disc. It is already apparent from Fig. 4 that larger discs enable percolation at lower values of area density so the results shown in Fig. 5 are much as expected. However, when there is a large size difference, the inclusion of even a small number of larger discs has a significant effect (e.g. less than 0.6% of discs by number with a 16× larger radius reduces critical area density by a factor of 10).
It has become commonplace in applications of DFN models to describe the variation in the size of fractures in the form of a truncated power-law probability density distribution of equidimensional or near equidimensional flat plates usually arranged in sets (Frampton and Cvetkovic 2010;Hartley et al. 2012;Follin et al. 2014). Use of the power law envisages large numbers of very small fractures and a diminishingly small number of very large fractures. If upper and lower limits are put on the radius then the probability that the radius is less than r is given by the cumulative distribution function: Variation with aspect ratio of the critical area density, Eq. (5), and the critical intersection density, Eq. (6), for ellipses of area 1 m 2 (Hartley et al. 2012;Follin et al. 2014) and use a minimum radius appropriate to an exploration borehole (e.g. 0.038 m) and a maximum value of 564 m (radius of a disc with an area of 1 km 2 ). They used values of the exponent, k, between about 2.5 and 3.5 and their impact is indicated in Fig. 6b. Bearing in mind that the addition of less than 1% of large discs can reduce the value of critical density of a mixture by factors greater than ten, the sensitivity of critical density to the choice of powerlaw exponent is evident in Fig. 6b.

Impact of network organization: orientations of flow-features
All of the impacts described above are based on randomly distributed and orientated discs or ellipses in space; in other words, no specific organization. All recent DFN models applied within proposals for radioactive waste repositories are based on structural geological interpretation that assigns fractures to specific orientations, or sets, five in the case of Forsmark, Sweden (Follin et al. 2014). Random spatial distribution within each fracture 'domain' is also assumed.
The impact of network organization was investigated using a box model in which the ensemble average critical area density of ellipses orientated orthogonally (equally in the three directions) was compared to the average for randomly orientated ones, all being same-size networks. The results show that orthogonally orientated flow-features require a slightly higher density (about 8%) than randomly orientated ones to reach critical density (Fig. 7). Organization in sets that are not orthogonal and not random could be expected to yield results between 0 and 8% denser than randomly organized.
Conclusions from evaluation of shape, size and organization on the critical density of networks of flow-features It is clear from the foregoing that the size and shape of flowfeatures included in any network are the most important factors in enabling percolation to occur at the lowest values of area density. Perhaps surprisingly, and excepting strongly anisotropically aligned fractures in sediments, network organization has b) c) a) minimal impact on network connectivity. These conclusions are important because they suggest that groundwater flow in BCRs could be characterized by measuring the size and shape distributions of the flow-features encountered in a field investigation; structure being relatively unimportant.

Measuring the key parameters of a groundwater flow network in BCR
What are the components of a representative network?
Since it is impossible to directly measure every feature enabling flow in BCR, a probabilistic network model is required in order to estimate flow and transport at anything other than the largest scale. In broad terms, such a model requires information describing the shape, size and spatial organization of the individual flow-features together with how many there are per unit volume and how transmissive they are. Viewed in terms of the preceding ellipse-based assessment, parameter values are required to describe area density of the entire network plus size, aspect ratio and transmissivity of the contained flow-features. Since flow-feature organization has little impact, a network comprising randomly distributed and orientated flow-features will suffice as a first approximation.

Reported measurements of area density
The area density (A V ) of flow-features is not a commonly measured parameter in groundwater investigations except where DFN modelling or solute transport modelling is intended. It can be measured either using packer tests with a short straddle interval or detailed flow meter measurements such as the Posiva flow logging (PFL) technique.
A few results are available in the literature, mostly from repository investigation projects, and show some noteworthy trends (Fig. 8a). It is to be expected that as the lower measurement limit of the field test method declines, more inflows are recorded. The results infer a maximum value of area density of less than 1 m −1 and are all within the range 0.008-0.8 m −1 which only spans two orders of magnitude. The hydraulic results (Fig. 8a) are divided between measurements using packers and those using flow logging, and clearly packer tests tend to yield higher values of area density at equivalent measurement limits. This may be due to the packer tests generally being of shorter duration, having a smaller zone of influence and therefore potentially measuring small zones that are not connected to the large-scale flowing network (Follin et al. 2014) (Tables 1 and 2).
Results are also divided into surface boreholes and underground because the higher hydraulic gradients underground enable testing with a lower measurement limit but are often coupled with more uncertain boundary conditions. A further difference between surface and underground is the predominance of sub-horizontal boreholes and therefore the more frequent sampling of vertical fractures underground. The flow log results from Sellafield are reported in two forms: either as the sum of all the observed inflow points or, where they are close together, multiple inflows are treated as single flow-features. It should be noted that they include a variable lower measurement limit. Values of area density derived from core logging and fracture mapping in the same projects as the hydraulic testing are included in Fig. 8b to show how observable fractures, even open ones, exceed the density of flow-features often by at least one order of magnitude.  The conversion of I L to A V is given in Appendix B. T transmissivity, p packer tests, f flow tests, pfl Posiva flow logging method, masl metres above sea level, SCV site characterization and validation experiment. a Denotes where results have been adjusted using the Terzaghi (1965) approach to compensate for the unequal sampling of fractures by boreholes orientated in a limited range of directions b Denotes authors' conversion of I L to A V by a factor of 2× Reported measurements of the size and shape of flow-features/channels Plates, flow-features and channels Barker (2018) describes his network components as plates.
Because this paper concerns features of unknown shape participating in groundwater flow, the plates have been termed 'flow-features' here; however, laboratory and field experiments in fractured rock show that flow and solute transport often occur along flow channels (Figueredo et al. 2016). The term 'channel' is retained in the following sections.

Channel size: extensiveness or 'length'
It is impossible to directly observe the length of individual channels in situ. Existing studies for repository sites use surface mapping to constrain the dimensions of fracture size distributions albeit through the assumption of equidimensional fractures and full-face flow (Follin et al. 2014). This seems inconsistent with the concept of flow within channels and high values of area density of mapped fractures (Fig. 8b) imply little relevance to channel-based flow. However, the nature of BCR situated between observable features such as major faults and fracture zones sets natural limits on the extensiveness of the active flow-features it could contain. With domain sizes in the order of several hundred metres, channel lengths greater than say 100 m seem unlikely. At the other end of the scale, any channel shorter than 0.5 m is likely to be contained on a single fracture surface. Since Black et al. (2017) propose that groundwater flow in BCR occurs in channels that extend across several fracture intersections without bifurcating, a length of 0.5 m seems a reasonable lower limit in the context of fracture frequency in BCR.

Channel width
It is widely acknowledged that groundwater flow occurs within a limited proportion of a fracture's surface (Figueiredo et al. 2016). Bourke (1987) (Fig. 9a). They calculated the areas of active flow and 'open-but-inactive' at three different values of cross-fracture stress and found that the area engaged in flow is in the order of 10%, virtually regardless of stress. Watanabe et al. (2009) also defined the form of the active flow network on the fracture surface and found it to consist of a 'braided channel' (Fig. 9a) including a 'choke' point and a maximum 'width' of about 50 mm. Abelin et al. (1994) conducted a detailed tracer experiment in a 'single' fracture in an underground laboratory and also reported a complex braided channel system (Fig. 9b,c). They concluded that channels occur in clusters ranging in width from 0.05 m up to 1 m with individual channels having widths between 2 or 3 mm and 100 mm. Cvetkovich and Cheng (2008) and Neretnieks (2006) developed credible interpretations of tracer tests by assuming channels are 0.2 m wide.
Ultimately, the main constraint on channel width is the finite strength of the rock in which the fractures occur. Considering background fractured rock at a typical repository depth (i.e. say 200-1,000 m below surface), it seems unreasonable to expect fractures to sustain open channel widths in excess of 0.5 m except possibly in extremely anisotropic stress conditions. At the other end of the size scale, channels of the order of a few millimetres in width appear likely to form part of a braided network (Fig. 9a) where the 'width' of the braided system is the 'hydraulic width'. Therefore, a reasonable range for channel width, likely to be seen in hydraulic testing, is 0.05-0.5 m.

Transmissivity distribution
The same programme of field testing that yields a value for area density will also provide a distribution of values of transmissivity. However, it is apparent that the distribution of transmissivity values as well as that of area density (Fig. 8a) depends on the test method. Not only does packer testing identify more flow-features (Fig. 10a), but the form of the transmissivity distribution is also different (Fig. 10b,c). The reasons for the difference are considered by Hartley et al. (2012) to be that packer tests have a much shorter duration; thereby investigating a small region that is more prone to 'round-packer' leakage and usually has a lower measurement limit (as in Fig. 10). A general conclusion to be drawn from Fig. 10 is that the larger-scale flow system, as measured by the PFL method, has a median value of transmissivity about an order of magnitude larger than the group of tests thought to represent little more than individual fractures.

Constraining the range of size-shape combinations
It is apparent in Fig. 4 that a specific value of critical density can be attained by an infinite set of combinations of shape and size of individual flow-features. However, this range of combinations can be reduced considerably by applying three (reasonable) assumptions in conjunction with field and laboratory measurements. The assumptions are: & Groundwater flow in BCR occurs within a sparse network of flow-features that is at or just above critical density: an assumption that agrees with the views expressed by Chelidze (1982) and Renshaw (1996). Thus, a measurement of area density in situ is likely to yield a value of area density close to critical density. & PFL results are a measure of the area density of the active flow system plus any dead-end clusters that connect to the test borehole. & Detailed packer test results are a measure of the area density of all possible flow-features.
In addition to these theoretical assumptions, the practical limits on channel width of 0.05-0.5 m already assumed are very important. A further limit on possible flow-feature dimensions derives from the intention here to understand flow in sparse networks. If area density is significant; say for instance, three orthogonal cross-cutting fractures per cubic metre (i.e. A V = 3 m −1 ) then the system would be amenable to analysis using a dual-porosity, porous-medium model (depending on the scale of interest). An area density of 10 m −1 is therefore set as a maximum area density for BCR.
Based on these theoretical assumptions and size and density constraints, it is possible to construct relatively limited sets of combinations of same-size flow-features that would just percolate. Figure 11 is based on the results in Fig. 4 limited by the maximum (0.5 m) and minimum (0.05 m) widths of channel considered reasonable in the preceding text. Figure 11a uses the values of critical density appropriate to a network that includes all flow-feature types (e.g. isolated, dead-end and active). These should be used with values of area density derived from detailed hydraulic testing where the zone of investigation is small (i.e. short-interval packer tests). In contrast, the relationships in Fig. 11b) should be used for datasets from detailed flow testing (i.e. PFL measurements) which are expected to include only 'dead-end' and 'active' flow-features. The difference between the density of all flow-features at critical density and that of only dead-end and the active flow path is discussed in section 'Actively flowing part of the network' and shown in Fig. 3a.
The fixed aspect ratio lines in Fig. 11 can be applied to ellipses of distributed sizes by using the equivalent diameter Þ; the fixed minor axis lines can also be generalized but no simple (closed-form) expression can be given for D e . The range of combinations is significantly limited by the constraints of channel width and field results todate and would appear to virtually exclude the possibility of equidimensional (i.e. disc-like) flow-features in flow systems in BCR.

Methodology and precepts
The nature of flow in BCR precludes a deterministic approach since an analyst cannot ever hope to identify and parameterize every transmissive flow-feature. The best that can be hoped for is to construct a numerical model that broadly represents reality and can be run sufficiently often to predict a distribution of outcomes that encompasses field observations. The approach proposed here involves a two-step process involving first a 'percolation' model to broadly define network geometry, followed by an equivalent lattice model to calculate flow. The percolation model is based on four 'descriptors' of its component flow-features (i.e. density, size, shape and organization). It relies on three precepts that follow from Black et al. (2017) and the ellipse investigations described above, namely that: 1. Groundwater flow in BCR occurs at, or just above, critical density (see the preceding section 'Constraining the range of size-shape combinations').

Flow-features (channels) do not necessarily form junctions
at the intersection of two or more transmissive fractures (the essence of the 'long-channel' concept of Black et al. 2017). 3. In fractured rock with at least three joint sets, the difference in the value of critical density between highly organized and randomly organized is no more than 10% (see preceding section 'Impact of network organization: orientations of flow-features').
The effect of precept 1 is that field measurements of area density represent measurements of critical area density that restrict the number of combinations of size and shape that could form an appropriate network (see Fig. 11). The combined effect of precepts 2 and 3 is to make a flow-feature model independent of any form of structural fracture model. A network involving random flow-feature organization is perfectly adequate to simulate realistic outcomes. Thus, the 'fourdescriptor' percolation model becomes a 'three-descriptor' model. Furthermore, in the absence of relevant flow-feature information from structural interpretation, this three-descriptor version has to rely solely on hydrogeological measurements to parameterize flow-feature density, size and shape.
The proposed approach to modelling flow involves a version of the two-parameter, lattice-bond model introduced previously in Black et al. (2017). The model can be run quickly to produce distributions of outcomes but requires the introduction of a transmissivity distribution from field measurements and the retention of density and size equivalence with the underlying percolation model.

Field measurements
The following methods are proposed. Some are well established and some are speculative and require further investigation. Area density (A V ) is readily estimated from measurements of intersection frequency (I L ) in boreholes using either straddle packer testing or the PFL approach (A V ≈ 2I L ). It is suggested that both methods are employed since the difference is indicative, to a limited extent, of flow-feature shape (Fig. 11). The conversion of a measurement of I L to A V is not particularly sensitive to network organization (i.e. ffiffi ffi 3 p I L ≤ A V ≤ 3 I L ), provided that area density is an average of three orthogonal lines of measurement (boreholes/underground scanlines).
The main problem for the long-channel concept, as it is for any four-descriptor model, is the shape and size of the flowfeatures. It is assumed that these are ultimately best described by distributions of shape and size but at this early stage of the concept, averages must suffice. Perhaps fortunately, existing measurements of area density of BCR cover a relatively limited range of about two orders of magnitude (Fig. 11) which, when combined with laboratory and in situ measurements of flow and tracer 'width', yield a relatively limited range of length and width combinations. Of the geometric variables, only flow-feature width can be measured directly. The approach used by Abelin et al. (1994) indicated two possibilities: measuring the inflows/outflows to fracture coincident boreholes and measuring the length of inflows to the lower half of a drift. The observations were probably modified by the proximity of the drift and any associated excavation damage. Although flow-feature length cannot be measured directly, it may be possible to measure it indirectly. The lattice network modelling of Black et al. (2017) indicates that flow convergence towards a drift can produce a positive hydraulic skin but only in a configuration where the channel length is significant relative to the dimension of the sink to which the flow is convergent. This suggests that a measure of the size of the region exhibiting hyper-convergent flow (i.e. flow convergence to a cylindrical sink at a value of flow dimension greater than 2 might reflect the average channel length (Fig. 12a,b). Such measurements would need to be relatively numerous and would also require a rigorous programme of lattice network modelling to substantiate the experimental concept and determine any correlation between size of skin thickness and flow-feature length.
Methods measuring inflow to an underground surface, either as inflow points, or better as lengths of inflow, are particularly useful since, based on the definition of I A in Appendix B, these could provide an indication of average ellipse perimeter, 〈P〉. The approach adopted by Forsmark and Rhén (2005), using 'nappy liners' around the perimeter of a large diameter borehole might be suitable in very low permeability environments but the recently developed infrared method of Neretnieks et al.
(2018) appears more versatile. They used infrared images to locate seepages and measure their flowrates on the surface of an 80-m-long drift at Äspö, URL Sweden (Fig. 12c).
There may be other indirect methods to derive an indication of channel length such as the 'fracture network testing' described by Sutton (1996), in which an attempt was made to measure the frequency of zone-to-zone responses within a multiple packer string in a single deep borehole. However, all these 'indirect' methods would require significant premodelling and development in the light of the values of area and intersection density reported here.
Constructing an appropriate lattice model Rather than construct a DFN type model incorporating high aspect ratio flow-features, the proposal here is to use an equivalent lattice model of the type described in Black et al. (2017). Such a model would require almost no structural knowledge and would not require complex computations to identify intersections; therefore, larger numbers of realizations could be generated to give better statistics. Although less visually realistic than a DFN model, the example in Black et al. (2017) managed to reproduce behaviours such as skin and compartments that are often observed in URLs. That two-parameter bond-lattice model differs from the 'classic' (one parameter) lattice model (e.g. Margolin et al. 1998) in that the average number of linear contiguous bonds (L) (i.e. 'channels') and linear contiguous bond gaps (G) are independent of each other (Fig. 13): in the 'classic' model 1/ L + 1/G = 1. In both models, the channel and gap lengths follow a geometric distribution and the bond probability (p b ) is: Through Monte Carlo simulation, it was found (Black et al. 2017) that at percolation the two (critical) length parameters are approximately related through: (The value 2.27 ensures that the well-known critical probability of p bc ≈ 0.2488 is obtained for the classic model.) It has been found possible to construct a lattice model equivalent to a measurement-based ellipse model at critical density by imposing the conditions: 1. The models have the same number densities (n; channels or ellipses per unit volume). 2. For a given number density, the average number of intersections ('bonds' to other channels; B), per channel or ellipse is the same. 3. Equation (8) should apply when the number density is at its percolation value (n c ).
Examples of the resulting equivalence are shown in Table 3 where parameters are labelled with subscript 'c' to emphasize that the results apply only at the critical density for percolation. Starting with the semi-major axes of the ellipses, a and b, the other parameters in Table 3 are evaluated from left to right as detailed in Appendix C. The quantity B c is the average number of intersections ('bonds') with a single channel, and should not be confused with the number of bonds per lattice node. The last column gives the bond spacing (Δ c ) in the same units as, a and b. The last entry in Table 3 is included to show that the classical lattice model at percolation can be represented by ellipses with a major axis about 2.5 times the bond length and with an aspect ratio of about 3.6. The relevant equations are given in Appendix C.
At densities away from the percolation threshold, the percolation condition, represented by Eq. (9), is lost and a substitute needs to be found which agrees at percolation. Having thus constructed the geometry of the lattice model, based on structural measurements interpreted using the ellipse model, the final step is to assign values of conductance to the channels in the lattice model. This can be performed as, for example, in calibrating a finite difference model against transmissivities from the field, preferably measured using the PFL method.

Implications of ellipse-based investigation for DFN modelling
The key fact gained from the application of the general percolation formula of Barker (2018)-Eqs. (1) and (2)-is that flow-feature shape, in addition to flow-feature size, is a significant factor in determining the critical density of a network of planar flow-features in 3D space. In contrast, the most complex DFN models of BCR to date, those of Forsmark, Sweden (Follin et al. 2014) and Olkiluoto, Finland (Hartley et al. 2012), assume that flow occurs within a network of equidimensional planar fractures (e.g. low aspect ratio rectangles and discs). Hartley and Roberts (2013) justify the use of equidimensional flow-features (their term 'fractures') within their DFN model of Forsmark on the basis of the claim by de Dreuzy et al. (2000) that the percolation threshold is relatively insensitive to the eccentricity of fracture ellipses. Figures 2, 4 and 11 indicate that this is not correct.
The implications of this assumption are considerable because equidimensional flow-features require the highest values of area density to percolate compared to more eccentric shapes. Hence, faced with low values of area density measured in situ (Fig. 8), the other significant factor affecting critical density, size is used to produce percolation at low values of area density. Referring to the proposed repository zone at Forsmark, Follin et al. (2014) claim that fractures around 10-100 m (radius) form the connected network. Bearing in mind that this refers to BCR at 500 m depth, open 20-200-m-diameter fractures with full-face flow seem extremely unlikely.
Overall, regardless of the value of area density measured in situ, a network of non-equidimensional flow-features will have a higher number of flow-features per unit volume than current DFN models. A minor implication arising from the 'box modelling' of Barker (2018) is the lack of impact of network organization. Providing the rock has at least three sets of fractures and they do not imply extreme anisotropy; then entirely random orientations are slightly more efficient at forming a connected network than one configured in sets (Fig. 7).

Conclusions
The 'long-channel concept', introduced by Black et al. (2017) is at odds with recent applications of DFN models to BCR that are based on equidimensional flow-features. This dichotomy has been addressed by assessing the percolation behaviour of networks of high aspect ratio ellipses versus equidimensional planar objects (e.g. discs) using the general result for percolation recently reported by Barker (2018). Applying the formulae shows that the value of critical area density for networks of planar flow-features in 3D space is most sensitive to flowfeature extensiveness and flow-feature shape. In addition, Monte Carlo simulation shows that network organization has a negligible impact on the value of critical density in rocks with at least three sets of roughly orthogonal fractures. It also shows that less than 10% of the flow-features at critical density actually form the active flow path and that this percentage decreases with increasing flow-feature aspect ratio (Fig. 3a).
Assembling reported field measurements of area density of flow-features reveals a relatively limited range of about two orders of magnitude with an upper limit less than 1 m -1 (Fig.  8a). The values for potential repository sites at the lower end of the range require the average diameter of equidimensional flow-features, as used in recent DFN models, to exceed 100 m (Fig. 4). Indeed, this conclusion is confirmed by Follin et al. (2014) who, concerning their DFN representation of BCR at Forsmark, concluded that: Bfractures around 10-100 m (radius) are the ones that typically form the connected network giving inflows in the simulations.P ractical use of the general result of Barker (2018) requires the assumption that flow is occurring within a network that has a density at or just above the critical value for percolation. This assumption appears reasonable in the light of the low values of flow-feature area density reported to date for such rocks. It also means that the value of area density measured in situ is the critical density for the network under investigation and that a limited range of combinations of size and shape of flow-features could percolate at that value. The range is limited by practical considerations, primarily the likely range of widths of the included channels/ellipses (Fig. 11). It is concluded that it is unlikely that equidimensional features form the groundwater flow system in BCR.
The ultimate aim of understanding groundwater flow in BCR is to be able to make reasonable predictions of flow and transport under various boundary conditions over very long times. A two-stage process is proposed involving the construction of a measurement-based, ellipse-type interpretation of percolation followed by a suitably linked, twoparameter lattice model (Black et al. 2017) to calculate flow and possibly transport. The value of the lattice model approach lies in its ability to produce multiple iterations quickly thereby facilitating the required probabilistic method. A testing strategy concentrating on the measurement of area density by two different methods, more comprehensive measurements of channel width and the possible development of size-related testing is suggested.
It is concluded that this investigation, based on ellipses, further establishes the 'long-channel network' concept as a credible conceptual model in understanding groundwater flow in background fractured crystalline rocks. This paper is an attempt, together with Black et al. (2017), to initiate a fresh approach to understanding and modelling flow in fractured rocks. Based almost entirely on hydrogeological measurements, it conforms to the entreaty of Neuman (2005) that Bsite-specific hydrogeologic models of flow and transport in fractured rocks are most robust when based on directly measurable rock flow and transport properties, and less so when based on properties derived indirectly from fracture geometric models.Â cknowledgements The authors would like to thank the editor and reviewers for their assistance in making these complex ideas more accessible than in their original form. 5. The range of values for averaging over three orthogonal directions: ffiffi ffi 3 p I L ≤ A V ≤ 3 I L Appendix C: equivalence of ellipse and lattice models at percolation Consider an identical-ellipse model at any number density, n, where the ellipses have given semi-major axes of lengths a and b or, equivalently, area A and perimeter P. The aim is to find the three equivalent parameters of the lattice model, L, G and Δ. The number density is given by: And the number of intersection per channel or ellipse (the 'bond number') is given by: A third relation is available at percolation when the number density is at its critical value ( n c ): where γ ≈ 2.27; this was given as Eq. (9) of the main text.