Quantitative assessment of the flow properties and mechanical stability of naturally fractured rock is frequently practised across the mining, petroleum, geothermal, geological disposal, construction and environmental remediation industries. These fluid and mechanical behaviours are strongly influenced by the connectivity of the fracture system and the size of the intact rock blocks. However, these are amongst the more difficult fracture system properties to characterize and honour in numerical simulations. Nonetheless, they are still the product of interactions between fractures that can be conceptualized as a series of deformation events following geomechanical principles. Generating numerical models of fracture networks by simulating this deformation with a coupled and evolving rock mass and stress field is a significant undertaking. Instead, large-scale fracture network models can be ‘grown’ dynamically according to rules that mimic the underlying mechanical processes and deformation history. This paper explores a computationally efficient rules-based method to generate fracture networks, demonstrates how different types of fracture patterns can be simulated, and illustrates how inclusion of fracture interactions can affect flow and mechanical properties. Relative to methods without fracture interaction and in contrast to some other rules-based approaches, the method described here regularizes and increases fracture connectivity and decreases flow channelling.

Naturally occurring rock fractures can form fluid conduits and affect the mechanical properties of the subsurface. Characterizing their properties is therefore an important task in many engineering applications such as rock excavation for mines and tunnels, hydrocarbon production, carbon storage, mineral extraction, disposal of spent nuclear fuel, and deep geothermal reservoirs. Not only are the properties of individual fractures of interest but also how they interconnect, as this is a primary control on bulk flows, fluid–rock interactions and the strength of the rock over a range of scales relevant to engineering applications (NERC 1996; Tsang and Neretnieks 1998; Zimmerman and Main 2004; Neuman 2005; Hoek 2023). The discrete fracture network (DFN) approach provides a commonly used conceptual framework for the characterization of fractured rock (Long et al. 1985; Dershowitz and Einstein 1988; Cacas et al. 1990; Berkowitz 2002; Dershowitz et al. 2004; Lei et al. 2017; Selroos et al. 2022) whereby the properties of both individual fractures and the networks they form are described quantitatively.

In a simple form, the geometrical properties (location, orientation, size and shape) of individual fractures in a DFN can be described by probability distribution functions, with fracture centres typically located by a Poisson point process according to a characterized intensity (Dershowitz and Einstein 1988). These properties can be readily calibrated by combining data from borehole logs and exposed rock faces. Generating such a DFN can then be accomplished through a simple statistics-based algorithm whereby each fracture is generated independently at very low computational cost (subsequently referred to as a static model or static DFN). However, natural fracture patterns demonstrate fracture interaction such as geometrical termination and spatial variations in fracture intensity associated with other fractures. Such interactions and variability have significant implications for network connectivity (Englman et al. 1983; Pollard and Aydin 1988; Olson and Pollard 1989; Renshaw 1993; Du and Aydin 1995; de Dreuzy et al. 2000; Sanderson and Nixon 2015; Sævik and Nixon 2017; Laubach et al. 2019; Sherman et al. 2020). Honouring these characteristics in DFN models requires some form of dynamic simulation of the rock's deformation history. This can be either be an explicit simulation of the evolving rock-mass stresses and failure (Olson 2007; Paluszny and Matthäi 2009; Paluszny and Zimmerman 2013; Li and Zhang 2021), a genetic approach whereby the fracture system is grown according to geomechanics-based rules for fracture nucleation, propagation and interactions through termination or coalescence (Bonneau et al. 2013; Davy et al. 2013; Maillot et al. 2014; Libby et al. 2019), or a combination of the two (Lavoine et al. 2020). By doing so, the model can more reliably describe several important characteristics of observed fracture networks (Lei et al. 2017):

  • the connectivity of the system, affecting hydraulic communication;

  • the anisotropy and clustering of the system, affecting channelling of fluid migration and the mechanical anisotropy of the rock mass; and

  • the size and shape distributions of rock blocks between discontinuities, affecting the stability and deformability of the rock.

