# Scaling of fractured rock flow. Proposition of indicators for selection of DFN based flow models

Philippe Davy, Romain Le Goc, Caroline Darcel, Jan-Olof Selroos

► To cite this version:

Philippe Davy, Romain Le Goc, Caroline Darcel, Jan-Olof Selroos. Scaling of fractured rock flow. Proposition of indicators for selection of DFN based flow models. Comptes Rendus. Géoscience, 2023, Geo-hydrological Data & Models [GDM - Géo-hydrologie, données et modèles], 355 (S1), pp.667-690. 10.5802/crgeos.174. insu-03940371

HAL Id: insu-03940371

https://insu.hal.science/insu-03940371

Submitted on 17 Jan 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Comptes Rendus
## Géoscience
### Sciences de la Planète

Philippe Davy, Romain Le Goc, Caroline Darcel and Jan-Olof Selroos

**Scaling of fractured rock flow. Proposition of indicators for selection of DFN based flow models**

Published online: 12 January 2023

https://doi.org/10.5802/crgeos.174

**Part of Special Issue:** Geo-hydrological Data & Models

**Guest editors:** Vazken Andréassian (INRAE, France), Valérie Plagnes (Sorbonne Université, France), Craig Simmons (Flinders University, Australia) and Pierre Ribstein (Sorbonne Université, France)

<figure><img src="image_2.png" /></figure>

This article is licensed under the
CREATIVE COMMONS ATTRIBUTION 4.0 INTERNATIONAL LICENSE.
http://creativecommons.org/licenses/by/4.0/

<figure><img src="image_3.png" /></figure>

Les Comptes Rendus. Géoscience — Sciences de la Planète sont membres du
Centre Mersenne pour l'édition scientifique ouverte
www.centre-mersenne.org
e-ISSN : 1778-7025
# Scaling of fractured rock flow. Proposition of indicators for selection of DFN based flow models

Philippe Davy<sup>①, *</sup>, <sup>a</sup> Romain Le Goc<sup>②, b</sup>, <sup>b</sup> Caroline Darcel<sup>③, b</sup> and Jan-Olof Selroos<sup>③, c, d</sup>

**Abstract.** The objective of the paper is to better understand and quantify the flow structure in fractured rocks from flow logs, and to propose relevant indicators for validating, calibrating or even rejecting hydrogeological models. We first studied what the inflow distribution tells us about the permeability structure from a series of analyses: distribution of transmissivities as a function of depth, proportion of flowing sections as a function of section scale, and scaling of the arithmetically-averaged and geometrically-averaged permeability. We then define three indicators that describe few fundamental characteristics of the flow/permeability, whatever the scale: a percolation scale $l_s$, the way permeability increases with scale above $l_s$, and the variability of permeability. A 4th indicator on the representative elemental volume could in principle be defined but the data show that this volume/scale is beyond the 300 m investigated. We tested a series of numerical models built in three steps: the geo-DFN based on the observed fracture network, the open-DFN which is the part of the geo-DFN where fractures are open, and a transmissivity model applying on each fracture of the open-DFN (Discrete Fracture Network). The analysis of the models showed that the percolation scale is controlled by the open-DFN structure and that the percolation scale can be predicted from a scale analysis of the percolation parameter (basically, the third moment of the fracture size distribution that provides a measure of the network connectivity). The way permeability increases with scale above the percolation threshold is controlled by the transmissivity model and in particular by the dependence of the fracture transmissivity on either the orientation of the fractures via a stress-controlled transmissivity or their size or both. The comparison with data on the first two indicators shows that a model that matches the characteristics of the geo-DFN with an open fraction of 15% as measured adequately fits the data provided that the large fractures remain open and that the fracture transmissivity model is well selected. Most of the other models show unacceptable differences with data but other models or model combinations has still to be explored before rejecting them.

<footnote><sup>a</sup> Univ Rennes, CNRS, Geosciences Rennes, UMR 6118, 263 Avenue General Leclerc, Rennes 35042, France</footnote>

<footnote><sup>b</sup> Itasca Consultants S.A.S., 29 Avenue Joannes Masset, Lyon 69009, France</footnote>

<footnote><sup>c</sup> Swedish Nuclear Fuel and Waste Management Company (SKB), Evenemangsgatan 13, Box 3091, Solna SE-169 03, Sweden</footnote>

<footnote><sup>d</sup> Department of Sustainable Development, Environmental Science and Engineering, Royal Institute of Technology (KTH), Sweden</footnote>

<footnote>E-mails: davy@univ-rennes1.fr (P. Davy), r.legoc@itasca.fr (R. Le Goc), c.darcel@itasca.fr (C. Darcel), jan-olof.selroos@skb.se (J.-O. Selroos)</footnote>

<footnote>* Corresponding author.</footnote>
The third indicator on model variability is still problematic since the natural data show a higher vari-
ability than the models but the open fraction is also much more variable in the data than in the
models.

**Keywords.** Fracture network, Crystalline rocks, Permeability, Scaling, Indicator, Percolation.

*Published online: 12 January 2023*

**1. Introduction**

In naturally fractured geological reservoirs, *fluid-flow structure* is a major question for hydrogeological modeling [Neuman, 2005]. Especially for the crystalline rocks targeted for energy recovery exploitation, or geological storage of industrial wastes like CO₂ or spent nuclear fuel, fluid flow only takes place in a tiny part of the whole rock, through the connected network of fractures between hydraulically active boundaries, while intact rock matrix permeability is considered as negligible. Fluid flow structure is hence a subpart of the fractured system.

While *fractures are ubiquitous* in crystalline rocks, the complexity of their spatial organization precludes the use of classical methods of homogenization to derive hydraulic properties for equivalent continuous media. Fracture systems result from long-standing geological histories and fracturing processes, leading to a complex multiscale structure. Fracture-size scales commonly cover several orders of magnitudes from fractures smaller than micrometers to millimeters to tectonic faults larger than hundreds of meters to kilometers with fracture size distributions adequately modeled by power-law distributions [Bonnet et al., 2001, Davy et al., 2013]. Their capacity to be hydraulically active is dependent on their topological structure (connectivity, intersections), their openness and transmissivity. These factors together contribute to present a very large variability—in space and intensity—of the inflows that can be measured at depth from flow logging [Follin, 2008]. In parallel, usual data acquisition capacities are far from what would be needed to easily reduce modeling uncertainties for typical site conditions. These aspects contribute to the difficulty to define relevant models of the fluid-flow structure for naturally fractured rocks and emphasizes the need to model the discrete nature of the flow before deriving any homogenized properties.

*Discrete Fracture Network (DFN) approaches*
emerged in the late eighties [Cacas et al., 1990a,b,
Long and Witherspoon, 1985, Long and Billaux,

1987]. Since then, they are widely used and devel-
oped, most often for applications that involve hy-
drogeological and hydro-mechanical modeling, e.g.
[Davy et al., 2006, 2018, De Dreuzy et al., 2013, Lei
et al., 2017, Park et al., 2002, Selroos et al., 2022]. As
they realistically reproduce the discrete nature of the
flows—typically spacing between inflows ranging
from a few meters to tens of meters along wellbore—
they are the most suitable for understanding flow
structures in-situ. They also permit aggregation of
various types of data (geological, mechanical, hy-
draulic, etc.) in a multidisciplinary approach where
the DFN is the unifying element, as part of a strategy
to minimize modeling uncertainties [Selroos et al.,
2022].

The major theoretical questions for hydrogeologi-
cal modeling of fractured media, which directly drive
and impact the flow structure modeling, are scale de-
pendence and upscaling or homogenization. The is-
sue of the *scale dependence* of the hydraulic prop-
erties is widely debated. One of the first compila-
tions on that subject is the seminal paper of Clauser
[1992] which compiled permeability measurements
from lab to regional scales, emphasizing a permeabil-
ity increase with scale increase up to a hundred me-
ters. This conclusion has been questioned by Hunt
[2003b], who posited that although experiments of-
ten indicate an increase in the hydraulic conduc-
tivity with increasing scale, this effect may be ex-
plained by sampling issues rather than intrinsic scale
effect. Other studies, based on field observation [Ill-
man, 2006, Neuman, 2003, Neuman and Di Federico,
2003, Ren et al., 2021] or modeling [de Dreuzy et al.,
2001b, 2002] are in line with the observations from
[Clauser, 1992] and a scale effect emerging from the
fracture network connectivity structure and trans-
missivity distribution. Illman [2006] also observed a
directional permeability scale effect in cross-hole test
results and hypothesizes that the effect is controlled
by the connectivity of fluid conducting fractures. A
scale increase of effective permeability with scale is
also shown in Martinez-Landa and Carrera [2005]
and Ren et al. [2021] from single hole and cross-hole
pumping tests data and identify a scale dependence
likely site-dependent too. de Dreuzy et al. [2001b,
2002] use a DFN model (in 2D) for which the fracture
length and fracture aperture distributions are power
laws and study the resulting hydraulic properties as
the flow structure and equivalent permeability. They
determine the equivalent permeability scale effects,
characterize the flow structure by a channeling in-
dicator, and establish the relationship with the frac-
ture distributions parameters. Upscaling permeabil-
ity is an active research area in hydrogeology both
for fractured and unfractured media. Upscaling con-
sists of deriving the equivalent permeability from the
smaller scale permeability distribution [de Dreuzy
et al., 2010]. Wen and Gómez-Hernández [1996] and
Renard and Marsily [1997] reviewed the state of the
art on upscaling conductivities in heterogeneous me-
dia. Oda [1985] developed an analytical approach to
define equivalent permeability from a DFN descrip-
tion based on geometric characteristics of fractures
with arbitrary orientations. de Dreuzy et al. [2001b,
2002] develop the relationship between equivalent
permeability and the parameters of multiscale DFN
models in 2D and hence the scale dependency. Chen
et al. [2015, 2018] develop numerical approaches to
compute the equivalent permeability from a DFN de-
scription.

The question investigated in this paper is how to
best use available hydro data for capturing the multi-
scale nature of in-situ flow structures and to derive
manageable and relevant metrics for site modeling
purposes. At the Forsmark area in Sweden, site in-
vestigations have been performed by Svensk Kärn-
bränslehantering AB (SKB), the Swedish Nuclear Fuel
and Waste Management Company, for more than two
decades in view of the future deep repository for
spent nuclear fuel [Follin, 2008, Follin et al., 2007, Sel-
roos et al., 2022]. Follin et al. [2014] use flow logs and
metrics based on the relative proportion of sealed,
open, and flowing fracture frequencies together with
equivalent transmissivity distributions to calibrate
the parameters of DFN models for flow simulations.
Maillot et al. [2016] compare the performance of dif-
ferent DFN models from the comparison of chan-
neling and equivalent permeability indicators. Zou
and Cvetkovic [2020, 2021] also use DFN models
and varying hydraulic boundary conditions to sim-
ulate steady-state pumping tests and derive indica-
tors of cumulative and of complementary cumulative

