Basically, the approach used in this work to generate a pipe network differs from other Discrete Fracture Network (DFN) methods. Usually, DFN approaches drop planar objects (or bonds) in a Euclidean space with a possible conditioning of the objects on various features, e.g., length, density, orientation, etc. Intersections between the objects are in fact the consequence of the ‘bombing’. In the present case, the method for generating a pipe network is inspired from results of previous works at the scale of a single rough-walled fracture plane [16]. It has been shown that flow can be strongly channelled with a major part of water fluxes (say roughly above 75%) passing along very few narrow paths. In addition, it is often observed that some locations are crux nodes (designated as nodes hereafter) through which an important fraction of the flux passes, whatever the general orientation of the head gradient.

The main assumption triggering this work is the existence of such key nodes spatially invariant irrespective of the global gradient direction when flow channelling occurs. In addition, it is assumed, as for several DFN approaches, that channelled flow can be mimicked by a network of 1D pipes connecting these key nodes. To argue in this sense, a synthetic example reported below focuses on a 2D test at the scale of a rough-walled fracture plane.

A fictitious square fracture plane is discretized into a regular percolation network of 90 × 90 bonds (each bond is of size 1 [L]). An aperture sampled in a lognormal distribution (μ, σ in decimal log) is assigned to each bond. ‘Permeameter-type’ boundary conditions are applied at the edges of the regular network (no flow through the east and west limits, constant prescribed head of 1 and 0 [L] for the north and south limits, respectively). The steady-state flow equation is then solved in the fracture plane, and the equivalent macroscopic hydraulic conductivity is computed according to the following expression: K_eq = Q_in/|ΔL ΔH|, with K_eq [LT⁻¹] the equivalent hydraulic conductivity, Q_in [L³T⁻¹] the total inflow (outflow) through the north (south) limit, ΔL [L] the size of the system, ΔH [L] the difference between prescribed heads at upstream and downstream boundaries.

The flow rate through each bond is calculated to locate highly conductive bonds in terms of fraction of the total flow passing through the bond. It is arbitrarily considered that a bond is conductive if more than 1% of the total flow passes through it. A map of the major flowing paths is obtained, assimilated to the network of channels. The latter are kept, whereas the others are eliminated. The resulting network is assumed to correspond to channels in the fracture plane representative of flow for a north–south-oriented general head gradient. Another threshold of 30% of total inflow can be applied over this network for a north–south-oriented general head gradient and then repeated with an east–west gradient. Bonds over 30% of total inflow for both north–south and east–west gradients correspond to invariant key ‘nodes’ through which flow must pass. All the above procedure is repeated for different ‘fields’ of local hydraulic conductivity. The latter follows lognormal distributions of constant mean (μ = −5, in decimal logarithm) and standard deviation σ = 1, 1.5, 2, 2.5, 3, 3.5 and its spatial distribution is either fully random or correlated according to a spherical variogram of one-tenth the size of the system as correlation length. Maps of conductive bonds (above the 1% threshold) for σ = 1, 2 and 3.5 with correlated and uncorrelated conductivity are reported in Fig. 1. Bonds above the 30% threshold (i.e. indicating the location of key nodes) are in bold, whereas the rest, corresponding to flowing channels for a north–south-oriented general head gradient, are in light. A few comments can be drawn from a first look at the maps in Fig. 1a, corresponding to correlated conductivities: an increase in σ makes the most conducting bonds become connected and spatially concentrated, which results in flow channelling. As expected from intuition, this channelling with a few paths holding flow only occurs with correlated conductivities (Fig. 1b). Correlation in space yields clustered values of conductivity that are used as flow paths. Uncorrelated fields generate flowing bonds almost evenly distributed over the domain and non-channelled flows. It must be emphasised that this notion of channelling is empirical here, since it results from an appraisal ‘at first glance’ of the spatial distribution of flowing bonds. An attempt to define indexes of channelling can be done. For instance, invariant key nodes (above the 30% threshold) only appear in the case of correlated fields, which could sign the presence of channelled flow. Other measures are reported in Table 1. When conductivities are correlated in space, a more important part of flow is transferred by a less important part of bonds [23,31] (see the proportion of bonds above the 1% threshold) and the range spanned by local flow rate values is about ten times that from uncorrelated conductivities.