Simulating the growth of a fracture network and its effects on rock-mass properties coupled to an evolving stress field is computationally demanding, which has limited the approach to 2D trace maps (Renshaw 1996; Olson 2007; Paluszny and Matthäi 2009), 3D DFNs with fracture numbers in the hundreds (Paluszny and Zimmerman 2013; Welch et al. 2019; Paluszny et al. 2020; Thomas et al. 2020) or hybrid rules-based approaches with components of stress modelling (Lavoine et al. 2020). It is also a forward-modelling approach in which it is difficult to determine the model parameters required to reproduce observed fractures and their statistics. Such models are therefore best used as an aid to understand the basic processes and results of fracture–fracture and fracture–rock interactions.

Generating DFNs using rules for the nucleation, propagation and termination or coalescence has therefore been considered a suitable compromise to generate large DFNs displaying more geologically realistic fracture patterns (Gringarten 1998; Bourne et al. 2001; Hoffmann et al. 2004; Srivastava et al. 2005; Davy et al. 2010, 2013; Bonneau et al. 2013; Maillot et al. 2016; Lavoine et al. 2020). By doing this without explicitly modelling the full physics of the stress evolution can achieve a significant reduction in computational cost. This is important because whilst DFNs with less than 100 000 fractures are useful for improving our understanding of rules-based behaviours, many practical applications of subsurface flow models require simulation of rock volumes measuring hundreds of metres in each dimension whilst including fractures on the metre scale. Even with fracture intensities at the percolation threshold, such applications can demand more than a million discrete fractures. The associated computational cost is exacerbated by the stochasticity of the DFNs requiring iterative Monte Carlo analyses.

The algorithm first described in Libby et al. (2019) uses mechanistic rules to facilitate better replication of natural fracture patterns, which can produce more geologically realistic DFNs. This algorithm grows fractures during a series of discrete time steps and the resulting DFNs are termed grown DFNs (GDFNs). This is implemented in the FracMan software (WSP UK Ltd 2023). A key feature of the GDFN algorithm is its computational efficiency, allowing it to readily produce DFNs with fracture counts of more than 10 million on desktop hardware.

Here, we describe some of the mechanistic rules that comprise the algorithm with comparison to other rules-based fracture generation models. Example outputs are then presented that illustrate how the GDFN can be configured to match a range of connectivity characteristics, before some of the hydrological and mechanical properties of a GDFN model are compared to an equivalent static DFN. Results that contrast with another published model are explored and a potential explanation for the discrepancy offered. Some representative fracture generation timings are then presented, recognizing that such timings are highly dependent on the host hardware and software configuration.

Nucleation and time steps

GDFN fractures are nucleated and grown in a series of time steps to allow interaction during their formation, as per other DFN generation methods governed by mechanistic rules (see the discussion in Lei et al. 2017). Each fracture set in a GDFN model typically corresponds to a deformation phase, and each subsequent set can interact with those fractures that precede it, mimicking natural deformation processes. GDFN fractures grow from a nucleation point stochastically located in space (similarly to Davy et al. 2010, 2013). ‘Rays’ grow away from this point on a plane orientated as per the assigned orientation of the fracture (Fig. 1), similar to Paluszny and Zimmerman (2012). The ends of these rays define the vertices of the polygon that comprise the fracture geometry.

Fractures are assigned a target size when they nucleate, as per a user-input size distribution. Both the total number of computed time steps required to reach a target intensity (termed total time steps) and the number of time steps each fracture takes to reach its target size (termed growth time steps) may be adjusted according to the desired mode of fracture growth. Increasing the number of total time steps reduces how many fractures are growing simultaneously, decreasing the probability of interaction between growing fractures without changing the probability of interaction with fractures that have finished growing (and vice versa). Increasing the number of growth time steps increases the probability of interaction with other fractures that are still growing (and vice versa).