distributions to evaluate the impact of fracture in-
plane heterogeneities on the flow structure at the
DFN scale.

The first objective of this study is to develop flow-
based metrics adapted to in-situ flow testing capaci-
ties and suited to emphasize the multiscale nature of
the flow structure. To address these issues, the data-
base of fractures and discrete inflows at Forsmark,
which includes flow logging and steady-state pump-
ing tests, was analyzed with metrics revisited from
the former work of Maillot et al. [2016]. The second
objective is to test the capacity of the latter to dis-
criminate between DFN models otherwise calibrated
to similar geological and geometrical observations.
A modeling environment is set up in the software
“DFN.lab” to generate DFN realizations from a pre-
defined range of DFN models and perform numeri-
cal flow simulations in steady-state pumping config-
uration simulations.

The paper is organized as follows: the data and
conditions of the Forsmark site are recalled in Sec-
tion 2. The framework to define the metrics is pre-
sented in Section 3. Data analyses are presented in
Section 4. In Section 5 the DFN models are defined
and finally, several a priori equiprobable DFN mod-
els are compared with each other and with the field
data and the interpretations and comparisons be-
tween data and models are presented. Results and
outcomes of the study are finally discussed in Sec-
tion 6 and summarized in the conclusion. To facil-
itate ease of reading this study, specific parts have
been moved to appendices.

Note that, while the analytical tools and method-
ology are generic, the results are specific to the stud-
ied area (mostly the Nordic granites). It is not pre-
sumed that thresholds and trends are generalizable
to all crystalline rocks [e.g., Dewandel et al., 2006].

## 2. Site and available data overview

The Forsmark site in Sweden, located in crys-
talline bedrock in south-eastern Sweden and lo-
cated roughly 130 km north of Stockholm, is selected
by the Swedish regulatory authorities to build a nu-
clear waste repository for spent nuclear fuel at 470 m
depth [SKB, 2011]. Svensk Kärnbränslehantering AB
(SKB) is in charge of site investigations, future con-
struction, and long-term safety assessment of the
final repository. The database built over the course
of the site investigations encompasses data from ge- ological, mechanical, hydrogeological, geochemical, etc., aspects. Several tens of boreholes each reach- ing lengths up to about 1 km have been drilled and exhaustively logged.

The presented analyses relies on a detailed de- scription of the fracture transmissivity structure in the different boreholes that have been extensively carried out at the Forsmark site by SKB with the Po- siva Flow Logs (PFL). PFL was originally developed by Posiva Oy to meet the demand on flow measuring techniques adapted to sparsely fractured rocks and low-permeability environments [Öhberg and Rouhi- ainen, 2000]. The method couples a very low thresh- old for flow detection (down to 0.1 l/min) and a pre- cise positioning of the device that allows them to de- tect single fractures inflow with a depth resolution of a couple of centimeters. The hydraulic properties of every cored borehole at the Forsmark site are hence tested [Rouhiainen and Pöllänen, 2003, Rouhiainen et al., 2004] by performing a standard hydraulic test where the applied pressure head and resulting flow are recorded. The derivation of a transmissivity is based on the Thiem (or Dupuits) solution for steady- state radial flow without skin effect, so that the flow rate Q₁ at a borehole test section of size L is equal to:

$$Q_1 = K \frac{2\pi L}{\ln\left(\frac{R}{r_0}\right)} (h_0 - h_1) \quad (1)$$

with

Q₁: measured flow rate in the test section (m³/s)

K: hydraulic conductivity of the test section (m/s)

L: length of the test section limited by packers (m)

R: radius of influence (500 m) (m)

r₀: radius of the borehole (m)

h₀: undisturbed hydraulic head far from test sec- tion (m)

h₁: pump-induced or natural borehole hydraulic head (m).

The two unknowns K and h₀ can be deduced from two pairs of hydraulic heads (h₁ and h₂) and resulting flow rates (Q₁ and Q₂) provided that K is not depen- dent on the induced pressure (no hydromechanical effects).

A combination of sequential and overlapping flow logging allows identifying water conductive fractures down to a spatial resolution of 0.1 m and flow rates down to 30 ml·h⁻¹. This leads to a measurement limit

of approximately T = 10⁻¹⁰ m²/s when expressed in equivalent transmissivity as defined below.

The transmissivity of a test section, or the fracture identified by the sequential mode, is defined as:

$$T = KL = \frac{\ln\left(\frac{R}{r_0}\right)}{2\pi} \cdot \frac{(Q_2 - Q_1)}{(h_2 - h_1)} \quad (2)$$

with R ~ 500 m it comes that T ~ ΔQ/Δh. The ratio ΔQ/Δh is also called a specific capacity (in (m²/s)). Martin and Follin [2011] rightly note that: “The PFL transmissivity value (T in (2)) reported to SICADA (the SKB database) will be close to the local value at the borehole (i.e., the transmissivity of the intersect- ing flowing fracture), if the local value is less than the overall transmissivity of the network to which it is connected. If the local value at the borehole is greater than the overall transmissivity of the network to which it is connected, the PFL transmissivity value represents a “hydraulic choke” (bottleneck) phenom- enon, which means that it is not the transmissivity of the intersecting flowing fracture that is determined”. The PFL transmissivity T (Equation (2)) is not the ac- tual fracture transmissivity, but rather the transmis- sivity of the closely connected structure. Other limi- tations of (1) and (2) are that R in (1) and (2) varies with time [Perrochet, 2005] and, in case of several fractures, T or K tends towards the geometric mean [Meier et al., 1999]. Thus T must be taken as an esti- mate rather than an exact value of the fracture trans- missivity, which is correct considering the high vari- ability of transmissivities (i.e., the distribution of the estimate gives a faithful measure of the transmissiv- ity distribution).

In practice, the same set of flow rate measure- ments is done twice all along a tested borehole: the first pass is for natural flow rate conditions (no pumping), and the second one is for a typically im- posed drawdown of Δh = 10 m.

The distribution of PFL transmissivities for all wells is given in Figure 1b. The log average is about 3 × 10⁻⁸ m²·s⁻¹ with values ranging from 2.5 × 10⁻¹⁰ to 1.3 × 10⁻³ m²·s⁻¹. The average intensity of PFL log-transmissivities does not vary significantly with depth (Figure 1c) in contrast to the density of PFL fractures, i.e., the number of fractures where a trans- missivity has been detected, and to the transmissiv- ity variability (Figure 1a). From these observations, we identify three main domains: the near surface (0–200 m), the domain of interest for the repository
<figure><img src="image_4.png" /></figure>

**Figure 1.** (a) (top left) Number of PFL fractures (i.e., fractures whose transmissivity is detectable with the PFL flow log) per 20 m depth increment. (b) (top left), Distribution of the decimal logarithm of PFL transmissivities in m²·s⁻¹; the bars are for the whole domain (100–1000 m) and the solid blue dots for the target domain between 220 and 520 m. No value smaller than 10⁻⁹.⁵ are observed; the impact of this threshold, although very low, is discussed in Section 4.3 and in Figure 6. (c) (bottom left), PFL transmissivity as a function of depth for all the boreholes. (d) (bottom right), Example of a borehole logging (KFM02A) showing the measured density of fully intersecting fractures (number per borehole length with sealed and open densities indicated by the solid and yellow dot, respectively) on the left and the measured PFL transmissivity on the right.

(200–500 m), and a deeper domain (>500 m). The near surface is highly fractured and highly permeable; the density of PFL fractures is 1 PFL every 5 m on average (Figure 1a). It is characterized by horizontal sheet joints, which makes them distinct from the deeper parts. The intermediate domain is rather homogeneous with a density of PFL fractures about half the near surface. The deep domain is very dry with only a few PFL fractures, on average 1 every 2T50 m.

The spatial distribution of transmissive fractures varies from one borehole to another but with

common features as clustering in flow zones, presence of large domains without detectable flowing structures, and positive correlation between zones of high fracture density and zones of large PFL transmissivities (Figure 1d). Note that the correlation between fracture density and PFL transmissivity is far from perfect and many counterexamples exist where PFL fractures are observed corresponding to zones of normal fracture density.

In the following, we focus on the domain of interest for the waste repository. We arbitrarily set the upper and lower depths to 220 m and 520 m,
respectively, in accordance with the density of the PFL fractures observed in Figure 1a. The distribution of PFL transmissivity between 220 and 520 m is not different from that observed on the whole domain (100–1000 m), despite the variability in the density of PFL fractures (Figure 1b).

### 3. Scaling metrics for flow log interpretation

In what follows, an interpretation framework based on the PFL data, including the transmissivity values and their spatial distribution is introduced. This is done through the prism of scale analysis, where the in-situ scale-dependent permeability is defined and evaluated for domains of increasing size.

The primary component of the analysis is the total inflow through a given core log section:

$$K_i = \frac{1}{L_i} \sum_{j=1}^{n_i} T_{i,j} \qquad (3)$$

where $K_i$ is the equivalent hydraulic conductivity (in (m/s)) of the core section i of length $L_i$, into which there are $n_i$ positive inflows interpreted as transmissivities noted $T_{i,j}$, from (2). The core section length $L_i$ will be used to analyze the scale dependency of $K_i$ in order to highlight fundamental properties of the site permeability, and in particular, the scaling laws that govern it. The scale analysis of hydraulic conductivity relies on the continuous records of all the discrete inflows (the PFL transmissivities as recalled in Section 2) over an entire borehole. The maximum investigated scale of the analysis is the length of the borehole in the investigated domain. For a given scale L (i.e., all section length $L_i = L$ in (3)) and a given borehole, the arithmetic average of $K_i(L)$ is exactly the hydraulic conductivity of the borehole if there is no overlap between the different sections, or a very close value, if there is some. In Appendices A and B, we argue from numerical simulations that this borehole conductivity is representative of the permeability of a 3D system with a typical size equal to the section length. Hereafter, we use the term permeability to refer to $K_i$.