Fig. 2 compares the evolution of the equivalent hydraulic conductivity K_eq with the increase of the standard deviation for correlated and uncorrelated conductivity fields [11,28]. K_eq is quite constant for uncorrelated conductivities, whereas it increases significantly (about three times between σ = 1 and σ = 3.5) when they are spatially correlated. In other words, a network made of a few highly conductive paths would be on average more permeable than a network with the same statistical distribution of local conductivity, but with evenly distributed flow paths [23]. Thus, getting an increased value of macroscopic conductivity compared with a theoretical mean calculated from the statistical distribution could also become an index of channelling. This index is interesting, since it does not need to ‘count’ flowing bonds or to calculate local flow rates. It can be calculated blindly, but requires however a prior guess on the statistical distribution of local conductivity values. Further work is needed to see whether this type of index is relevant and applicable in the field.

Even if this type of evidence is not clearly reported at the scale of 3D fracture networks, one may suggest the existence of such nodes. The method proposed hereafter is first to distribute nodes over space; second to draw pipes between these nodes. Doing this gives the priority to a conditioning on node locations; additional prescribed conditioning on pipe features (length, orientation, etc.) will be used, but with possible bias in the generated network, since node locations are conditioned first. A few additional comments may argue in this sense:

• (1) nodes can be generated randomly but according to prescribed local densities (the number of nodes per unit volume). Of course, this node density is hardly available in field data, but it is expected intuitively, when fracture apertures are spatially correlated, that a densely fractured area will contain more nodes than a scarcely fractured one;
• (2) knowing that each flowing channel is included in a fracture plane, the orientation and the length of the channel are partly (indirectly) correlated with the dip, strike and length of the fracture. Thus, in a synthetic network, the pipe generation may inherit from information on the fracture field when available. However, conditioning is not to be followed strictly, which makes it reasonable that the conditioning of pipe features is partly biased when generating the network, since node locations are imposed first;
• (3) by fixing nodes first and then drawing pipes, the resulting network should be simpler, avoiding for instance dead ends and isolated non-flowing clusters. Further calculations are simplified and the network can be inverted on flow (and/or transport) measurements by making the hydraulic conductivity of the pipes vary, but also by deforming the network. This deformation can be performed by simply moving the node locations. On a numerical standpoint, the graph represented by the network does not change as well as the shape of matrices and vectors used to solve the discrete flow equations; this would obviously speed up the computations. We do not discussed hereafter on inversion capabilities of the networks, but this study is envisioned in the near future based on slug tests and interference pumping tests available from the Hydrogeological Experimental Site (HES) in Poitiers (France) [3,5,12,14].

Given the above general considerations, it is expected from the MODFRAC3D code presented hereafter that it is able to mimic the macroscopic hydraulic behaviour of a 3D fracture network by keeping only the channels where significant flow occurs. It must be pointed out however that even with a ‘strict’ conditioning, the simulated pipe network may be of a different geometry than the assumed real one. Actually, this feature is consistent with the philosophy of homogenisation inspiring this work: local geometry can be partly biased, provided that the global behaviour is respected.

The generation of a network is initiated by means of a Poissonian spatial distribution to locate randomly the nodes, while honouring the constraint of volumetric node density in a cubic domain. The previous example of a fracture plane proved that nodes exist only if the distribution of local conductivities is spatially correlated. Then, for a given threshold on flux, the number of nodes increases with the standard deviation of the conductivity distribution and with the clustering effect separating in space highly conductive from almost impervious sub-areas. The consequence is that node density may vary in space and to allow this, a 3-D network can be divided into entities designated hereafter as sub-networks. A sub-network is spatially located in a parallelepiped volume defined by the min and max coordinates along the x, y and z directions: {x_min, x_max, y_min, y_max, z_min, z_max} and is characterised by its associated volumetric node density (or by a number of nodes).