In the case where fracture growth is not affected by other factors (discussed subsequently), each fracture will grow a proportion of its target size such that each ray on the fracture grows an equal length during each time step (Fig. 1). Consequently, a fixed number of growth time steps for a range of target fracture sizes means that fractures with a larger target size will grow faster than one with a smaller target size. This is a simplification of the growth rate used in some DFN growth models that include a growth-rate exponent parameter (Davy et al. 2010, 2013; Maillot et al. 2016).

Nucleation of new fractures continues until the target fracture intensity is reached. Individual fracture growth continues until the fracture reaches its target size, enough of its rays has terminated that growth on the whole fracture is arrested or the target fracture intensity is reached. As only one of these criteria is met when growth on the individual fracture finishes, the number of growth time steps is a target value.

As growth on individual fractures can arrest at an earlier time step than specified, it may take more time steps than specified to grow enough fractures to reach the target intensity. This means the total time steps value input by the user is also a target value. The number of new fractures nucleated in each time step is estimated from the number required to reach the target intensity given the target size distribution and the termination rate.

Shadow zones

The stress release on the flanks of growing fractures can inhibit the growth or nucleation of other fractures (Scholz et al. 1993; Cowie and Shipton 1998; Vermilye and Scholz 1998; Kim et al. 2004; Bonneau et al. 2013). Lavoine et al. (2020) describe a development of the rules-based fracture model in Davy et al. (2013) that preferentially nucleates fractures away from the flanks and towards the tips of existing fractures. However, this approach requires a matrix stress model that is a function of existing fractures.

The GDFN algorithm includes behaviour to inhibit fracture nucleation on the flanks of existing fractures without incurring the computational cost of modelling matrix stress. This is achieved with regions on either side of GDFN fractures that inhibit nucleation of new fractures of the same set. These zones, which are termed shadow zones, are represented as pyramidal polyhedrons, the base of which is centred on the nucleation point of the fracture (Fig. 2). The maximum height of the polyhedron is proportional to the equivalent radius of the fracture surface, and so increases as the fracture grows.

Figure 3 shows a trace map of a cross-section through an example GDFN with shadow zones, illustrating that the regions around the larger, older fractures have reduced intensity. Figure 4 shows a pair of trace maps that illustrate the regularization of fracture spacing due to shadow zones. These behaviours can be calibrated to observed fracture trace maps.

Fracture terminations

To mimic the cross-cutting relationships that may be observed in natural fractures, GDFN fractures can be configured to terminate on pre-existing fractures according to a specified probability. Termination is assumed to only occur against fractures larger than the growing fractures, a behaviour that is an oversimplification of the natural processes but one that can provide a reasonable approximation of the natural equivalent (Davy et al. 2010, 2013; Maillot et al. 2016; Lavoine et al. 2020).

As terminations may restrict fracture growth, the target fracture size may not be achieved. Consequently, the size distribution of a GDFN is an emergent parameter that is sensitive to the termination probability and target intensity.

Termination of a given fracture occurs on a ray-by-ray basis. This means that growth may cease on some rays whilst it continues on others. The decision to terminate a particular ray happens when a ray, or the growing edge between rays, intersects with the surface of another fracture. When this potential collision is detected, an initial check is made to see if a collision between the two fractures has already been detected (in this or a previous time step). If a collision has already occurred, the decision to terminate or not terminate in the first collision is honoured. If a collision has not already occurred, the decision to terminate or not is made using the user-input termination probability. If this probability calculation determines that termination will occur, the terminated rays grow the amount required so that the edge of the terminating fracture aligns with the plane of the fracture it is terminating against.

Once a given ray is terminated, the amount of growth that would be added to that ray in each subsequent time step is divided equally among the remaining non-terminated rays. Once a user-specifiable proportion of rays (i.e. a sufficient proportion of the fracture circumference) is terminated, growth on the whole fracture is arrested.