To reveal the flow network organization, we focus on the distribution of “wet” sections (i.e., sections containing at least one inflow), as it is done when determining the fractal properties of a network [Tél et al., 1989]. Whether a section is dry or wet depends on the resolution of the inflow measurement. Even if

<figure><img src="image_5.png" /><figcaption>**Figure 2.** Schematic evolution of the geometric average of the permeability K as a function of the measurement scale L (solid blue line). The thin red segments indicate the three main regimes: below the percolation threshold scale $p_c$ (I), between $p_c$ and the representative elementary volume (REV) (II), and above the REV (III) [Charlaix et al., 1987, de Dreuzy et al., 2001b, 2002]. See text for an explanation of the three regimes.</figcaption></figure>

the PFL instrument is extremely sensitive, it is clear that the transmissivity distribution is not complete and truncated at the PFL resolution (Figure 1b). This issue will be discussed further down in the text.

By taking the geometric average of the section permeability, as in de Dreuzy et al. [2002], we expect to identify some important characteristics of the flow structure (Figure 2): the percolation threshold scale $p_c$ (transition between regime I and II in Figure 2), the permeability increase above $p_c$ (regime II), and possibly a Representative Elementary Volume (REV) that indicates the scale above which the permeability no longer varies. The correspondence between the percolation threshold and a scale was established for fracture networks with a power-law size distribution [Berkowitz et al., 2000, Bour and Davy, 1997, 1998, Darcel et al., 2003b, de Dreuzy et al., 2000]. Compared to networks of constant size, the percolation threshold reflects the connectivity of small fractures and the probability of encountering fractures larger than the size of the systems that connect the networks on their own. It is this dual connectivity that makes connectivity scale-dependent and gives a correspondence between the percolation threshold and a given scale. Below the percolation threshold, the wet sections are likely containing one main fracture/channel
with a permeability that decreases as $\overline{K} = T_o/L$ where $T_o$ is the channel transmissivity and $L$ the section length. The permeability increase above the threshold has been demonstrated by [Charlaix et al., 1987, de Dreuzy et al., 2002]. It reflects the distribution of channel transmissivities. Above $p_c$, the scale increase results in a merging of different channels with different transmissivities, and the permeability increase with scale reflects the broadness of the transmissivity distribution [de Dreuzy et al., 2001b, 2002]. Depending on the transmissivity distribution, it may exist a REV over which systems can be handled as homogeneous media. The REV is a correlation length of the transmissivity distribution but its existence is questioned if the distribution of fracture transmissivities and lengths is too broad [de Dreuzy et al., 2001b].

To calculate $K(L)$, we subsample each borehole with section lengths smaller than the total borehole length obtaining a set of equivalent hydraulic conductivities $K_j$ (Equation (3)). We note $N_{>0}(L)$ the number of core sub-sections containing at least one positive inflow and hence a positive $K_j$. The first proposed metric, the function $K_a(L)$, is the arithmetic average of the above-defined equivalent permeabilities derived from the ensemble of flowing core sections (non-flowing core sections are discarded) for a size scale $L$, as summarized in the equation below:

$$K_a(L) = \frac{1}{N_{>0}(L)} \sum_{i=1}^{N_{>0}(L)} K_i(L). \quad (4)$$

Since the permeability arithmetic average including zero-flow sections $\overline{K}_a$ is independent of $L$ (see section following (3)), $K_a(L)$ varies as the inverse of the percentage of flowing sections:

$$K_a(L) = \frac{N(L)}{N_{>0}(L)} \overline{K}_a. \quad (5)$$

The second proposed metric, the function $K_g(L)$, is the geometric average of the same set of elementary bricks:

$$K_g(L) = \exp \left( \frac{1}{N_{>0}} \sum_{i=1}^{N_{>0}} \log(K_{(i)}(L)) \right). \quad (6)$$

This metrics has been proposed by de Dreuzy et al. [2002] as the most representative mean for a permeability distribution that tends to be log-normally distributed (Figure 1).

Both $K_a(L)$ and $K_g(L)$ are likely to be dependent on the inflow detection limit estimated at $10^{-10} \text{ m}^2 \cdot \text{s}^{-1}$ for PFL measurements. Even if this

value is very low, the transmissivity distribution in Figure 1b suggests that smaller values exist. By artificially increasing the detection limit, we test in the following section that the metrics do not depend too much on its value as long as it remains much in the lower end of the transmissivity distribution.

To characterize the flow channeling, we also calculate the average number of inflows in all flowing sections of length $L$, $n(L)$, and a channeling metric $n_Q(L)$ derived from [Maillot et al., 2016]:

$$n_{Q,i} = \frac{\left( \sum_{j=1}^{n_i} T_{ij} \right)^2}{\sum_{j=1}^{n_i} T_{ij}^2} \quad (7)$$

where $i$ refers to the core section $i$ that contains $n_i$ inflows noted $T_{ij}$, over the section size $L_i$. The quantity $n_Q(L)$ is then the arithmetic average over all flowing core sections of size $L$, and it can be viewed as the number of efficient channels in a flowing core section. If all inflows have the same intensity, $n_Q$ is equal to the number of inflow $n$, however, in the general case with a distribution of inflow intensity, $n_Q$ is smaller than $n$. The ratios $n/L$ and $n_Q/L$ have the same meaning as the number of fracture intersections per unit core length, classically referenced as $p_{10}$, but only for flowing (n)—or significantly flowing (n_Q)—fractures.

# 4. Hydraulic data interpretations

In this section, we compute the different metrics and indicators defined in the previous section for the datasets listed in Table 1 in the depth range 220–520 m (see Section 2 for an explanation of this depth range).

## 4.1. Permeability averages

The scale evolution of equivalent arithmetic and geometric permeabilities are plotted in Figure 3. From borehole to borehole and comparable size scale, the calculated values are spread over 4 orders of magnitude. Despite this large dispersion, all the plots display a similar shape with a $L^{-1}$ decrease over about 1 to 2 orders of magnitude when increasing scales until a stable (for $K_a$) or increasing (for $K_g$) regime is reached with large-scale permeabilities that range between $2 \times 10^{-10}$ and $4 \times 10^{-7} \text{ m} \cdot \text{s}^{-1}$. The transition between both regimes is at about 30 m for both permeability averages.
<figure><img src="image_6.png" /><figcaption>**Figure 3.** $K_a(L)$ (left) and $K_g(L)$ (right) computed for each borehole dataset listed in Table 1 (grey symbols) as a function of scale. For both, the yellow dots indicate the average of all boreholes and the red dots, the permeability average (i.e., $K_a$ for the $K_g$ graph and vice versa). The blue lines represent power-law fits at small and large scales, respectively. The dark yellow dots in the $K_g$ plot are for KFM08A, a borehole that will be discussed in the modeling section.</figcaption></figure>

<table><caption>**Table 1.** List of boreholes (IDCODES in the SICADA database) with PFL data including the total number of inflows, the depth of the uppermost and lowermost inflow, the smallest and the largest PFL transmissivity</caption><thead><tr><th>IDCODE</th><th>Number of inflows (T > 0)</th><th>Uppermost inflow (m)</th><th>Lowest inflow (m)</th><th>Min(T) (m²/s)</th><th>Max(T) (m²/s)</th><th>Total logged length (m)</th></tr></thead><tbody><tr><td>KFM01A</td><td>34</td><td>105.3</td><td>363.4</td><td>2.45 × 10<sup>-10</sup></td><td>5.31 × 10<sup>-8</sup></td><td>891</td></tr><tr><td>KFM01D</td><td>34</td><td>106</td><td>571.17</td><td>6.59 × 10<sup>-10</sup></td><td>2.30 × 10<sup>-6</sup></td><td>708</td></tr><tr><td>KFM02A</td><td>125</td><td>101.8</td><td>894</td><td>6.16 × 10<sup>-10</sup></td><td>4.21 × 10<sup>-5</sup></td><td>899</td></tr><tr><td>KFM03A</td><td>71</td><td>106.4</td><td>994</td><td>8.90 × 10<sup>-10</sup></td><td>9.21 × 10<sup>-5</sup></td><td>896</td></tr><tr><td>KFM04A</td><td>142</td><td>109.6</td><td>954.8</td><td>7.16 × 10<sup>-10</sup></td><td>3.54 × 10<sup>-5</sup></td><td>876</td></tr><tr><td>KFM05A</td><td>27</td><td>108.9</td><td>720</td><td>4.45 × 10<sup>-10</sup></td><td>1.23 × 10<sup>-3</sup></td><td>897</td></tr><tr><td>KFM06A</td><td>99</td><td>102.4</td><td>770.8</td><td>2.40 × 10<sup>-10</sup></td><td>1.92 × 10<sup>-5</sup></td><td>895</td></tr><tr><td>KFM07A</td><td>23</td><td>110.8</td><td>261.4</td><td>3.57 × 10<sup>-9</sup></td><td>7.67 × 10<sup>-5</sup></td><td>892</td></tr><tr><td>KFM07C</td><td>15</td><td>98.4</td><td>279.8</td><td>8.99 × 10<sup>-10</sup></td><td>4.81 × 10<sup>-5</sup></td><td>413</td></tr><tr><td>KFM08A</td><td>41</td><td>107.6</td><td>687</td><td>2.48 × 10<sup>-10</sup></td><td>2.20 × 10<sup>-6</sup></td><td>846</td></tr><tr><td>KFM08C</td><td>21</td><td>102.4</td><td>683.6</td><td>7.51 × 10<sup>-10</sup></td><td>2.95 × 10<sup>-6</sup></td><td>848</td></tr><tr><td>KFM10A</td><td>56</td><td>60.3</td><td>484.4</td><td>1.64 × 10<sup>-9</sup></td><td>2.79 × 10<sup>-5</sup></td><td>437</td></tr></tbody></table>

The evolution of $K_a(L)$ and $K_g(L)$ reflects the organization of flows along the boreholes. In the first regime, as long as the length of the core sections is smaller than the minimum distance between two inflows, every element contains only one inflow and hence the equivalent permeability is inversely proportional to the core section length:

$$K_j = \frac{T_j}{L}. \qquad (8)$$

As expected, the arithmetic average of $K_j$ is larger than the geometric average but the scaling trends are similar for this small-scale regime. We note that, for $K_a$, the average of all boreholes (yellow dots) seems to decrease slightly less than $L^{-1}$, $\propto L^{-0.8}$. This reflects a difference between the least and most permeable boreholes; the latter (KFM02A, KFM03A, KFM04A, KFM10A) having a less steep decrease with scale than the former. A basic explanation is that a high
<figure><img src="image_7.png" /><figcaption>**Figure 4.** Plot of the average number of inflows per section (left) and of the channeling metrics defined in (7) (right) as a function of section scale.</figcaption></figure>