Once the seeding of nodes in the parallelepiped volume of the sub-network according to the volumetric density is achieved, the next step is to draw pipes connecting these nodes. As stated above, a pipe is entirely located in a fracture plane. Thus, its orientation in space can be conditioned by the strike (direction of the steepest slope of the fracture plane) and dip (angle of this steepest slope with the horizontal plane) of the fracture plane. Moreover, the length, section, and conductivity of a pipe can also be correlated with the fracture aperture [27]. All this conditioning is achieved by associating a pipe with a family. A pipe family is defined by the statistical distribution of strike and dip of the fracture family to which it belongs and by the possible statistical distributions of pipe length, section, and conductivity. These distributions can be of three types: power law, lognormal, or negative exponential. Since pipes are grouped into families, it must be noted that a sub-network is also characterized by the percentages of each pipe-family that it may enclose. A sub-network may contain several ‘isolated’ clusters not connected with each other, but connected to the boundaries of the system (Fig. 3). Note also that isolated clusters of two different subnetworks may be connected between each other if specified in the input file, or if a vertical well is put in the system, establishing hydraulic connection between clusters of different subnetworks located at different depths.

Before describing the cluster generation (Fig. 3), we first focus on the generation of a single pipe. Two possibilities are available for the length of a pipe: either controlled by a power law, a lognormal, a negative exponential distribution or not controlled at all (not statistically constrained).

3.1 Pipe generation with a length obeying a statistical distribution