The greater the number of rays specified for each fracture, the higher the resolution of the intersection geometry.

Groups of GDFN fractures with common orientations can be generated in a sequence. Each group can be configured to terminate on pre-existing groups, or fractures in their own group, according to user-specified probabilities. This can be calibrated according to fracture patterns observed in outcrop and/or the deformation history of the modelled rock.

The primary aim of the GDFN model is to produce DFNs that better mimic the connectivity characteristics of natural fracture networks. To quantify how effectively this has been achieved the Sanderson and Nixon (2015) fracture topology metric may be used (Libby et al. 2019). This topology model has been linked to permeability and characterizes the network using metrics that are independent of scale and orientation (Sævik and Nixon 2017). It involves identifying the locations where traces terminate or cross other traces, subsequently termed nodes, and are classified as isolated (I), abutting (Y) or crossing (X), as shown in Figure 5.

By considering the proportions of the three different node types in a fracture network, the connectivity characteristics of a fracture network may be represented as a single point on a ternary plot, facilitating efficient comparison with other fracture networks. Figure 6 demonstrates how a wide variety of fracture patterns can be generated using the GDFN algorithm and how these span the full range of connectivity characteristics, in contrast to static models that do not provide the same controls on termination behaviour. The connectivity characteristics of a DFN trace map are influenced by every parameter input to the DFN that impacts its geometry. In addition, the sensitivity of the connectivity characteristics to some parameters is highly sensitive to the configuration of other parameters. For example, a high fracture intensity network may be very sensitive to termination probability and vice versa; however; either scenario may not be the case if the network is of an extreme size distribution or if the orientation distribution is such that it encourages or inhibits terminations. This complex interconnectivity of DFN parameters as controls on connectivity characteristics inhibit generalized statements on the calibration of these parameters. Consequently, proper calibration of the connectivity characteristics is achieved by iterative parameter adjustment supported by parameter sensitivity analyses. For information on the inputs to the models shown on Figure 6 and additional discussion of GDFN connectivity characteristics, please refer to Libby et al. (2019).

To illustrate some potential effects of the GDFN algorithm, a series of GDFNs were generated that included a single conjugate fracture set with the parameters shown in Table 1. GDFNs were generated for a range of intensities, detailed subsequently, to produce 50 realizations with total fracture counts of approximately 17 000–55 000.

The resulting DFNs were compared with an equivalent static DFN with the same properties but without terminations or shadow zones. To produce a static DFN with an identical size and orientation distribution, the GDFN fractures were stochastically relocated throughout the fractured region, as per a similar comparison in Maillot et al. (2016). As there are fewer fractures at the edge of the generation region for fractures to terminate against, the GDFN fractures at the edges of the region tend to be larger. To preclude this influence, fractures were generated in a region 10 m larger in all dimensions than the region used for the following analysis. Because the GDFN terminations modify the shape of the fractures, the stochastically relocated fractures were also restored to (approximate) circles. Comparative analyses were then run on these DFN pairs. Example trace maps of one of these DFN pairs are shown in Figure 7.

Here, ‘blocks’ are regions fully bounded by fractures on a 2D regular grid with 25 cm resolution (Elmo et al. 2014). For a given intensity, the GDFN produces a more uniform distribution of block sizes than the static DFN. This is illustrated qualitatively in Figure 7 for a single realization and demonstrated quantitatively for all 50 realizations in Figure 8. For example, when the P32 fracture intensity is 1.5 m2 m−3, the proportion of blocks with a cross-section greater than 10 m2 is c. 6% for the static DFN but c. 9% for the GDFN. A c. 50% increase in the proportion of large blocks would have implications for the selection and design of ground support strategies (e.g. rock bolting and shotcreting) and excavation sequences. The greater separation between the GDFN curves relative to the static DFN curves arises from the increasing number of terminations that occur as the intensity increases. The cross-cutting of the different curves for the GDFN arises from interaction between different aspects of the GDFN algorithm: for example, higher intensities increases the number of potential termination opportunities, resulting in fewer fractures reaching their target size, which creates smaller shadow zones, reducing their regularization effect. There are many such interactions, all of which in turn affect each other. Consequently, a full description of the precise cause of each behaviour shown in Figure 8 is beyond the scope of this paper.