permeability is also associated with a higher density of inflows. This is not observed for the borehole-averaged $K_g$ because the geometric average gives a more distributed weight to all boreholes whatever their permeability.

When the section length is large enough to gather several inflows, the permeability stabilizes with scale for the arithmetic average $K_a$, or grows as ~$L^{0.7}$ for the geometric average $K_g$ as the result of the transmissivity distribution. The transition between both regimes spans over one order of scale magnitude between ~6 and ~60 m with a critical value at 30 m determined from the intersection of both end-member regimes (Figure 3). There is no limit to the increase of $K_g$ with scale in this regime for the observed scale range (30–300 m), which means that a REV such as anticipated in Figure 2 is beyond 300 m, if it exists.

The average number of inflows per section confirms the trends observed on the permeability averages (Figure 4, left). For sections smaller than 4–5 m, only one inflow is observed on average per wet section, and then the number increases as $L^{0.75}$. The maximum number of inflows for the largest section of 300 m varies between 10 and 80 with an average of 20, except for the two boreholes KFM01A and KFM07C where only 1 and 2 inflows are recorded, respectively. We also calculate the channeling metrics $n_Q$ (Equation (7)), which quantifies the number of main channels in each section. $n_Q$ shows the same two-regime scaling trend as the number of inflows per section but with smaller values, as expected. It increases above 4–5 m approximately as the square root of the number of inflows. For the largest sec-

<figure><img src="image_8.png" /><figcaption>**Figure 5.** Evolution of the wetting ratio, i.e., the proportion of “wet” sections, i.e., containing at least one inflow, with section length scale $L$ for all the boreholes (grey symbols) and for the average over all the boreholes (yellow dots). The scaling is $L$ (dashed line) and $L^{0.3}$ (solid blue line).</figcaption></figure>

tion (300 m), only 3 main channels on average are detected.

## 4.2. Flow structure organization

We plot in Figure 5 the evolution with $L$ of the ratio of the proportion of “wet” sections—i.e., number of sections that have at least one inflow divided by the total number of sections—with scale. This ratio, hereafter called the *wetting ratio* $r_w$, controls the evolution of $K_a(L)$ (see discussion above and (5))
<figure><img src="image_9.png" /><figcaption>**Figure 6.** Evolution of the permeability averages $K_a$ (left) and $K_g$ (right) as a function of section scales for different transmissivity thresholds (see text).</figcaption></figure>

and it can be used to identify a potential fractal nature of the flowing structure as it is classically done in box-counting methods [Mandelbrot, 1982, Pavón-Domínguez and Moreno-Pulido, 2022]. In a fractal of dimension $D_{1D}$ in 1D, the number of wet sections decreases as $(L/H)^{D_{1D}}$ with $H$ the investigated borehole length (here 300 m), and the wetting ratio increases as $r_w \sim (L/H)^{1-D_{1D}}$. $D_{1D} = 0$ is a point-like structure in 1D corresponding to a plane in 3D. Figure 5 shows the two regimes already described: for scales smaller than ~30 m, the wetting ratio increases as $L$ corresponding to individual flow planes. Above 30 m, $r_w = (L/H)^{0.3}$ indicating a fractal structure with a 1D dimension of 0.7 (2.7 in 3D). This complex structure and the distribution of transmissivity are likely to be responsible for the increase in the geometric average of permeability.

### 4.3. Dependency on the detection threshold of transmissivity

The previous analysis is potentially dependent on the flow detection resolution to identify wet and dry (no-flow detected) sections. Both $K_a$ and $K_g$ rely on the permeability of “wet” sections defined as the dual of no-flow sections. If inflows smaller than the detection limit of transmissivity exist, the number of wet sections will increase and the average permeability decreases. Even if the detection limit is very small, it does not seem to represent the lower bound of the transmissivity distribution (Figure 1b).

Considering that we cannot decrease the measurement resolution, we first apprise the reader that

the previous scaling analysis is performed with a transmissivity threshold of $10^{-10} \text{ m}^2 \cdot \text{s}^{-1}$ and that the results may depend on this threshold. Then, we performed the same analysis as above but with different detection thresholds of transmissivity (Figure 6). The curves are identical for detection threshold ($T_d$) up to $10^{-10} \text{ m}^2 \cdot \text{s}^{-1}$ and similar, i.e., presenting the same log-shape but with a higher permeability for $T_d$ up to $10^{-8} \text{ m}^2 \cdot \text{s}^{-1}$, which is the mode of the probability distribution of log-transmissivity. We thus conclude that the detection threshold is likely to be small enough to allow for a relevant scaling analysis of the flow structure, but we will be certain of this conclusion only if analyses were carried out with an even lower detection limit.

### 4.4. Flow structure indicators

From the previous analysis, we identify a series of flow structure indicators which could serve as a basis for testing, and calibration of validating hydrogeological models. The indicators are best chosen to be characteristics of the fractured media that are likely to be controlling the flow properties, that can be measured, and whose measure is not too much affected by the conditions of the experiment. The 4 indicators listed below and in Figure 7 are chosen to be the measures of transmissivity distribution and flow structure.

For each borehole, the equivalent geometric permeability average (6) with scale $L$ can be characterized by a V shape (Figure 7) with a scattering from one borehole to another. These trends are basic to the definition of our 4 indicators.
<figure><img src="image_10.png" /><figcaption>**Figure 7.** Representation of the 4 indicators (number in circles) described in the text with a graph which superposes the data of Figure 3 and the schematic scaling evolution of Figure 2. The first is the percolation threshold scale, the second is the permeability increase above the percolation threshold, the third is the distribution of transmissivities, and the fourth is the REV scale.</figcaption></figure>

The first indicator is the scale at which network percolates, which is about 30 m in the studied area. This scale also represents the average distance between inflows or inflow clusters. The evolution of $K_g(L)$ below the percolation is representative of sections in which the number of inflows does not increase with scale, which explains the $L^{-1}$ scaling.

The second indicator is the increase of the geometric permeability average with scale above the percolation threshold. It reflects the spatial distribution of inflows and the probability to encounter a flow path of large intensity when increasing the section scale. The slope of this increase is characteristic of the flow organization, a steeper slope indicating a higher heterogeneity of flows. Note that the increase of permeability with scale is long debated in the hydrogeological literature [Hunt, 2003a,b, Illman, 2006, Neuman, 2005, 1994, 2003, Neuman and Di Federico, 2003].

The third indicator is the dispersion of the borehole curves around the $L^{-1}$ decrease. From borehole to borehole and comparable size scale, the calculated values are spread over 3 orders of magnitude (Figure 7). At sizes $L$ smaller than the minimum distance between two inflows $d_{min}$, $K_g(L)$ is simply equal to

the geometric average of the set of individual transmissivities divided by $L$:

$$K_g(L)_{L<d_{min}} = \frac{1}{L} \times \left( \prod_{i=1}^{N} T_i \right)^{\frac{1}{N}} = \frac{\overline{T_g}}{L} \quad (9)$$

A borehole of total investigated length $H$ is thus characterized by a set of $N$ transmissivities $T_i$ drawn from the full transmissivity distribution of the whole site with a geometric average $\overline{T_g}$. Both $N$ and the distribution of $T_i$—and thus $\overline{T_g}$—are indicators of the full transmissivity distribution. Note that the distribution of $\overline{T_g}$ must vary with the investigated borehole length $H$ around the geometric average of $\overline{T_g}$ (yellow dots in Figure 3).

In theory, it exists a fourth indicator which is the representative elementary volume (REV), above which the permeability becomes scale independent. Figure 3 shows that the REV scale, if it exists, is larger than the investigated section length of 300 m.

## 5. Calibration/validation of DFN models

This section aims at testing the relevance of the indicators to calibrate, validate or reject DFN models. It should be taken as an illustration of the potential of the indicators and not as an in-depth analysis of the modeling.

The DFN modeling process presented here is in line with the DFN methodology defined in Selroos et al. [2022], where the definition of a DFN model for a hydrogeological application involves several steps. In the first, geological and geometrical data are used to define the DFN model of all the fractures, whatever their internal properties may be (by convenience hereafter referred to as geo-DFN). The second one is an in-between between geometric and hydraulic modeling steps. It consists in defining the DFN model of all open and partly open fractures, as a subpart of the geo-DFN model (by convenience hereafter referred to as open-DFN). Before even assigning transmissivity properties to the fractures of the open-DFN model, and hence defining the hydro-DFN model, the structure of the open-DFN model is critical for the flow: having a connected path of open fractures between hydraulically active boundaries is the prerequisite to define hydraulic properties. In each step of the modeling process, data are integrated with prior models and hypotheses and evaluated in a cycle
<table><caption>**Table 2.** List of the different fracture sets for the FFM01 fracture domain at depths between –200 and –400 m with their intensity *P*<sub>32</sub> and the Fisher parameters that define their orientation distributions</caption><thead><tr><th rowspan="2">Set name</th><th rowspan="2">P<sub>32</sub> (m<sup>2</sup>/m<sup>3</sup>)</th><th colspan="3">Fisher distribution parameters</th></tr><tr><th>Trend (°)</th><th>Plunge (°)</th><th>Kappa (-)</th></tr></thead><tbody><tr><td>NS</td><td>1.292</td><td>292</td><td>1</td><td>17.8</td></tr><tr><td>NE</td><td>1.733</td><td>326</td><td>2</td><td>14.3</td></tr><tr><td>NW</td><td>0.948</td><td>60</td><td>6</td><td>12.9</td></tr><tr><td>EW</td><td>0.169</td><td>15</td><td>2</td><td>14.0</td></tr><tr><td>HZ</td><td>0.624</td><td>5</td><td>86</td><td>15.2</td></tr></tbody></table>

*P*<sub>32</sub> values are taken from Tables 6–33 in Glamheden et al. [2010] and Fisher parameters from Table C-1 in Follin [2008].

of Sensitivity Analyses, Calibration, and Rejection tests Selroos et al. [2022].

In this section, we present only a small part of this modeling process which focuses on the integration of hydrological data analyzed in Section 4. The other parts are only briefly described. The analysis focuses on the fracture domain FFM01 [Follin et al., 2014, Olofsson et al., 2007], which is the target geological formation for the spent nuclear fuel repository [SKB, 2011]. Models are compared with the borehole KFM08A (Figure 3b, dark yellow dots) that is representative of FFM01.