Knowing the location of a node (node #n), the program seeks another node (node #n + 1) to be connected with the first one by a pipe (Fig. 4). A pipe family is selected randomly according to its chosen proportion (or probability) among the available families. A dip θ and strike ϕ are also drawn from the corresponding statistical distributions. Moreover, a length l is randomly selected in the corresponding length distribution of the family. A ‘virtual’ plane P of dip θ and strike ϕ and passing by node #n is then drawn. On P, a ‘virtual’ circle C is drawn centred on node #n and of radius l. At node #n, the straight line (D) laying on P and perpendicular to the projection of the global head gradient divides C into two semicircles C⁺ and C⁻. They correspond to the positive and negative directions with respect to the global head gradient, respectively. According to a prescribed direction, e.g. along C⁻, the algorithm seeks the node closest to C⁻ and located in a prism of axis (D) and angle α, a specified tolerance parameter. The prism allows us to look for a node #n + 1 yielding a pipe [#n–#n + 1] with a few fluctuations on length strike and dip around the prescribed values. Otherwise, it is very unlikely to find a node #n + 1 giving a pipe able to match up strictly the prescribed values. If a node (#n + 1) is found, the pipe is generated. If no node is found on C, other attempts can be done from node #n after redraw of the pipe length. If the third attempt fails in seeking node #n + 1, the pipe is automatically generated as the line passing by node #n and intersecting the boundary of the system along the direction parallel to the global head gradient. Of course, generating a pipe in this latter case will partly bias the possible prescribed pipe length distribution. But it must be reminded that the priority is to honour first the node locations with the aim to generate a channel network as simple as possible, without dead ends and for which nodes could then be moved for inversion purpose (see Section 2). Note also that the risk to generate by default a pipe intersecting the boundary is lowered when the density of nodes is high and when the search angle α (see above) is well chosen (the weaker the node density, the larger the angle).

3.2 Pipe generation with uncontrolled length

3.3 Cluster generation

Since the generation process is the same for each subnetwork, for the sake of simplicity, only a single sub-network will be considered below. A node is selected randomly among all the nodes belonging to the subnetwork (node #0, Fig. 5). The algorithm consists in generating two ‘half-clusters’ from node #0 and connecting the upper constant head boundary (nodes are sought in the positive direction, see above) and the lower head boundary (nodes are sought in the negative direction). The generation of the first half-cluster begins with the generation of the first bond in the negative direction connected to node #0 (see above for the conditioning–unconditioning of the pipe length). After node #1 (i.e. the tip of the first pipe) has been found, the algorithm is then iterated from node #1 until the half-cluster is connected to the lower constant head limit (nodes #0, #1, #2, Fig. 5). The second half-cluster starting at node #0 is generated in the same way, but in the positive direction (nodes #0, #3, #4, Fig. 5). Finally, the two half-clusters form the entire first cluster. The next cluster is built in the same way. A node is selected randomly among the remaining available free nodes, i.e. the nodes that do not belong to the previous cluster (e.g., node #12, Fig. 5), and a new cluster is generated according to the previous algorithm. However, the generation of a half-cluster stops either when it intersects the boundary of the block (node #14 in the positive direction), or when it intersects a node belonging to a previous cluster (node #5). In the latter case, this yields to the connection of the two clusters. The cluster generation is iterated until there are no longer available free nodes. Note that a node shared between several clusters cannot be connected to more than six nodes, potentially: four pipes representing the four channels of two intersecting fractures and two pipes for the fracture intersection [22].

3.4 Boreholes

3.5 Flow calculation

The code is asked to mimic a fully known ‘synthetic’ network made up of four fractures generated by the code FRACMAN [15] (Fig. 6a). The space is discretized in 100 × 100 × 100 square meshes. Each fracture plane is characterized by a spatially correlated hydraulic conductivity (with a spherical variogram of range the half-length of the fracture plane) following a lognormal distribution (μ = –5, σ = 1 in decimal log). As shown in Section 2, correlation + highly variable values of conductivity are assumed to yield channelled flow with invariant ‘crux nodes’.

‘Permeameter-type’ conditions are applied at the boundary and the flow is then computed for the three possible gradient directions (along the three principal directions of the block). A threshold is then applied as follows: meshes in which flow is higher than 1.5‰ of the total inflow for the three gradient directions are marked, and they will represent the location of the crux nodes for the pipe network generated with MODFRAC3D (black points, Fig. 6b). For a given gradient direction, the aim is to generate a pipe network that best mimics the behaviour of the reference four fractures. The latter are separated into two families (two fractures in each family) of identical dip and strike. The length of the pipes is not controlled, since there are not enough nodes to get a significant representation of any prescribed pipe length distribution. This information as well as the location of the crux nodes is given to MODFRAC3D. As shown in Fig. 6b, the simulator generates the equivalent channelled flow field with a very small number of pipes and one just needs to set the local product KS (conductivity × section) of each pipe at 9.63 × 10⁻³ [L³T⁻¹] to obtain a similar equivalent hydraulic conductivity.

The two networks are then compared in terms of channelling effect by comparing local flux distributions. For the set of four fractures and the given gradient direction, the meshes that have a flow higher than 1.5‰ of the total inflow are kept and assumed to represent the reference channel network. Since ‘crux nodes’ of the pipe network are also determined with the same threshold but for the three gradient directions, they are included in the reference network of channels. Fig. 7 displays the flow distribution for both the reference network of channels extracted from the four fractures and the pipe network generated by MODFRAC3D. The correlation between both distributions is good, the correlation coefficient is 0.97, and the relative errors on mean and variance between both distributions do not exceed 13%. Finally, the pipe network generated by MODFRAC3D has the same equivalent hydraulic conductivity and almost the same flux distribution. One may therefore consider that this pipe network is a good representation of channelled flow in fractures. Of course, the local geometry is not preserved, since the pipe network is a ‘homogenised’ and much simpler image of the real channel network. Nevertheless, this facilitates further modelling exercises and computations.

This work is a first approach to flow and transport in fractured rocks for which the principal conditioning is based on general considerations stating that (i) flow is highly channelled both in each fracture plane and at the scale of the entire fracture network; (ii) there exist some invariant crux nodes through which flow must pass whatever the general flow direction.

The features mentioned above have motivated a new method to build a discrete fracture network. Usually, planar objects (or bonds) are ‘dropped’ in a Euclidean space and intersections of objects are the consequence of the dropping. In the present case, crux nodes are dropped first according to a varying or uniform density in space, and then, bonds are generated between these nodes to build a channel network. Length and orientation of bonds can be conditioned on statistical distributions, but the latter may be slightly biased on the final network, since node locations are honoured first. The bias is a priori all the more important if the node density is small.

In the end, the channel network assumed to represent the major flow paths in a fracture field (or a rough fracture plane) is quite simple in terms of topology; for instance, it does not enclose dead ends and the number of bonds is small. This should facilitate further flow and transport computations. Because of this simplicity, it can be envisioned to inverse flow and transport measurements by making the properties of bonds to vary, but also by deforming the network by moving the node locations. This deformation would not change the topology of the graph representing the network, thus the shape of the matrices of the linear systems to solve the flow (and transport) equations would remain unchanged. Of course, this would speed up the inversion in simplifying computations at each convergence iteration.