To compare the permeability of the GDFN and static DFNs, their flow properties were upscaled to a regular grid with 5 m cells. The equivalent permeability of the portion of the DFN within each grid cell was calculated by meshing the subset of the DFN inside each cell and solving for steady-state flow by applying a linear head gradient across the cell (WSP UK Ltd 2023). To ensure that geometrical connectivity was the dominant factor controlling permeability, a constant transmissivity value of 10−6 m2 s−1 was assigned to all fractures. Figure 9 summarizes the resulting distribution of flows for 10 fracture realizations at each intensity. It may be observed that, in comparison to the static equivalent, the GDFN produces a flow distribution with a more uniform range of permeabilities and a higher average permeability. This has implications wherever DFNs are used to support the characterization of groundwater flow.

The more uniform range of permeabilities is due to the regularization of fracture spacing caused by the shadow zones. This decreases the likelihood of cells having a significantly higher or lower fracture count, leading to more consistent permeability values. The mean permeabilities are 5–10% higher for the GDFNs relative to the static DFNs. This is in contrast to an existing rules-based approach (Davy et al. 2013), where the permeability of the rules-based approach was 1.5–10 times smaller than the static equivalent (Maillot et al. 2016). The increase in average permeabilities is due to the fracture terminations reducing the area of ‘dead-end’ fracture area, causing more of the available fracture area to contribute to flow through better connectivity.

The flow channelling indicator detailed in Maillot et al. (2016) is a measure of the density of the parts of a fracture network that transport a significant part of flow, and is formulated as:
where f denotes the components of the network on which channelization is assessed, usually fractures but modified later in this paper; and dQ is the flow channelling density, V is the volume of the domain over which the flow channelling is calculated, Sf is the component surface area and Qf is the component flow.

Figure 10 shows how this flow channelling indicator changes with increasing intensity for the GDFN and static DFN pairs. The greater values for the GDFNs indicate reduced channelling relative to the static DFNs, with the disparity increasing as intensity increases. This is consistent with the more regularized fracture spacing of the GDFNs, the effect of which increases with intensity. This reduced channelling will increase the surface area of the matrix in contact with percolating fractures, potentially enhancing fracture–matrix interactions such as matrix diffusion, sorption and thermal conduction.

Although more fundamental properties of the DFN such as intensity remain the primary control on the overall connectivity of a DFN, these findings demonstrate that, given identical size, orientation and transmissivity distributions, the GDFN rules-based approach can produce a more regular and connected DFN than the static equivalent.

The reduced channelling and increased permeability of GDFNs relative to their static equivalent contrasts with the rules-based DFN generation model applied in Maillot et al. (2016) and originally described in Davy et al. (2013). This is despite the Davy et al. (2013) model using similar rules to the GDFN. A key difference between the GDFN and the Davy et al. (2013) model is the handling of terminations, and we suggest that this may be the primary cause of the differences in behaviour.

In the Davy et al. (2013) model, when termination occurs, continued fracture growth either ceases or moves the centre of the growing fracture away from the one it has terminated on. This means that the length of the terminated edge does not increase. In contrast, when GDFN rays terminate, fracture growth can continue by extension of the unterminated rays (Libby et al. 2019). Consequently, as more rays terminate, the length of the terminated edge increases. Figure 11 illustrates two end members of this behaviour. Figure 11a shows a fracture with a small proportion of rays terminated and the resulting short length (relative to the size of the fracture) of terminated fracture edge. Figure 11b shows a fracture with a large proportion of rays terminated and the resulting long length (relative to the size of the fracture) of terminated fracture edge. This is significant because the length of the terminated edge determines the size of the pathway available to flow between the two fractures.