## 5.1. *Description of the selected models*

The geo-DFNs rely on site investigations and data interpretations at the Forsmark site [Darcel et al., 2009, Fox et al., 2007, Olofsson et al., 2007]. The fracture density and fracture orientation distributions are measured from borehole logging and core mapping. The number of fractures fully intersecting the boreholes per unit core length *P*<sub>10</sub> is converted into total surface of fracture per unit volume *P*<sub>32</sub> by means of stereological rules [Dershowitz and Herda, 1992, Darcel et al., 2003a, Davy et al., 2006, Piggott, 1997]. The values of *P*<sub>32</sub> for the different orientation poles recorded in Forsmark are given in Table 2 with references therein.

The multiscale organization of fracture networks is a critical attribute for fracture connectivity and flow [Bonnet et al., 2001, Bour and Davy, 1998, Davy et al., 2006]. It results in fractal properties and power law fracture size distributions, which are confirmed by the analysis of outcrop and lineament trace maps [Darcel et al., 2009]. Davy et al. [2010] have shown

that the outcrop data in Forsmark are consistent with a kinematic model of fracture nucleation, growth and arrest (hereafter called the UFM model). The resulting fracture-size distribution is a double power law with a scaling exponent close to –4 above a critical fracture size *l*<sub>c</sub> and –3 below. The geo-DFN relies on the genetic UFM model with rules given in Davy et al. [2013]; it is characterized by the critical fracture size *l*<sub>c</sub>, which fixes the fracture density *P*<sub>32</sub>; the larger is *l*<sub>c</sub>, the smaller is *P*<sub>32</sub>.

The open-DFN is a subset of the geo-DFN where only open fractures—i.e., fractures with a non-nil open aperture—are represented. The only data available is the surface ratio of open fractures (hereafter referred as the open fraction *f*<sub>op</sub>), which is varying between 0.15 and 0.25 in the geological formation [Follin et al., 2014, Olofsson et al., 2007] with a value around 0.15 for the target fracture domain FFM01 [Doolaeghe Wehowsky, 2021]. *f*<sub>op</sub> can be measured from the core logging by identifying open and partially-open fractures and measuring the ratio of fractures presenting open voids, but this is a bulk ratio that does not indicate how the fractures are opened individually or depending on their characteristics. We suspect that the relationship between the open fraction and the fracture attributes (size or orientation) is a key parameter to define the hydraulic properties of the open-DFN.

We analyze different models for the distribution of the open fraction in the network, whose parameters are given in Table 3. In all the models, we assume that fractures are either completely open or sealed. The first two models—further referred as “*l*<sub>c</sub>-open”—are built as genetic UFM models [Davy et al., 2013]
<figure><img src="image_11.png" /><figcaption>**Figure 8.** Density distribution (left) and open fraction (right) as a function of the fracture size for the Geo-DFN, *l*<sub>c</sub>-open, α-open and *pp*-open models, respectively.</figcaption></figure>

<table><caption>**Table 3.** Parameters of the open-DFN models</caption><thead><tr><td>Model</td><td>Open fraction</td><td><i>p</i>P<sub>32</sub> borehole scale</td><td>Fracture size distribution</td></tr></thead><tbody><tr><td>Geo-DFN</td><td>0</td><td>4.6</td><td>UFM size model, <i>l</i><sub>c</sub> = 3.2 m</td></tr><tr><td><i>l</i><sub>c</sub>-open 25%</td><td>0.25</td><td>1.2</td><td>UFM size model, <i>l</i><sub>c</sub> = 17.6 m</td></tr><tr><td><i>l</i><sub>c</sub>-open 21%</td><td>0.21</td><td>1</td><td>UFM size model, <i>l</i><sub>c</sub> = 21.7 m</td></tr><tr><td>α-open</td><td>0.25</td><td>1.2</td><td>UFM size model, <i>l</i><sub>c</sub> = 3.2 m</td></tr><tr><td><i>pp</i>-open</td><td>0.21</td><td>1.02</td><td>Power-law, exponent −3.4</td></tr></tbody></table>

The open fraction and $P_{32}$, which are scale-dependent parameters, are calculated for a fracture size of 0.076 m equivalent to the borehole diameter. Both parameters are thus comparable to field measurements. For all the models, the smallest fracture for flow calculation is 2 m, and the system is a cube of 150 m side. The number of fractures is between 37,000 and 70,000 depending on the open-DFN model.

with a larger critical fracture size $l_c$ than the geo-DFN to account for the decrease of fracture density $P_{32}$. The first and second models differ by the proportion of open fraction. The second model—further referred to as α-open, is a geo-DFN for which a number of fractures have been removed to account for the open fraction independently of their size (Figure 8b). The last model, originally presented in Follin et al. [2014], relies on a different approach where the geo-DFN step was skipped and a open/hydro DFN was directly generated with only a single power-law fracture size distribution and no genetic base but only Poisson-point processes (no spatial correlations between fractures), with the orientation of Table 2 [Follin et al., 2014]. For simplicity, this last model is named *pp*-open (*pp* for poisson-point). The size distribution and the open fraction of the Geo-DFN, $l_c$-open, α-open and *pp*-open models are given in Figure 8 and the other model characteristics summarized in the legend of Table 3.

Although the 4 DFN models are consistent with the same data, they are very different in terms of connectivity, as highlighted by the analysis of the percolation parameter *p* as a function of system size *L* (Figure 9). This parameter measures the degree of connectivity of the networks and is expressed as [Charlaix et al., 1984, de Dreuzy et al., 2000]:

$$p(L) = \int_{l_{\min}}^{L} \frac{\pi^2}{8} l^3 n(l) dl \quad (10)$$

where *n*(*l*) is the density distribution of fracture size *l* per unit volume, and $l_{\min}$ is the smallest value of the distribution *n*(*l*). Bour and Davy [1998] demonstrate that the percolation parameter increases with scale and that the percolation threshold marks the limit scale above which statistically networks are connected by several fractures. The larger the percolation threshold, the larger the number of connected paths in a system. Figure 9 shows that the connection scale ranges from 10 to 50 m for the studied models,
<figure><img src="image_12.png" /><figcaption>**Figure 9.** Evolution of the percolation parameters as a function of the system size for the Geo-DFN, $l_c$-open, $\alpha$-open and $pp$-open models. The expression of $p$ is given by (10). The dashed rectangle indicates the percolation threshold for the prescribed orientations.</figcaption></figure>

and that $l_c$-open models are better connected than $\alpha$-open and $pp$-open model.

The Hydro-DFN is defined by assigning a transmissivity for each fracture. We define 4 different transmissivity models:

* all fractures have the same transmissivity ($T = 1$)—in this case, the resulting permeability is called the effective connectivity [de Dreuzy et al., 2001a];

* fracture transmissivity is proportional to the fracture size ($T = l$), which is the standard model studied by Frampton and Cvetkovic [2010] for similar geological formations;

* fracture transmissivity is a function of the normal stress applied to the fracture ($T = f(\sigma_n)$); it mimics the closure of fractures by the normal stress to the fracture plane $\sigma_n$; in practice, $\sigma_n$ is calculated from the projection of the remote stress tensor $\bar{\sigma}$ given in Glamheden et al. [2007]; the transmissivity-stress function is derived from Follin and Stigsson [2014] as:

$$f(\sigma_n) = \exp\left(-\frac{\sigma_n}{\sigma_c}\right)$$

with $\sigma_c = 4$ MPa;

* fracture transmissivity is assumed to be the product of the two previous relationships, i.e., fracture size and stress: $T = l f(\sigma_n)$.

This list does not cover the whole range of model possibilities but it gives good examples of the types of transmissivity relationships that could be investigated for a site modeling. Note that, in these relations, the transmissivities are relative with respect to a reference value that must be calibrated to the data. The following analysis is independent of this reference.

## 5.2. Comparison models/data, and sensitivity analysis

Flow simulations are performed with the computational suite “DFN.lab” [Le Goc et al., 2019], according to the setup described in Appendices A and B. Results are averaged from 10 realizations for $K_g$ (6) and the value is normalized by the permeability calculated first for a scale of 0.5 m: $K_g^*(L) = K_g(L)/K_g(L = 0.5)$.

Figure 10 shows the effects of the transmissivity model on the evolution of $K_g^*$. The first part of the curve is similar for all the transmissivity models with a $L^{-1}$ decrease. After the percolation threshold, $K_g^*$ stabilizes around a constant value for $T = 1$ (constant fracture transmissivity), and it grows for all other transmissivity models. The permeability increase after threshold is primarily induced by a size dependency of fracture transmissivity (see the models $T = l$ and $T = l \cdot f(\sigma)$); it is also induced, even though to a lesser degree, by a dependence of the transmissivity on stress ($T = f(\sigma)$) and thus on orientation. The combination of fracture size and stress/orientation dependency of fracture transmissivity ($T = l \cdot f(\sigma)$) has an effect on permeability scaling that is more than each of the property effects taken separately. This illustrates that the flow variability inducing the permeability increase is not a linear juxtaposition of effects.

The transition scale between both regimes has been estimated from the intersection of the end-member regimes at small and large scales: $K \propto L^{-1}$ at small scales, and $K \propto L^\alpha$ at large scales with $\alpha$ dependent on the transmissivity model. The transition is about independent of the transmissivity model around 4–4.5 m for the $l_c$-open model with an open fraction of 21% (Figure 10a) and around 10 m with 15% (Figure 10b). It is slightly larger than—but still
<figure><img src="image_13.png" /></figure>

**Figure 10.** Evolution of the geometric averaged permeability as a function of scale for the $l_c$-open model 21% (left) and 15% (right) with the various transmissivity models (see the list above). An approximate power law fit for the different regimes of scaling is provided for all the models. The fit exponent is an indication of the intensity of the permeability with scale for large scales. The permeability scaling measured at KFM08A (yellow dots) and the fit for large scaling (grey dashed lines) are also indicated for comparison.

<figure><img src="image_14.png" /></figure>

**Figure 11.** Evolution of the normalized geometric-averaged permeability for the different open-models and two transmissivity models: $T = 1$ (left) and $T = l \cdot f(\sigma)$ (right). The symbol colors and line types are similar to the one used in the Figure 9. The yellow dots on the right plot indicate the KFM08A scale analysis.

consistent with—the percolation scale for this model estimated from Figure 9.