To test the influence of this behaviour, a series of GDFN realizations were generated with random orientations, a constant target size of 10 m, no shadow zone and a termination probability of 100%. The random orientations of these GDFNs match the similar configuration used for the channelling assessment of Davy et al. (2013). The DFNs were generated to a P32 intensity of 0.2 m2 m−3 in a domain measuring 200 × 100 × 100 m. As for the DFNs described earlier in this paper, the static equivalent of the GDFNs was created by stochastically locating near-circular fractures with orientations and sizes identical to the GDFNs. By varying the proportion of rays that can be terminated before fracture growth ceased on the GDFN fractures, it is possible to quantify the impact of the termination behaviour. The result is shown in Figure 12.

The main cause of the positive trends in the flow channelling density of both the grown and static DFNs is an increase in fracture size as the termination percentage increases. This occurs as the lower termination probabilities result in earlier termination of growth. By comparing the GDFN to the static equivalent, these changes in fracture size can be accounted for. It may be observed that reducing the proportion of rays that can be terminated before growth ceases on the GDFN fractures results in a less channelized DFN, with termination percentages of less than 15% resulting in GDFNs more channelized than the static equivalent. Consequently, we speculate that the termination behaviour of the Davy et al. (2013) model causes flow within the DFNs it produces to become more channelized.

We favour a large (e.g. >40%) proportion of rays able to be terminated before growth ceases on GDFN fractures. As exposed bedrock typically represents a mostly 2D view of the 3D natural fracture network, we argue that the large number of terminations that may be observed in some natural fracture networks (Fig. 6) demonstrates that the fracture configuration illustrated in Figure 11b occurs more frequently in nature than that in Figure 11a. This observation is backed up by consideration of the mechanical forces: in the scenario illustrated in Figure 11a there are narrow sections of unfractured rock either side of the terminated edge that would be mechanically weak and therefore preferentially broken during ongoing fracture growth, resulting in something more like what is illustrated in Figure 11b. Simulation of these scenarios with 3D full physics models could offer more insight.

Figure 13 illustrates the flow channelling indicator calculated for the same DFN pairs on a per-flow-mesh-cell basis as opposed to a per-fracture basis. In comparison to Figure 12, a significantly larger increase in flow channelling for the GDFN relative to the static DFN may be observed. As observations of channelling behaviour in natural fracture networks will commonly be limited to borehole and tunnel data that bisect fractures, we suggest that quantifying channelling with a method that subdivides individual fractures may allow better comparisons with this observed data. Further work is required to explore how channelling in rules-based DFNs may be best quantified.

The stochastic aspects of DFNs mean that their effective and scientifically rigorous application typically requires probabilistic assessments involving multiple iterations of fracture generation. This compounds with the parametric sensitivity analyses often required for proper model calibration. Consequently, computational efficiency is a primary component of what differentiates the GDFN algorithm from conceptually similar alternatives and makes it suitable for technical applications. To achieve this computational efficiency, it employs a multi-threaded approach using localized interactions rather than solving global systems. For illustration, some indicative timings are given below, recognizing that many factors, including both hardware and parameter configuration, influence the generation time:

  • Each GDFN realization described in the ‘Termination behaviour’ section contains 40 000–140 000 fractures that were generated on a six-core laptop in 30–120 s, with the range of fracture count and computation time dependent on the specified fracture intensity.

  • A GDFN with 7 million fractures each with full 3D geometries, individual properties, terminations and shadow zones can be generated in 25 min on a six-core laptop, whilst using a maximum of 11 GB of RAM.

This enables probabilistic assessments to be undertaken using large GDFNs.