These simulations are merely indications of the possible effects associated with the transmissivity models. A more complete study is needed to derive general laws. Note also that the normalization by the first point of the curve removes the variability from one model to another and from one realization to another. This point is discussed later.

How the permeability scaling reveals the open-model is presented in Figure 11 for two transmissivity models: $T = 1$ (left) and $T = l f(\sigma)$ (right). The

trends are similar to those presented in Figure 10, and the differences between the different open models mostly lie in the scale at which the transition between both two end-member regimes occurs. The increasing order of the different transition scales is in agreement with the evolution of the percolation parameter with the scale: $l_c$-open models first, then $pp$-open model and $\alpha$-open. Even the crossing of the curves of the last two models observed on $K_g^*$ is predicted by the percolation parameter. A comparison with the borehole data of KFM08A is made in Figure 11 (right). None of the presented open models
<figure><img src="image_15.png" /><figcaption>**Figure 12.** Evolution of the standard deviation for 10 realizations of the lognormal distribution of permeability as a function of scale. Left: comparison of the different open-DFN models, all with the transmissivity model $T = l$; right: comparison of the different transmissivity models, all with the $l_c$-open DFN model.</figcaption></figure>

correctly reproduce the increase of permeability with scale after the transition scale for the two transmissivity models presented, $T = 1$ and $T = lf(\sigma)$. The $\alpha$-open model is too weakly connected at large scale irrespective of the fracture transmissivity model, while *pp*-open model is slightly too connected at small scales. The $l_c$-open model with an open fraction of 15% and fracture transmissivity that scales with the fracture size is a good candidate to reproduce the permeability scaling at small and large scales (Figure 10b).

In Figure 12, we present an attempt to estimate the model variability (third validation indicator in Section 4.4). We calculate the lognormal standard deviation of the geometric average permeability $\sigma_{\ln(K_g)}$ for 10 realizations for models, and in-between boreholes for data. $\sigma_{\ln(K_g)}$ is dimensionless and does not depend on the unit of the fracture transmissivity, which makes it a useful indicator for comparing models and data. $\sigma_{\ln(K_g)}$ has been calculated for 10 realizations for each combination between an open model and a transmissivity model. Figure 12 shows the evolution of $\sigma_{\ln(K_g)}$ with scale for different transmissivity models and $l_c$-open 15% (left), and for different open models and $T = lf(\sigma)$ (right). Not surprisingly, the transmissivity model with the largest variability is $T = lf(\sigma)$, i.e., when the transmissivity depends both on stress and fracture size. This is the only model that shows an increase in variability with scale above 10–20 m. This trend is also observed in the data but to a lesser degree. For all other

transmissivity models, the variability is about constant with scale and it is not very different from one model to another except for the model with a constant transmissivity which has a lower variability (Figure 12a). This result must be confirmed by a larger number of simulations.

The impact of the open model on $\sigma_{\ln(K_g)}$ is presented in Figure 12b. The $\alpha$-open model has the largest variability, probably because it is largely controlled by a few big fractures. The *pp*-open model is the only one where the variability does not increase with scale for this transmissivity model.

Note that the comparison between model and data is not fair because the boreholes mix domains that can be quite heterogeneous, in particular, in terms of open fraction, while models are built with a constant open fraction. The $\sigma_{\ln(K_g)}$ values are therefore likely to be overestimated for the data compared to the models and a correction may be necessary to make a true comparison.

## 6. Discussion

The four indicators defined in the Section 4.4 turn out to be useful to validate or calibrate models. All the models presented are consistent with the data, but no all of them can reproduce the data indicators. We discuss here the ability of the indicators to constrain both the structure that carries the flow and the transmissivity model that is applied to that structure.
All the models—again, calibrated on the data from
Forsmark—reproduce the V-shape scaling curve of
the permeability geometric average observed in data.

* The first indicator is the scale which marks the transition between the $L^{-1}$ scaling at small scale and the transmissivity scaling at large scale. This scale is closely related to a connectivity threshold of the open-DFN structure. Indeed, for open-DFN, whose distribution of the largest fractures follows a power law, the connectivity increases with scale and the percolation threshold is reached at a given scale which can be predicted by the evolution of the percolation parameter with the scale. For all the tested models, the transition/percolation scale seems to be independent of the transmissivity model.

* The second indicator is the increase of the geometrically-averaged permeability with scale. It can be modeled by a power-law scaling whose exponent gives the increase in intensity. This indicator reveals the transmissivity model and a large increase is likely associated with a dependency of fracture transmissivity with fracture size. A high variability in fracture transmissivity, even without this transmissivity scaling effect, can also create an increase in permeability but to a lesser extent.

* The third indicator is the variability of geometrically-averaged permeability from one well to another (data) or from one realization to another (models). It has been quantified as the standard deviation of the logarithm of permeability, which is a dimensionless parameter that relies on the likely lognormal permeability distribution. This indicator is much more sensitive to the structure model (open-DFN) than the transmissivity model with the exception of the extreme transmissivity model $T = l \cdot f(\sigma)$, where the variability is high and scale dependent above the transition scale (i.e., 2nd indicator) as it is the case for the data. We do not want to draw definitive conclusions on this indicator because it would require many more models and realizations, and a

comparison taken on more homogeneous
data.

Although the analysis is far from exploring all the
structure and transmissivity models of DFN, we can
draw few preliminary conclusions on the capacity of
the chosen models to reproduce the indicators.

* In the DFN modeling workflow, an under-constrained critical point is the way in which the open fracture fraction is distributed in the network. Only the total open fraction cannot be retrieved from flow logs. If the open fraction is uniformly distributed over all fractures regardless of their size or orientation, as it is for the α-open model, the network is less connected than other models, and probably the natural networks, whatever the assigned transmissivity model may be (see Figure 8a). On the other hand, if the open fraction is large at small scale, as it is for the pp-open model, the network is more connected than other models and data at these scales. This highlights the primary importance of properly defining what fractures are sealed and the proportion of open fractures per fracture size.

* The transmissivity model $T = 1$ (i.e., constant transmissivity for each fracture) cannot not reproduce the V-shape of $K_g$ vs scale. This model is obviously unrealistic, however we suspect that any transmissivity model, whose variability from one fracture to another is too low could not reproduce the permeability increase with scale as observed in data above ~10 m. On the other hand, the transmissivity model combining fracture size and normal stress $T = l f(\sigma)$ predicts permeability increase with scale too sharp compared to data. It is not certain to find an unequivocal solution between the transmissivity model and the observed permeability scaling by playing on the fracture size and orientation/stress dependence, but this indicator still reduces the field of possibilities.

* It is difficult to know at present whether a
model should be rejected altogether (e.g.,
the α-open or pp-open models) or whether
it is still possible to modify some of their
parameters to fit the data. This exercise is
however necessary to limit the field of possible models.

* The best-fit model among those tested is $l_c$-open with a transmissivity model $T = l$ but it is not totally satisfactory, particularly in terms of the variability it induces (3rd indicator), which is too low and does not depend on the scale.

Note that the fourth indicator, which is the existence of a representative elementary volume, is not observed in the data where the permeability geometric average increases continuously with the scale.

## 7. Conclusion

The objective of the paper is to better understand and quantify the flow structure in fractured rocks from PFL flow logs, and also to propose relevant indicators for validating, calibrating or even rejecting hydrological models. The data consists of a series of inflows along each borehole with an equivalent transmissivity measured from pumping tests. We first studied what the inflow distribution tells us about the permeability structure from a series of analysis: distribution of transmissivities as a function of depth, proportion of flowing sections as a function of section scale, scaling of the arithmetically-averaged and geometrically-averaged permeability. We argue that the permeability scaling evolution provides more information on the flow structures than macroscopic values (i.e., large scale permeability), the latter being too integrative.

We thus define three indicators of the flow structure from the scale evolution of the permeability geometric average that highlight intrinsic properties of the flow structure:

* a percolation scale $l_s$,

* the way permeability increases with scale above $l_s$,

* the permeability variability, measured as the standard deviation of the permeability logarithm for different boreholes or model realizations, as a function of scale.

A 4th indicator on the representative elemental volume could in principle be defined but the data show that this volume/scale is beyond the 300 m investigated.

We tested a series of numerical models by running flow simulations on generated DFNs to compute the same indicators on synthetic data. Most of the models are built from the DFN methodology with three steps: the geo-DFN based on the observed fracture network, the open-DFN which is the part of the geo-DFN where fractures are open, and a transmissivity model applying on each fracture of the open-DFN. The analysis of the models showed that the percolation scale is controlled by the open-DFN structure and that the percolation scale can be predicted from a scale analysis of the percolation parameter (basically, the third moment of the fracture size distribution that provides a measure of the network connectivity). For the same geo-DFN, the distribution of the open fraction with fracture size is the critical parameter that controls the percolation parameter and the percolation threshold. The way permeability increases with scale above the percolation threshold is controlled by the transmissivity model. We have tested several of them with a dependence of the fracture transmissivity on either the orientation of the fractures via a stress-controlled transmissivity or their size or both. Although size dependence induces the greatest effect on permeability increase, the variability due to orientation (i.e., stress) dependence also has a non-negligible effect that is enhanced when the two dependencies are combined.

The comparison with data on the first two indicators shows that a model that matches the characteristics of the geo-DFN with the measured open fraction of 15% adequately fits the data, provided that the large fractures remain open and that the fracture transmissivity model is well selected. Most of the other models shows unacceptable differences with data, however other models or model combinations still have to be explored before rejecting them.

The third indicator on model variability is still problematic since the natural data show a higher variability than the models but the open fraction is also much more variable in the data than in the models.

By characterizing the flow structure in terms of permeability scaling and variability, and by defining three indicators that describe some fundamental characteristics of the flow/permeability whatever the scale may be, we bring a new way to calibrate and/or validate models or at least to reduce the field of model possibilities.
**Conflicts of interest**

Authors have no conflict of interest to declare.

**Acknowledgements**

The authors acknowledge the two anonymous ref-
erees who helped improve the manuscript. Philippe
Davy is particularly grateful to Professor Ghislain de
Marsily for the help and encouragement to develop
the fundamental research on the flow properties of
fractured media.

**Appendix A. DFN geometrical model definition**

We consider the fractured rock mass as a set of frac-
tures represented by discrete elements embedded in
a linear elastic media, as commonly described by
the Discrete Fracture Network (DFN) methodology.
In 3D, fractures are represented by surfaces (usually
disks), embedded in a volume. We hereafter refer to
two types of DFN models, the *Poisson Point process-based* and the *genetic (Universal Fracture Model) UFM* ones.