Libby et al. (2019) demonstrated that an appropriate configuration of the GDFN model can produce DFNs that can honour the full range of connectivity characteristics that may appear in natural fracture networks (Fig. 6). An appropriately configured GDFN can therefore more closely honour observed natural fracture patterns, giving greater reliability in predicting aspects of fracture networks that are affected by these connectivity characteristics, such as rock-mass stability and deformation, hydraulic responses, and fluid flow (Sævik and Nixon 2017).

This paper demonstrates that these connectivity characteristics have implications for the mechanical behaviours of the rock hosting the fracture network and the flow properties of the network. Figures 8 and 9 illustrate how, for a given fracture intensity relative to an equivalent DFN generated using established static algorithms, the GDFN algorithm can increase the hydraulic and mechanical connectivity of the generated DFN and decrease flow channelling. This is mostly achieved through the regularization of fracture spacing and a reduction in the dead-end (hydraulically stagnant) fracture area.

These changes also increase the surface area of the matrix in contact with percolating fractures, potentially enhancing fracture–matrix interactions such as matrix diffusion, sorption and thermal conduction. The regularized and increased connectivity also results in an increased large-block count and a corresponding reduction in the average block volume. This can be critical for predicted rock-mass behaviour and ground support response in the application of continuum or discontinuum numerical tools in geotechnical modelling. This has implications for the selection and design of ground support strategies (e.g. rock bolting and shotcreting), excavation sequences, and estimates from rock-mass classification systems (Palmstrom 2005). DFNs have also supported developments of standard rock-mass metrics such as rock quality designation and geological strength index (Marinos and Carter 2018; Schlotfeldt and Carter 2018). Further work is required to assess how rules-based DFNs such as GDFNs may be used to inform or expand this work.

The decreased channelization of fluid flow through DFNs generated with a rules-based approach in comparison to an equivalent static DFN is in contrast to the findings of Maillot et al. (2016). We demonstrate flow channelling is sensitive to the termination behaviour of the GDFNs and hypothesize that the findings of Maillot et al. (2016) could be explained by differences in how terminations are handled. In addition, the channelization metric is shown to be sensitive to the scale on which channelization is considered. Further work is required to explore both the impact of fracture orientation distribution on channelling, and the method of characterizing channelization.

We provide examples to illustrate that despite inclusion of mechanistic rules, the GDFN algorithm is computationally efficient enough that probabilistic assessments of large (>1 million fractures) GDFNs may be carried out in practical timescales on commonly available hardware. This supports many of the practical applications of DFNs.

We thank Martin Stigsson (Svensk Kärnbränslehantering AB) for his detailed review of the manuscript.

SL: conceptualization (supporting), data curation (lead), formal analysis (lead), investigation (lead), methodology (lead), validation (lead), visualization (lead), writing – original draft (lead), writing – review & editing (supporting); LH: conceptualization (equal), funding acquisition (equal), methodology (supporting), project administration (lead), software (supporting), supervision (lead), validation (supporting), writing – review & editing (lead); RT: methodology (supporting), software (lead), writing – review & editing (supporting); MC: conceptualization (supporting), funding acquisition (supporting), methodology (supporting), software (supporting), writing – review & editing (supporting); TB: software (supporting), validation (supporting); NJ: software (equal); RM: conceptualization (equal), formal analysis (supporting), methodology (supporting), writing – review & editing (equal); J-OS: conceptualization (supporting), funding acquisition (lead), writing – review & editing (supporting); DMI: writing – review & editing (supporting).

This work was financially supported by Svensk Kärnbränslehantering AB, the Swedish Nuclear Fuel and Waste Management Company, and Posiva Oy, the Finnish radioactive waste management company.

The authors declare that they have no known competing personal relationships that could have appeared to influence the work reported in this paper; financial interests are declared under the funding section.

The datasets generated during and/or analysed during the current study are not publicly available due to their large size but are available from the corresponding author on reasonable request.

This is an Open Access article distributed under the terms of the Creative Commons Attribution 4.0 License (http://creativecommons.org/licenses/by/4.0/)