The Poisson Point process-based DFN model—
named Poisson DFN model—relies only on statistics-
based fracture generation. The DFN properties are
predefined by statistical distributions (density distri-
butions of fracture sizes and orientations) and frac-
tures are generated from a random sampling of ge-
ometrical properties (Poisson Point process-based),
where each fracture position, size, and orientation
are independent of the other fractures, but the aver-
age density is predetermined.

The UFM DFN model introduced by Davy et al.
[2013] belongs to the range of genetic DFN models,
introducing a fracture hierarchy between fractures by
mimicking the fracturing process. The UFM model
produces fracture size distributions and fracture in-
tersections that are consistent with observations, us-
ing simple rules (N–G–A) for:

* Nucleation: fractures are introduced in the system at a constant nucleation rate *n*̇*<sub>N</sub>* with a nuclei size of *l*<sub>n</sub>.

* Growth: once created, fractures grow follow-
ing a power-law relationship that describes
the crack tip velocity v(l) = Cl^{a_G}, with l the

fracture size, *C* the growth rate, and *a*<sub>G</sub> the
growth exponent.

* Arrest: while growing, fractures stop if they encounter a larger one, which is a first-order approximation of mechanical interactions.

While the DFN fracture size distribution is imposed
from the beginning in a Poisson DFN model, it nat-
urally arises from the growth and arrest processes in
the UFM, as part of a spatially organized structure.
For constant nucleation rate and growing fractures,
there is a stationary solution for the fracture size dis-
tribution (referred to as the dilute regime) described
by $n_G(l) = n_N/C \cdot l^{-a_G}$, while the arrest rule is re-
sponsible for the emergence of a quasi-universal self-
similar fracture size distribution (referred to as the
dense regime) described by $n_A(l) = D\gamma^D l^{-D+1}$, with
*D* the topological dimension associated to fracture
centers, and *γ*, a geometrical parameter dependent
on fracture orientations. The model thus results in a
two power-law size distribution, where the transition
size $l_c$ between $n_G$ and $n_A$ is both the scale at which
the network is connected and the average size of frac-
ture blocks. Moreover, the arrest rule is responsible
for the emergence of T intersections, where a fracture
abuts on another one. The fracture spatial organiza-
tion introduced by the UFM framework has an im-
pact on both the network topology [Hope et al., 2015]
and its hydrogeological behavior [Maillot et al., 2016].

**Appendix B. Numerical set-up**

The DFN and flow simulation volume is a generic
3D cubic domain. Boundary conditions are assigned
all over the domain sides, either fixed heads or no
flow and possibly along a generic borehole (like a
1D line) added in the cube. Steady-state flow tests
are defined with PFL-like conditions to mimic the in-
situ Posiva Flow Log (PFL) tests done at the Forsmark
site and described in Section 2. A head drawdown of
-10 is imposed at a vertical borehole at the center of
the domain to simulate the constant head pumping,
while zero heads are imposed on the vertical sides
and the top and no-flow on the bottom. The steady-
state flow-field was calculated using DFN.lab, and
the inflows to the borehole at each connected open
fracture intersected by the borehole were retrieved.
The borehole itself was represented by a vertical line.

The first step of the flow simulation is to iden-
tify the set of fractures of the DFN that constitute a
<figure><img src="image_16.png" /><figcaption>**Figure 13.** Clipped view along a vertical central plane going through the vertical borehole for PFL-like boundary conditions (a) flow rate (in m³/s) and (b) head (in m) fields.</figcaption></figure>

connected path between the defined active hydraulic boundaries. If there is no connected path, there is simply no flow simulation. Once the set of connected fractures is identified, it is cleaned up from its dead-ends (we assume that a fracture with only one intersection to the connected cluster or hydraulic boundary does not contribute to flow). This cleaned fracture cluster can be called a flow backbone. Finally, fractures of the flow backbone are meshed, and flow is solved using a mixed-hybrid finite-element scheme in steady-state conditions. (Figure 13).

For each intersection of the borehole with a flowing fracture f, one defines the observed fracture transmissivity $T_{PFL,f}$ as:

$$T_{PFL,f} = \frac{Q_{PFL,f}}{2\pi\Delta h} \cdot \ln\left(c \cdot \frac{R}{r_m}\right) \quad (11)$$

with $Q_{PFL,f}$ the flow rate at the borehole intersection, $r_m/c$ (with c equal to 10) the equivalent borehole radius and R the radial distance between the borehole and the boundary conditions. c is a correction factor that accounts for numerical artifacts relative to the borehole type condition and meshes element shape (triangle). With this definition, the transmissivity of an isolated fracture connected to the borehole and

the system limits can be inferred from the simulated values of $Q_{PFL,f}$ and $\Delta h$. If the intercepted fracture is not isolated, measured equivalent transmissivity rather reflects the flow capacity of the closely connected cluster of fractures around the intercepted one rather than the transmissivity of the intercepted fracture itself.

The transmissivity $T_{PFL}(s)$, of a borehole section of size s, including potentially several flowing fractures. It is defined from (11) with $Q_{PFL}(s)$ the sum of the flow rates over the individual fractures of the section.

The equivalent hydraulic conductivity of the fracture network (in m/s) estimated from a borehole section is equal to the sum of inflow transmissivities divided by the section length.

## References

* Berkowitz, B., Bour, O., Davy, P., and Odling, N. (2000). Scaling of fracture connectivity in geological formations. *Geophys. Res. Lett.*, 27(14), 2061–2064.

* Bonnet, E., Bour, O., Odling, N. E., Davy, P., Main, I., Cowie, P., and Berkowitz, B. (2001). Scaling of fracture systems in geological media. *Rev. Geophys.*, 39(3), 347–383.
Bour, O. and Davy, P. (1997). Connectivity of random fault networks following a power law fault length distribution. *Water Resour. Res.*, 33(7), 1567–1583.

Bour, O. and Davy, P. (1998). On the connectivity of three-dimensional fault networks. *Water Resour. Res.*, 34(10), 2611–2622.

Cacas, M. C., Ledoux, E., de Marsily, G., Barbreau, A., Calmels, P., Gaillard, B., and Margritta, R. (1990a). Modeling fracture flow with a stochastic discrete fracture network: Calibration and validation: 2. The transport model. *Water Resour. Res.*, 26(3), 491–500.

Cacas, M. C., Ledoux, E., de Marsily, G., Tillie, B., Barbreau, A., Durand, E., Feuga, B., and Peaudecerf, P. (1990b). Modeling fracture flow with a stochastic discrete fracture network: calibration and validation: 1. The flow model. *Water Resour. Res.*, 26(3), 479–489.

Charlaix, E., Guyon, E., and Rivier, N. (1984). A criterion for percolation threshold in a random array of plates. *Solid State Commun.*, 50(11), 999–1002.

Charlaix, E., Guyon, E., and Roux, S. (1987). Permeability of a random array of fractures of widely varying apertures. *Transp. Porous Media*, 2, 31–43.

Chen, T., Clauser, C., Marquart, G., Willbrand, K., and Hiller, T. (2018). Upscaling permeability for three-dimensional fractured porous rocks with the multiple boundary method. *Hydrogeol. J.*, 26(6), 1903–1916.

Chen, T., Clauser, C., Marquart, G., Willbrand, K., and Mottaghy, D. (2015). A new upscaling method for fractured porous media. *Adv. Water Resour.*, 80, 60–68.

Clauser, C. (1992). Permeability of crystalline rocks. *Eos, Trans. Am. Geophys. Union*, 73(21), 233–238.

Darcel, C., Bour, O., and Davy, P. (2003a). Stereological analysis of fractal fracture networks. *J. Geophys. Res.*, 108(B9), 13–1–13–14.

Darcel, C., Bour, O., Davy, P., and de Dreuzy, J. R. (2003b). Connectivity properties of two-dimensional fracture networks with stochastic fractal correlation. *Water Resour. Res.*, 39(10), article no. 1272.

Darcel, C., Davy, P., Le Goc, R., Bour, O., and de Dreuzy, J. R. (2009). Statistical methodology for discrete fracture models – including fracture size, orientation uncertainty together with intensiycertainty and variability. Technical Report R-09-38. Retrieved from Stockholm: http://www.skb.com/

publication/1983611/R-09-38.pdf.

Davy, P., Bour, O., De Dreuzy, J.-R., and Darcel, C. (2006). Flow in multiscale fractal fracture networks. In *Fractal Analysis for Natural Hazards*, volume 261 of *Geol. Soc. London, Spec. Publ.*, pages 31–45. Geological Society of London.

Davy, P., Darcel, C., Le Goc, R., Munier, R., Selroos, J.-O., and Mas Ivars, D. (2018). DFN, why, how and what for, concepts, theories and issues. In *Paper presented at the 2nd International Discrete Fracture Network Engineering Conference*. OnePetro, Richardson, TX.

Davy, P., Le Goc, R., and Darcel, C. (2013). A model of fracture nucleation, growth and arrest, and consequences for fracture density and scaling. *J. Geophys. Res. Solid Earth*, 118(4), 1393–1407.

Davy, P., Le Goc, R., Darcel, C., Bour, O., de Dreuzy, J.-R., and Munier, R. (2010). A likely universal model of fracture scaling and its consequence for crustal hydromechanics. *J. Geophys. Res.*, 115(B10), 1–13.

de Dreuzy, J. R., Davy, P., and Bour, O. (2000). Percolation parameter and percolation-threshold estimates for three-dimensional random ellipses with widely scattered distributions of eccentricity and size. *Phys. Rev. E*, 62(5), 5948–5952.

de Dreuzy, J. R., Davy, P., and Bour, O. (2001a). Hydraulic properties of two-dimensional random fracture networks following a power law length distribution 1. Effective connectivity. *Water Resour. Res.*, 37(8), 2065–2078.

de Dreuzy, J. R., Davy, P., and Bour, O. (2001b). Hydraulic properties of two-dimensional random fracture networks following a power law length distribution 2. Permeability of networks based on log-normal distribution of apertures. *Water Resour. Res.*, 37(8), 2079–2096.

de Dreuzy, J. R., Davy, P., and Bour, O. (2002). Hydraulic properties of two-dimensional random fracture networks following power law distributions of length and aperture. *Water Resour. Res.*, 38(12), 12–1–12–9.

de Dreuzy, J. R., de Boiry, P., Pichot, G., and Davy, P. (2010). Use of power averaging for quantifying the influence of structure organization on permeability upscaling in on-lattice networks under mean parallel flow. *Water Resour. Res.*, 46(8), 1–11.

De Dreuzy, J. R., Pichot, G., Poirriez, B., and Erhel, J. (2013). Synthetic benchmark for modeling flow in 3D fractured media. *Comput. Geosci.*, 50,
59–71.

Dershowitz, W. S. and Herda, H. H. (1992). Interpretation of fracture spacing and intensity. In *Paper presented at the 33th US Symposium on Rock Mechanics (USRMS)*. OnePetro, Richardson, TX.

Dewandel, B., Lachassagne, P., Wyns, R., Marechal, J. C., and Krishnamurthy, N. S. (2006). A generalized 3-D geological and hydrogeological conceptual model of granite aquifers controlled by single or multiphase weathering. *J. Hydrol.*, 330(1–2), 260–284.

Doolaeghe Wehowsky, D. (2021). *Colmatage des réseaux de fractures, modèles et conséquences hydrologiques*. PhD thesis, University of Rennes 1, Rennes.

Follin, S. (2008). Bedrock hydrogeology Forsmark, Site descriptive modelling, SDM-Site Forsmark. Technical Report R-08-95, Svensk Kärnbränslehantering AB (SKB), Stockholm, Sweden.

Follin, S., Hartley, L., Rhén, I., Jackson, P., Joyce, S., Roberts, D., and Swift, B. (2014). A methodology to constrain the parameters of a hydrogeological discrete fracture network model for sparsely fractured crystalline rock, exemplified by data from the proposed high-level nuclear waste repository site at Forsmark, Sweden. *Hydrogeol. J.*, 22(2), 313–331.

Follin, S., Leven, J., Hartley, L., Jackson, P., Joyce, S., Roberts, D., and Swift, B. (2007). Hydrogeological characterisation and modelling of deformation zones and fracture domains, Forsmark modelling stage 2.2. Technical Report R-07-48, Svensk Kärnbränslehantering AB (SKB), Stockholm, Sweden.

Follin, S. and Stigsson, M. (2014). A transmissivity model for deformation zones in fractured crystalline rock and its possible correlation to in situ stress at the proposed high-level nuclear waste repository site at Forsmark, Sweden. *Hydrogeol. J.*, 22(2), 299–311.

Fox, A., La Pointe, P., Hermanson, J., and Öhman, J. (2007). Statistical geological discrete fracture network model. Forsmark modelling stage 2.2. Technical Report R-07-46, Svensk Kärnbränslehantering AB (SKB), Stockholm, Sweden.

Frampton, A. and Cvetkovic, V. (2010). Inference of field-scale fracture transmissivities in crystalline rock using flow log measurements. *Water Resour. Res.*, 46(11), article no. W11502.

Glamheden, R., Fälth, B., Jacobsson, L., Harrström, J., Berglund, J., and Bergkvist, I. (2010). Counter-

force applied to prevent spalling. Technical Report TR-10-37, Svensk Kärnbränslehantering AB (SKB), Stockholm, Sweden.

Glamheden, R., Fredriksson, A., Roeshoff, K., Karlsson, J., Hakami, H., and Christiansson, R. (2007). Rock mechanics Forsmark. Site descriptive modelling Forsmark stage 2.2. Technical report, Swedish Nuclear Fuel and Waste Management, Sweden.

Hope, S. M., Davy, P., Maillot, J., Le Goc, R., and Hansen, A. (2015). Topological impact of constrained fracture growth. *Front. Phys.*, 3, article no. 75.

Hunt, A. G. (2003a). Reply to comment by S. P. Neuman on "Some comments on the scale dependence of the hydraulic conductivity in the presence of nested heterogeneity". *Adv. Water Resour.*, 26(11), 1215.

Hunt, A. G. (2003b). Some comments on the scale dependence of the hydraulic conductivity in the presence of nested heterogeneity. *Adv. Water Resour.*, 26(1), 71–77.

Illman, W. A. (2006). Strong field evidence of directional permeability scale effect in fractured rock. *J. Hydrol.*, 319(1–4), 227–236.

Le Goc, R., Pinier, B., Darcel, C., Lavoine, E., Doolaeghe, D., De Simone, S., de Dreuzy, J.-R., and Davy, P. (2019). DFN.lab: software platform for Discrete Fracture Network models. In *Paper presented at the AGU Fall Meeting 2019*. American Geophysical Union, Washington, DC.

Lei, Q., Latham, J.-P., and Tsang, C.-F. (2017). The use of discrete fracture networks for modelling coupled geomechanical and hydrological behaviour of fractured rocks. *Comput. Geotech.*, 85, 151–176.

Long, J. and Billaux, D. M. (1987). From field data to fracture network modeling: an example incorporating spatial structure. *Water Resour. Res.*, 23(7), 1201–1216.

Long, J. C. S. and Witherspoon, P. A. (1985). The relationship of the degree of interconnection to permeability in fracture networks. *J. Geophys. Res.*, 90(B4), 3087–3098.

Maillot, J., Davy, P., Le Goc, R., Darcel, C., and de Dreuzy, J. R. (2016). Connectivity, permeability, and channeling in randomly distributed and kinematically defined discrete fracture network models. *Water Resour. Res.*, 52(11), 8526–8545.

Mandelbrot, B. B. (1982). *The Fractal Geometry of*
*Nature*. W.H. Freeman, New-York.

Martin, D. and Follin, S. (2011). Review of possible correlations between in situ stress and PFL fracture transmissivity data at Forsmark. Technical Report R-08-69, Svensk Kärnbränslehantering AB (SKB), Stockholm, Sweden.

Martinez-Landa, L. and Carrera, J. (2005). An analysis of hydraulic conductivity scale effects in granite (Full-scale Engineered Barrier Experiment (FEBEX), Grimsel, Switzerland). *Water Resour. Res.*, 41(3), article no. W03006.

Meier, P. M., Carrera, J., and Sanchez-Vila, X. (1999). A numerical study on the relationship between transmissivity and specific capacity in heterogeneous aquifers. *Groundwater*, 37(4), 611–617.

Neuman, S. (2005). Trends, prospects and challenges in quantifying flow and transport through fractured rocks. *Hydrogeol. J.*, 13(1), 124–147.

Neuman, S. P. (1994). Generalized scaling of permeabilities. *Geophys. Res. Lett.*, 21(5), 349–352.

Neuman, S. P. (2003). Comment on “Some comments on the scale dependence of the hydraulic conductivity in the presence of nested heterogeneity” by A. G. Hunt. *Adv. Water Resour.*, 26(11), 1213.

Neuman, S. P. and Di Federico, V. (2003). Multifaceted nature of hydrogeologic scaling and its interpretation. *Rev. Geophys.*, 41(3), 4–1–4–31.

Oda, M. (1985). Permeability tensor for discontinuous rock masses. *Géotechnique*, 35(4), 483–495.

Öhberg, A. and Rouhiainen, P. (2000). Posiva ground-water flow measuring techniques. Technical Report 2000-12, Posiva Oy, Helsinki, Finland. Retrieved from Helsinki.

Olofsson, I., Simeonov, A., Stephens, M., Follin, S., Nilsson, A., Röshoff, K., Lindberg, U., Lanaro, F., Fredriksson, A., and Persson, L. (2007). Site descriptive modelling Forsmark, stage 2.2. Technical Report R-07-15, Svensk Kärnbränslehantering AB (SKB), Stockholm, Sweden.

Park, B. Y., Kim, K. S., Kwon, S., Kim, C., Bae, D. S., Hartley, L. J., and Lee, H. K. (2002). Determination of the hydraulic conductivity components using a three-dimensional fracture network model in volcanic rock. *Eng. Geol.*, 66(1), 127–141.

Pavón-Domínguez, P. and Moreno-Pulido, S. (2022). Sandbox fixed-mass algorithm for multifractal unweighted complex networks. *Chaos Solitons Fractals*, 156, article no. 111836.

Perrochet, P. (2005). A simple solution to tunnel or well discharge under constant drawdown. *Hydrogeol. J.*, 13(5), 886–888.

Piggott, A. (1997). Fractal relations for the diameter and trace length of disc-shaped fractures. *J. Geophys. Res.*, 102(B8), 18121–18125.

Ren, S., Zhang, Y., Jim Yeh, T. C., Wang, Y., and Carr, B. J. (2021). Multiscale hydraulic conductivity characterization in a fractured granitic aquifer: the evaluation of scale effect. *Water Resour. Res.*, 57(9), article no. e2020WR028482.

Renard, P. and Marsily, G. d. (1997). Calculating equivalent permeability: a review. *Adv. Water Resour.*, 20(5–6), 253–278.

Rouhiainen, P. and Pöllänen, J. (2003). Forsmark site investigation Difference flow logging of borehole KFM01A. Technical Report P-03-28, Posiva Oy, Helsinki, Finland.

Rouhiainen, P., Pöllänen, J., and Tec-Oy, P. (2004). Forsmark site investigation: Difference flow logging in borehole KFM06A. Technical Report P-05-15, Posiva Oy, Helsinki, Finland.

Selroos, J.-O., Ivars, D. M., Munier, R., Hartley, L., Libby, S., Davy, P., Darcel, C., and Trinchero, P. (2022). Methodology for discrete fracture network modelling of the Forsmark site. Part 1—concepts, data and interpretation methods. Technical Report R-20-11, Svensk Kärnbränslehantering AB (SKB), Stockholm, Sweden.

SKB (2011). Long-term safety for the final repository for spent nuclear fuel at Forsmark. Main Report of the SR-Site Project. Technical Report TR-11-01, Svensk Kärnbränslehantering AB (SKB), Stockholm, Sweden.

Tél, T., Fülöp, A., and Vicsek, T. (1989). Determination of fractal dimensions for geometrical multifractals. *Phys. A: Stat. Mech. Appl.*, 159(2), 155–166.

Wen, X.-H. and Gómez-Hernández, J. J. (1996). Upscaling hydraulic conductivities in heterogeneous media: An overview. *J. Hydrol.*, 183(1), ix–xxxii.

Zou, L. and Cvetkovic, V. (2020). Inference of transmissivity in crystalline rock using flow logs under steady-state pumping: impact of multiscale heterogeneity. *Water Resour. Res.*, 56(8), article no. e2020WR027254.

Zou, L. and Cvetkovic, V. (2021). Evaluation of flowlog data from crystalline rocks with steady-state pumping and ambient flow. *Geophys. Res. Lett.*, 48(9), article no. e2021GL092741.