Investigating the effects of intersection flow localization in equivalent-continuum-based upscaling of flow in discrete fracture networks

Predicting effective permeabilities of fractured rock masses is a crucial component of reservoir modeling. Its often realized with the discrete fracture network (DFN) method, whereby single-phase incompressible fluid flow is modeled in discrete representations of individual fractures in a network. Depending on the overall number of fractures, this can result in high computational costs. Equivalent continuum models (ECMs) provide an alternative approach by subdividing the fracture network into a grid of continuous medium cells, over which hydraulic properties are averaged for fluid flow simulations. While continuum methods have the advantage of lower computational costs and the possibility of including matrix properties, choosing the right cell size to discretize the fracture network into an ECM is crucial to provide accurate flow results and conserve anisotropic flow properties. Whereas several techniques exist to map a fracture network onto a grid of continuum cells, the complexity related to flow in fracture intersections is often ignored. Here, numerical simulations of Stokes flow in simple fracture intersections are utilized to analyze their effect on permeability. It is demonstrated that intersection lineaments oriented parallel to the principal direction of flow increase permeability in a process we term intersection flow localization (IFL). We propose a new method to generate ECMs that includes this effect with a directional pipe flow parameterization: the fracture-and-pipe model. Our approach is compared against an ECM method that does not take IFL into account by performing ECM-based upscaling with a massively parallelized Darcy flow solver capable of representing permeability anisotropy for individual grid cells. While IFL results in an increase in permeability at the local scale of the ECM cell (fracture scale), its effects on network-scale flow are minor. We investigated the effects of IFL for test cases with orthogonal fracture formations for various scales, fracture lengths, hydraulic apertures, and fracture densities. Only for global fracture porosities above 30 % does IFL start to increase the systems permeability. For lower fracture densities, the effects of IFL are smeared out in the upscaling process. However, we noticed a strong dependency of ECM-based upscaling on its grid resolution. Resolution tests suggests that, as long as the cell size is smaller than the minimal fracture length and larger than the maximal hydraulic aperture of the considered fracture network, the resulting effective permeabilities and anisotropies are resolution-independent. Within that range, ECMs are applicable to upscale flow in fracture networks.

1 Introduction

Discontinuities in rocks provide major pathways for subsurface fluid migration. Thus, fractured reservoirs are frequent targets for oil, gas, or water production, geothermal energy recovery, and CO₂ sequestration. In addition, the safety of nuclear waste disposal and subsurface contaminant transport crucially depends on the presence of fractures. Characterizing natural fracture networks across scales and modeling the fluid flow therein to predict their effective permeabilities has thus been a long-standing topic of research (e.g., Long et al., 1982; Dershowitz and Einstein, 1988; Cacas et al., 1990; Neuman, 2005; de Dreuzy et al., 2012).

Numerical modeling of fluid flow is most accurately based on the Navier–Stokes equations (Bear, 1972). For a single phase of incompressible and isoviscous fluid in an isothermal system, they simplify to the Stokes equations if laminar flow conditions are considered (i.e., Reynolds numbers below 1–10). Assuming an impermeable rock matrix, one can solve for the velocity distribution resulting from prescribed pressure boundary conditions, allowing for the determination of the rock effective permeability utilizing Darcy's law for flow through porous media (e.g., Andrä et al., 2013b; Osorno et al., 2015; Eichheimer et al., 2019, 2020; Kottwitz et al., 2020). Those so-called direct flow modeling approaches crucially rely on a digital representation of a rock that separates pore space from the matrix, which results from high-resolution X-ray computed tomographies (Andrä et al., 2013a; Cnudde and Boone, 2013). However, they are limited in maximum scannable size and respective trade-off to numerical resolution, making them applicable to small scales only (nanometers to a couple of centimeters at most). At larger scales (above a couple of centimeters), so-called continuum flow approaches serve to model fluid flow, usually based on the concepts for flow through porous media proposed by Darcy (Darcy, 1856). Instead of a representation of the medium's pore space, they require an initial hydraulic representation of the medium. This is given by prescribed effective permeabilities for certain control volumes within the medium, which upscale hydraulic properties from smaller scales to observation scales. Thus, the key of this so-called upscaling problem (e.g., Zhou et al., 2010; Hauge et al., 2012; Lie, 2019) is to adequately represent the rock structure with an appropriate model of effective permeabilities, which for fractured rock masses is often cumbersome due to their structural heterogeneity (Dershowitz and Einstein, 1988; Odling et al., 1999). The main problem is that acquiring detailed natural fracture data in 3D is intricate, as seismic imaging techniques suffer from resolution limits (Cartwright and Huuse, 2005; Malehmir et al., 2017), preventing a multi-scale structural assessment of individual features in fracture formations. Hence, outcrop (2D) and borehole (1D) studies are the only possibilities to acquire detailed natural fracture data, despite their reduced dimensionality (Lei et al., 2017), and acquiring deterministic knowledge of all individual structures in a fracture formation is impossible. Due to this, the discrete fracture network (DFN) method has been extensively used as a conceptual framework to provide statistically based approximations of real fracture networks for decades (Long et al., 1982; Cacas et al., 1990; Bogdanov et al., 2003; Darcel et al., 2003; Xu and Dowd, 2010; Davy et al., 2013; Maillot et al., 2016). In this approach, each fracture in a given network is represented by a reduced-order object (lines in 2D and disks or rectangles in 3D) with a prescribed location, size, and orientation. Naturally measured structural properties like size and orientation distributions (Odling et al., 1999; Healy et al., 2017) as well as fracture density and spacing (Ortega et al., 2006) serve as a quantitative basis to prescribe their geometrical properties (e.g., Hyman et al., 2015; Alghalandis, 2017). The hydraulic response to pressure changes in each individual fracture is then parameterized with the cubic law (Snow, 1969; Witherspoon et al., 1980), relating the fracture's effective permeability to its aperture. In reality, surface roughness, fracture closure, and fluid–rock interactions (e.g., erosion or crystal growth) cause deviations from the parallel-plate assumption (Brown, 1995; Oron and Berkowitz, 1998; Méheust and Schmittbuhl, 2000). Semi-empirical functions derived from numerical simulations in rough-walled fractures with quantified statistics of the aperture field (e.g., Patir and Cheng, 1978; Brown, 1987; Renshaw, 1995; Zimmerman and Bodvarsson, 1996; Mourzenko et al., 2018) serve as corrections to the cubic law if the fracture's internal correlation length scale is significantly smaller than the size of the considered fracture (e.g., Méheust and Schmittbuhl, 2003; Kottwitz et al., 2020).

A large number of numerical methods to compute effective permeabilities of fractured media have been developed (see reviews of Jing, 2003; Berre et al., 2019), all relying on (modified) cubic-law assumptions. Improved discretization techniques with individual fracture treatment (DFN method), inclusion of matrix properties in multi-dimensional meshes (discrete fracture and matrix – DFM – method), and multi-continuum methods come at the cost of high computational expenses. Discretizing the fractured media as equivalent single-continuum blocks significantly reduces the computational effort at comparable numerical accuracy (Hadgu et al., 2017).

According to Long et al. (1982) and Oda (1985), fractured rocks behave similarly to porous media and can be represented by a positive definite permeability tensor (Chen et al., 1999) as long as the considered system behaves like a representative elementary volume (REV) (Bear, 1972); i.e., its effective properties (permeability or porosity for example) are more or less homogenous at the reference scale of the system. Due to the multi-scale character of fracture systems (e.g., Bonnet et al., 2001; Davy et al., 2006), determining the required homogenization scale is difficult, as distinct larger features may dominate overall flow. Thus, a discrete representation of all fractures in a network given by the DFN method is essential to adequately capture that multi-scale character. However, La Pointe et al. (1995), Jackson et al. (2000), Svensson (2001), Leung et al. (2012), and Hadgu et al. (2017), among others, have shown that representing a DFN with a grid of equivalent continuum blocks of sizes lower than the REV yields similar flow results, if resolved sufficiently, and thus reproduces the overall flow behavior of the DFN method. This highlights the fact that continuum methods for flow modeling in fractured rocks are not restricted to REV scales and can thus be used equivalently to the DFN method.

Several techniques to generate equivalent continuum models (ECMs) of DFNs have been developed in 2D (Reeves et al., 2008; Botros et al., 2008; Rutqvist et al., 2013; Chen et al., 2015) and 3D (Hadgu et al., 2017; Sweeney et al., 2020), whereby the so-called Oda method (see Oda, 1985) is used to formulate permeability tensors of grid cells that intersect fractures. The permeability tensor is aligned with the orientation of the intersecting fracture, and the permeabilities of the individual fractures are summed up if multiple fractures intersect one cell, yielding a positive definite, fully anisotropic tensor (e.g., Chen et al., 1999). The groundwater flow equations for porous media (Bear, 1972), i.e., Darcy's law (Darcy, 1856), are then used to simulate laminar, steady-state, single-phase flow to compute effective permeabilities of the medium. A current issue in commonly used 3D flow solvers, such as PFLOTRAN (Lichtner et al., 2015), is a lack of a fully anisotropic permeability representation at the local cell level. So-called staircase patterns are the direct consequence of these simplifications, which introduce artificially prolonged flow paths, especially in transport simulations, which have to be compensated for (e.g., Reeves et al., 2008; Botros et al., 2008; Sweeney et al., 2020) when predicting effective permeabilities of fractured media. On the other hand, MODFLOW (McDonald and Harbaugh, 1988) introduced support for local permeability anisotropy but not within a massively parallelized framework, making it difficult to conduct large numbers of high-resolution simulations. However, assessing permeabilities in a Monte-Carlo-like framework (e.g., de Dreuzy et al., 2012) is necessary to explore the variance of hydraulic system properties induced by stochastically generated input data. Hence, a flow solver that combines the advantages of local permeability anisotropy and massive parallelization should be beneficial for numerical permeability assessments of fracture networks.

Next to these issues, this study focuses on an often ignored but potentially important aspect in fracture network modeling given by the complexity of fracture intersection flow. To our knowledge, only a few studies have presented 3D flow simulations within fracture intersections (Zou et al., 2017; Li et al., 2020), revealing the fact that flow velocities will increase within the fracture intersections compared to the fracture itself (shown by increasing Péclet numbers within the intersections). Theoretically, this effect should increase if the direction of the applied pressure gradient is aligned with the orientation of the intersection. As a consequence, the system's effective permeability should increase by a certain amount due to a local permeability increase within the intersection. To demonstrate that, we systematically conduct 3D numerical simulations of Stokes flow within differently oriented, planar fracture crossings to analyze the permeability increase caused by intersection flow localization (IFL). Using these results, we extend the current state-of-the-art methodology for equivalent continuum representations of DFNs to account for IFL in a quantitative manner and analyze its impact on effective permeability computations at fracture and network scales. It is still unclear at which level of detail the ECM has to be discretized to conserve the structural complexity of the DFN, as the aforementioned staircase patterns and artificial connectivity cause resolution dependencies. Subsequently, resolution tests are performed on two DFN test cases with a newly developed, massively parallelized, and high-performance computing (HPC) optimized finite-element Darcy flow solver capable of handling fully anisotropic permeability tensor cells. By that, we consistently investigate the upscaling capabilities of the ECM method, which is frequently used for effective permeability predictions in fractured porous media.

Figure 1 Panels (a, d) show the binary voxel models (impermeable matrix in transparent gray) for a fracture intersection that is oriented along and transverse to the flow direction, respectively. The red bottom face is the high-pressure boundary (0.02 Pa) and the blue top face the low-pressure boundary (0.01 Pa), forcing the fluid to flow in the z direction. The orientations (arranged as dip direction $/$ dip) for the fracture pair in (a) are $f_{\mathrm{1}}=\mathrm{100}/\mathrm{90}$, $f_{\mathrm{2}}=\mathrm{190}/\mathrm{80}$ and $f_{\mathrm{1}}=\mathrm{170}/\mathrm{90}$, $f_{\mathrm{2}}=\mathrm{260}/\mathrm{10}$ for the fractures in (d). The length of both cubes is 1 cm, and all fracture apertures are constant (1.25 mm). Panels (b, e) visualize flow velocity distribution in the void space. Panels (e, f) highlight velocity vectors within the intersections at slices indicated with green rectangles in (b, e), respectively.

Figure 2 2D sketch of the half-hypotenuse assumption in an idealized rectangular fracture crossing (gray regions indicate rock matrix, white regions fracture pore space). The hydraulic apertures (a_h1 and a_h2) of both intersecting fractures are indicated with solid blue lines. The hypotenuse of the right-angled triangle with the two hydraulic apertures as legs is given by the black dashed line. The hydraulic radius r_h (indicated by the solid red line) to approximate the radius of the pipe model is defined as half the length of the hypotenuse.

2 Fracture intersection flow modeling

Fluid flow in porous and fractured media is described by the well-known Navier–Stokes equations (Bear, 1972). It is commonly assumed that subsurface flow in fractures ranges in the laminar regime, i.e., Reynolds numbers below unity (Zimmerman and Bodvarsson, 1996). Assuming the flowing fluid to be incompressible and isoviscous and the impact of gravity to be negligible, steady-state flow at constant temperature is defined by the Stokes momentum balance (Eq. 1) and continuity (Eq. 2) equations (Bear, 1972):

$$\begin{array}{}\text{(1)}& {\displaystyle }& {\displaystyle }\mathit{\mu }{\mathrm{\nabla }}^{\mathrm{2}}\mathit{v}=\mathrm{\nabla }P,\text{(2)}& {\displaystyle }& {\displaystyle }\mathrm{\nabla }\cdot \mathit{v}=\mathrm{0},\end{array}$$

with the fluid's dynamic viscosity μ, pressure P, and velocity vector $\mathit{v}=(v_x,v_y,v_z)$. ∇, ∇⋅, and ∇² denote the gradient, divergence, and Laplace operator for 3D Cartesian coordinates, respectively.

Here, the 3D staggered-grid, finite-difference code LaMEM (Kaus et al., 2016) is used to solve the coupled system of Eqs. (1) and (2) utilizing PETSc (Balay et al., 2018) for HPC optimization. Applying different absolute pressures on two opposing sides of a 3D voxel model representing the fractured or porous medium (e.g., a or d in Fig. 1), while setting the other boundaries to no-slip (velocity component normal to the boundary is zero), enables the prediction of the medium's directional permeability. After obtaining the steady-state solution, the volume integral of the pressure-gradient-aligned velocity component v_z (e.g., Osorno et al., 2015) is computed according to

$$\begin{array}{}\text{(3)}& {\displaystyle }{\displaystyle }\overline{v}={\displaystyle \frac{\mathrm{1}}{V}}\underset{V}{\int }\left|v_z\right|\mathrm{d}z,\end{array}$$

with domain volume V. Using Darcy's law for flow through porous media (Darcy, 1856), which relates the specific discharge Q for a pressure drop ΔP along a distance L according to

$$\begin{array}{}\text{(4)}& {\displaystyle }{\displaystyle }Q=-{\displaystyle \frac{kA\mathrm{\Delta }P}{\mathit{\mu }L}},\end{array}$$

with intrinsic permeability k and cross-sectional area A in combination with the fact that $Q=\overline{v}A$, the directional permeability k_z is calculated by

$$\begin{array}{}\text{(5)}& {\displaystyle }{\displaystyle }{k}_z={\displaystyle \frac{\mathit{\mu }\overline{v}L}{\mathrm{\Delta }P}}.\end{array}$$

As demonstrated by Eichheimer et al. (2019), Kottwitz et al. (2020), and Eichheimer et al. (2020), the numerical resolution has to be sufficiently high to produce a converged result. Generating every model at different levels of detail (e.g., 128³, 256³, 512³, and 1024³ voxels) ensures that the most accurate solution is obtained (as will be shown later by a comparison of errors to the result at the largest resolution in Fig. 5b). Figure 1 presents Stokes flow in simple fracture intersections and highlights the IFL effect. If the fracture intersection is aligned with the principal flow direction (panels a–c), the velocity significantly increases within the intersection, resulting in higher directional permeabilities. In the opposite case, when the fracture intersection connects no-pressure boundaries (panels d–f) and is thus oriented oblique to the flow direction, the flow velocity slightly disperses around the intersection, and the overall impact on the directional permeability is minor.

3 Permeability parameterization concepts

As the two main structural features (fractures and intersections) composing a fracture network differ significantly in terms of their hydraulics (Fig. 1), they require independent concepts to parameterize their permeabilities for formulating their effective grid block permeability tensor. For fractures, it is usual practice to use the cubic-law parameterization (e.g., Snow, 1969; Long et al., 1982), relating the specific discharge Q through a void system between two parallel plates for a pressure drop ΔP along a distance L according to

$$\begin{array}{}\text{(6)}& {\displaystyle }{\displaystyle }Q=-{\displaystyle \frac{w{a}_{\mathrm{m}}^{\mathrm{3}}\mathrm{\Delta }P}{\mathrm{12}\mathit{\mu }L}},\end{array}$$

with the fracture width w and distance between the two plates, i.e., mechanical aperture, a_m. Comparing this analytical solution with Darcy's law (Eq. 4, cross-sectional area A=wa_m) leaves the intrinsic permeability of a fracture k_f defined by

$$\begin{array}{}\text{(7)}& {\displaystyle }{\displaystyle }{k}_{\mathrm{f}}={\displaystyle \frac{a_{\mathrm{m}}^{\mathrm{2}}}{\mathrm{12}}}.\end{array}$$

Natural fractures deviate from the assumptions of parallel plates, which is why a_m in Eq. (7) is commonly replaced with a hydraulic aperture a_h that corrects the parameterization for fracture closure and surface roughness (e.g., Patir and Cheng, 1978; Brown, 1987; Renshaw, 1995; Zimmerman and Bodvarsson, 1996; Kottwitz et al., 2020). Yet, there is no ready-to-use parameterization concept tailored for fracture intersections. The simulations shown in Fig. 1 suggest that the flow in the intersection is approximately pipe-like. Then, the specific discharge Q through a tube of radius r_t and length L is related by the Hagen–Poiseuille flow solution through a pipe (e.g., Batchelor, 1967) according to

$$\begin{array}{}\text{(8)}& {\displaystyle }{\displaystyle }Q=-{\displaystyle \frac{\mathit{\pi }{r}_{\mathrm{t}}^{\mathrm{4}}\mathrm{\Delta }P}{\mathrm{8}L\mathit{\mu }}}.\end{array}$$

Again, combining this equation with Darcy's law (Eq. 4, cross-sectional area A=πr²) results in the following expression for the intrinsic permeability of a pipe k_i:

$$\begin{array}{}\text{(9)}& {\displaystyle }{\displaystyle }{k}_{\mathrm{i}}={\displaystyle \frac{r^{\mathrm{2}}}{\mathrm{8}}}.\end{array}$$

The apparent pipe radius should then be modified based on the intersection shape to calculate an equivalent hydraulic radius r_h to compensate for the structural difference. As a first-order approximation, we use half the size of the hypotenuse in a right-angled triangle whose legs are given by the two intersecting apertures (called half-hypotenuse assumption in the following; see Fig. 2 for details). This delivers sufficiently good results, as will be demonstrated later (Fig. 6).

Figure 3 Workflow for generating an equivalent continuum model of a DFN. Panel (a) shows the input fracture network of four arbitrarily oriented fractures (gray) and their intersections (magenta). Panel (b) displays a grid of ellipsoids, each reflecting the shape of the permeability tensor in the equivalent continuum model in (a) with a resolution of 4³ voxels. The size of the ellipses is scaled to the norm of the permeability tensor of the cell such that larger ellipsoids denote higher permeabilities. The green plane in (b) indicates the location of the 2D slice displayed in (c). Different green intensities present the norm of the permeability tensor of each cell. Black lines denote fractures in 2D and yellow ellipses the x and y shape of the permeability tensor of each cell. Note how the shape of the ellipse changes from being planar if multiple fractures cross a cell.

Figure 4 Fracture-intersection-caused changes in permeability tensor characteristics. Panel (a) shows a simple DFN structure of two arbitrarily oriented fractures (gray) intersecting at a line (magenta). The cube length is set to 1 m and the system origin is at $(\mathrm{0},\mathrm{0},\mathrm{0})$. The center point of the first fracture is located at $\left(\mathrm{0.4899}\right|\mathrm{0.5685}\left|\mathrm{0.5110}\right)$ and its normal vector is given as $(-\mathrm{0.3195},\mathrm{0.7894},\mathrm{0.524})$. The second fracture's center point is located at $\left(\mathrm{0.7604}\right|\mathrm{0.5000}\left|\mathrm{0.5000}\right)$, whereas its normal vector is given by $(-\mathrm{0.9461},\mathrm{0.1715},\mathrm{0.2747})$. Both fractures have the same hydraulic aperture of $\mathrm{1}\times {\mathrm{10}}^{-\mathrm{3}}\mathrm{m}$ and both fully penetrate the system. The resulting intersection ranges from point $\left(\mathrm{0.6499}\right|\mathrm{0.3086}\left|\mathrm{1.000}\right)$ to $\left(\mathrm{0.8003}\right|\mathrm{1.000}\left|\mathrm{0.0505}\right)$, and its orientation is given by the unit vector $(\mathrm{0.1270},\mathrm{0.5839},-\mathrm{0.8018})$. The hydraulic pipe radius resulting from the half-hypotenuse assumption is $\mathrm{7.0711}\times {\mathrm{10}}^{-\mathrm{4}}$. Panel (b) visualizes the shape of the permeability tensor for an ECM model that considers only fracture permeability (gray, inside) and for the presented fracture-and-pipe model (transparent magenta, outside). The size of both ellipses is scaled with the norm of the resulting permeability tensor to provide comparability. Panel (c) presents the norm of the permeability tensor K_ijk as a function of the ratio between the ECM grid spacing δx and the initial point spacing δp for the discretization approach described in Sect. 4 (the fracture-and-pipe model) and an approach wherein we did not take the IFL parameterization into account (i.e., leaving out Eq. 14 in the discretization procedure, hence the name fracture-only). The dashed black line denotes the condition $\mathit{\delta }p/\mathit{\delta }x\ge \mathrm{16}$, which is used to provide a correct approximation of the fracture and intersection volume per cell.

Figure 5 Panel (a) displays the location of all 100 intersection lineaments considered in the flow benchmark. A total of 52 intersection configurations directly connect inlets and outlets of flow (upper and lower z face), whereas 48 connect non-boundary flow faces. Panel (b) compares the numerically estimated permeability at highest resolution (1024³ voxels) to the ones obtained at lower resolutions by calculating their error norms $\left|\right|{\mathit{\delta }}_k\left|\right|$ according to Eq. (15). Gray dots represent the average error norm for all considered intersection configurations at resolutions lower than 1024³ voxels, and the light gray area highlights the range between minimum and maximum error.

Figure 6 The left panel shows a comparison of directional permeabilities obtained from high-resolution Stokes flow simulations (k_s) and counterparts (k_ecm) analytically derived with the ECM approach described in the text as a function of the angle γ between the intersection and the principal flow direction. Magenta dots represent the mean permeability ratios (10 values per point) for the ECM approach described in Sect. 4 with the half-hypotenuse pipe radius parameterization. Gray dots present the mean permeability ratio for an ECM approach that ignores the effect of intersections. The right panel shows a box plot of the error norm $\left|\right|{\mathit{\delta }}_k\left|\right|$ computed according to Eq. (15) with k_ecm as k_r for all 100 fracture-and-pipe (magenta) and fracture-only models (gray).

4 Equivalent continuum representation of DFNs

The use of the ECM approach instead of the DFN method to predict the effective permeabilities of fractured media crucially depends on the capability to reflect the anisotropic flow properties at the scale of the continuum cells. Therefore, it is essential to integrate the geometry of a DFN into the generation procedure of the ECM instead of generating the grid cell conductivities in a stochastic manner (Hadgu et al., 2017). The accuracy of the ECM permeability prediction then depends on the resolution of the DFN-mapped continuum grid. Jackson et al. (2000) and Svensson (2001) already demonstrated that using cell sizes that are larger than the average fracture spacing of the network introduces artificial connectivity and hence overestimates effective permeabilities. Sufficient resolution of the continuum grid is therefore required to obtain comparable results with the DFN method (e.g., Botros et al., 2008; Leung et al., 2012).

To our knowledge, there is no approach to generate an ECM of a DFN that takes the effect of IFL (Fig. 1) into account. Thus, we will explain our new approach to generate continuum representations based on DFN structures – the fracture-and-pipe model.

Generally, the DFN approach offers a straightforward way to characterize structurally complex fracture networks. Most commonly, every fracture is modeled as a geometric primitive (here a disk) with a prescribed length l, center coordinate p₀, and unit normal vector $\overline{\mathit{n}}$ defining its orientation. Based on this, fracture intersections can be calculated to define the backbone of the network. Here, fracture intersections are approximated with a line defined by two points i₀ and i₁, whereas the unit vector $\overline{\mathit{i}}$ between the two points defines its orientation. The goal of the ECM method is to generate a 3D representation of computational cells that contains a symmetric, positive definite permeability tensor that is based on the fractures and their intersections. For simplicity, we prescribe a regular grid with constant x–y–z spacing δx instead of octree refined grids as utilized by Sweeney et al. (2020).

To map each individual fracture to its corresponding grid cells, we first assume a horizontal disk (normal vector $\overline{\mathit{g}}=[\mathrm{0},\mathrm{0},\mathrm{1}]$) at center point $p_g=(\mathrm{0},\mathrm{0},\mathrm{0})$ with corresponding fracture radius r ($r=l/\mathrm{2}$) and represent it with an equally spaced set of points in the x–y plane P_g, with the condition $\left|\right|{\mathbf{P}}_{\mathbf{g}}-p_g\left)\right||\le r$. By that, we obtain a constantly spaced grid of points representing the fracture in horizontal orientation, provided that the initial equal spacing of the points δp is a small fraction of the cell size δx to prevent gaps in the mapped 3D grid. Next, we seek the rotation matrix R_f that aligns the current normal vector of the x–y plane $\overline{\mathit{g}}=[\mathrm{0},\mathrm{0},\mathrm{1}]$ with the actual normal vector of the fracture $\overline{\mathit{n}}$. Utilizing the Rodrigues rotation formula (Rodrigues, 1840) around the rotation axis $\mathit{w}=(\overline{\mathit{g}}\times \overline{\mathit{n}})/\left|\right|(\overline{\mathit{g}}\times \overline{\mathit{n}})\left|\right|$ (unit vector orthogonal to $\overline{\mathit{g}}$ and $\overline{\mathit{n}}$) yields the rotation matrix R_f according to

$$\begin{array}{}\text{(10)}& {\displaystyle }{\displaystyle }{\mathbf{R}}_{\mathbf{f}}=\mathbf{I}+\left|\right|\overline{\mathit{g}}\times \overline{\mathit{n}}\left|\right|\mathbf{C}+(\mathrm{1}-\overline{\mathit{g}}\cdot \overline{\mathit{n}}){\mathbf{C}}^{\mathrm{2}},\end{array}$$

with ×, ⋅, and $\left|\right|\mathit{x}\left|\right|$ denoting the cross-product, dot product, and vector norm of x, respectively. I represents the 3 × 3 identity matrix and C the cross-product matrix of the rotation axis $\mathit{w}=[w_x,w_y,w_z]$:

$$\begin{array}{}\text{(11)}& {\displaystyle }{\displaystyle }\mathbf{C}=\left[\begin{array}{ccc}\mathrm{0}& w_z& w_y\\ w_z& \mathrm{0}& -w_x\\ -w_y& w_x& \mathrm{0}\end{array}\right].\end{array}$$

Following this, R_f is used to rotate the 3×n array of points representing the fracture plane P_g (n is the number of 3D points in P_g) around p_g and translate all points to the actual center point p₀ to produce a rotated set of points P_r representing the fracture in its actual 3D position:

$$\begin{array}{}\text{(12)}& {\displaystyle }{\displaystyle }{\mathbf{P}}_{\mathbf{r}}={\mathbf{P}}_{\mathbf{g}}*{\mathbf{R}}_{\mathbf{f}}+p_{\mathrm{0}},\end{array}$$

where * denotes matrix–matrix multiplication. By ensuring that the lower left corner coordinate of the rectangular grid's bounding box is initially located at $(\mathrm{0},\mathrm{0},\mathrm{0})$ (this may require a translation of all center points to incorporate all fractures), we obtain the grid indices (i, j, and k in the x, y, and z direction, respectively) of the fracture by dividing P_r with the cell size δx and rounding the results. Finally, we compute the individual anisotropic permeability tensor K_ijk for the cells by using a parameterized fracture permeability value (Eq. 7) and the rotation matrix R_f according to

$$\begin{array}{}\text{(13)}& {\displaystyle }{\displaystyle }{\mathbf{K}}_{ijk}={\displaystyle \frac{V_{\mathrm{f}}}{V_{\mathrm{c}}}}{k}_{\mathrm{f}}\left({\mathbf{R}}_{\mathbf{f}}\left[\begin{array}{ccc}\mathrm{1}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{1}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}\end{array}\right]{\mathbf{R}}_{\mathbf{f}}^{\prime }\right).\end{array}$$

V_c denotes the cell volume (δx³) and V_f the fracture volume per cell, which is approximated by counting the number of P_r points per individual cell, then multiplying it with the squared initial point spacing δp and the hydraulic aperture a_h of the fracture. Obviously, the accuracy of V_f crucially depends on the initial point spacing of P_g – the finer the spacing, the better the approximation of V_f. Figure 4c shows that the condition $\mathit{\delta }x/\mathit{\delta }p\ge \mathrm{16}$ delivers sufficiently constant permeability values. In the case that multiple fractures transect the same cell, the permeability tensors are summed, similar to Chen et al. (1999) or Hadgu et al. (2017). However, these cells need additional treatment as they incorporate fracture intersections. We follow the same workflow as presented for individual fractures to map all previously found intersections to the grid cells and formulate their permeability tensors. A horizontal line of the same length as the intersection ($\left|\right|{\mathit{i}}_{\mathrm{1}}-{\mathit{i}}_{\mathrm{0}}\left|\right|$), parallel to the x axis, is represented by a constantly spaced set of points (similar spacing as in the case of a fracture, i.e., δp). The mean point of the line is again located at $(\mathrm{0},\mathrm{0},\mathrm{0})$. We then calculate the rotation matrix R_i (Eq. 10) by using $\overline{\mathit{g}}=[\mathrm{1},\mathrm{0},\mathrm{0}]$ and $\overline{\mathit{n}}=({\mathit{i}}_{\mathrm{1}}-{\mathit{i}}_{\mathrm{0}})/\left|\right|{\mathit{i}}_{\mathrm{1}}-{\mathit{i}}_{\mathrm{0}}\left|\right|$. After identifying the corresponding grid i,j, and k indices as described above, their permeability tensors are increased by using a parameterized intersection permeability (Eq. 9):

$$\begin{array}{}\text{(14)}& {\displaystyle }{\displaystyle }{\mathbf{K}}_{ijk}={\mathbf{K}}_{ijk}+{\displaystyle \frac{V_{\mathrm{i}}}{V_{\mathrm{c}}}}{k}_{\mathrm{i}}\left({\mathbf{R}}_{\mathbf{i}}\left[\begin{array}{ccc}\mathrm{1}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}\end{array}\right]{\mathbf{R}}_{\mathbf{i}}^{\prime }\right).\end{array}$$

V_i represents the intersection volume per cell, which is again approximated by counting the number of P_r points per cell and multiplying it with point spacing δp and the term $\mathit{\pi }{r}_{\mathrm{h}}^{\mathrm{2}}$, whereas r_h denotes the hydraulic radius of the pipe approximating the intersection. Figure 3 shows the resulting ECM structure with 4³ cells of an arbitrary complex DFN, generated with the presented approach. For certain fracture systems (ideally no more than two fractures that fully penetrate the system, e.g., Fig. 4a), the presented approach can be used to derive an analytical solution for permeability by setting δx equal to the system size, resulting in a single permeability tensor for the whole system. Figure 4 demonstrates that incorporating the intersection as a pipe has a significant effect on the shape and absolute value of the permeability tensor at intersections, which could cause an overall permeability increase by almost 1 order of magnitude. However, the exact amount of permeability increase depends on the chosen hydraulic radius of the pipe, and the impact on the overall permeability of the network needs to be evaluated.

5 Fracture-scale intersection flow benchmark

To test the half-hypotenuse assumption (see Fig. 2 for details) as a first-order approximation for the hydraulic radius of the pipe, we conduct a benchmark study in the following. The directional permeabilities of simple fracture crossings with varying orientations calculated from high-resolutions Stokes flow simulations (e.g., Sect. 2) are compared to their analytically derived ECM single-cell counterparts (i.e., δx is equal to the full system size L) using (1) the half-hypotenuse parameterization for intersection flow (fracture-and-pipe model) and (2) omitting this intersection flow parameterization (fracture-only model). For each intersection model, two fully persistent fractures with constant hydraulic apertures of 1.25 mm are placed in a cube of length 10 mm. Two fractures with a dip angle of 90^∘ and dip directions separated by 90^∘ (i.e., 90 and 180^∘) are consecutively rotated counterclockwise by increments of 10^∘ around the center of the cube until a total rotation of 90^∘ is reached. This procedure is repeated nine times, whereas the dip angle of one of the two fractures is consecutively reduced by increments of 10^∘ for each iteration. The dip angle of the remaining fracture is kept constant (i.e., 90) to maintain connectivity in the z direction. This results in a total of 100 different intersection configurations (52 representing direct inlets and outlets of flow, 48 connecting non-boundary flow faces), producing a wide variety of intersection orientations within two opposing octants in the cube (see Fig. 5a for all generated intersection lineaments). For each configuration, we produce a binary voxel model (pore space and matrix) of two crossing parallel-plate fractures (similar to a and d in Fig. 1). Following the approach described in Sect. 2, different pressures at the bottom and top boundary are applied to numerically estimate the directional permeability (setting the remaining boundaries to no-slip yields the vertical permeability component of the permeability tensor, k_z). We systematically increased the numerical resolutions of the Stokes flow simulations (128³, 256³, 512³, and 1024³ voxels) for each intersection configuration (resulting in a total of 400 HPC flow simulations) to determine whether the result at the highest level of detail represents a sufficiently converged solution. This is done by calculating the L2 error norm $\left|\right|{\mathit{\delta }}_k\left|\right|$ according to

$$\begin{array}{}\text{(15)}& {\displaystyle }{\displaystyle }\left|\right|{\mathit{\delta }}_k\left|\right|=\left|{\displaystyle \frac{k_x-k_{\mathrm{1024}}}{k_{\mathrm{1024}}}}\right|,\end{array}$$

where k₁₀₂₄ represents the directional permeability obtained at the highest resolution (i.e., 1024³ voxels) and k_x the directional permeability from simulations with lower resolution (i.e., 128³, 256³, 512³ voxels). The resulting average error norms for all 100 intersection configurations are plotted in Fig. 5b, which demonstrate the convergence towards the numerical result at the highest resolution. With an average error norm of about 0.6 % and a maximum error of 2.4 % for simulations with 512³ voxels compared to the simulations at 1024³ voxels, we assume that the solution at 1024³ voxels represents a sufficiently accurate solution and can furthermore be used to benchmark the tensors generated with the ECM approach. Next, we follow the approach of Sect. 4 to generate a single-cell permeability tensor of each intersection model using a $\mathit{\delta }p/\mathit{\delta }x$ ratio of 16, extract the vertical permeability component of the tensor (k_zz), and compare it with the one resulting from the Stokes flow simulations. The results (Fig. 6) demonstrate that if the intersection connects the two pressure boundary faces (intersection-to-flow-direction angle $\mathit{\gamma }\le \mathrm{40}{}^{\circ }$), the actual permeability obtained from the Stokes simulations is reasonably well reproduced with a small underestimation by the fracture-and-pipe model and heavily underestimated by the fracture-only approach (e.g., Hadgu et al., 2017). Using the half-hypotenuse assumption sufficiently integrates the effect of IFL at the scale of a continuum cell. If intersections that connect no-pressure boundary faces are considered ($\mathit{\gamma }>\mathrm{40}{}^{\circ }$), both models fail to predict the accurate directional permeabilities, indicating that the effect of flow dispersion within the crossing fracture may play a more important role than previously thought. However, the cumulative error box plot in Fig. 6 indicates that both methods give statistically acceptable predictions of the directional permeabilities (median error of 2.7 % for the fracture-and-pipe model and a median error of 7.9 % for the fracture-only model). Thus, the systematic error observed for $\mathit{\gamma }>\mathrm{40}{}^{\circ }$ appears to be negligible.

Figure 7 Panel (a) visualizes the geometrical configurations of all 381 test DFNs as a function of their global P₃₃ value and their ratio between maximum intersection length l_i and system size L. Colors indicate the model size L. Note that high global fracture porosities are predominantly achieved for smaller system sizes for the considered test cases. Panels (b–g) show the underlying DFN structure of the models indicated in (a) with system sizes of 100 m. Fractures are approximated with green disks and intersections with magenta lines. Panels (b–d) have a low prescribed number of fractures ($\mathrm{10}\times L/\mathrm{2}l$), whereas (e–g) have a high prescribed number of fractures ($\mathrm{100}L/\mathrm{2}l$). The ratio of constant fracture length l and system size L for (b, e) is 0.25, 0.5 for (c, f), and 1 for (d, g).

Figure 8 The figure demonstrates the difference between the vertical component of the permeability tensor ΔK_zz resulting from continuum flow simulations of ECMs discretized with the fracture-and-pipe and fracture-only methods described in the text as a function of their respective global fracture porosities (i.e., their P₃₃ value). Colors indicate the constant hydraulic aperture of all fractures in the respective model. Note that high fracture densities are predominantly reached for high fracture apertures.

6 Intersection flow effects at network scales

In the previous section, we demonstrated the effects of IFL on the permeability of systems at the scale of a local ECM cell (i.e., system sizes near the fracture aperture) by comparing analytically derived cell permeabilities to the results of direct flow simulations (Stokes equations). If the intersection orientation aligns with the applied pressure gradient and connects inlet and outlet pressure faces, the permeability of the system is increased (i.e., case b in Fig. 1). To explore the effects of IFL on the permeability of systems at larger scales that cannot be fully resolved with current imaging techniques (i.e., above a few decimeters), we conduct continuum flow simulations (as described in the Appendix) of several test cases in which IFL potentially matters. Following the results of the previous section, this should be the case for fracture formations containing two fracture sets with perpendicular strike and steep dip angles. These so-called cross-joint patterns can be naturally observed (e.g., Gross, 1993; Ruf et al., 1998; Li and Ji, 2021) and are thought to mainly result from local stress field rotations in extensional tectonic settings (Bai et al., 2002; Boersma et al., 2018). Hence, we use the software ADFNE (Alghalandis, 2017) to generate several test DFNs with two orthogonally striking (dip directions are separated by 90^∘) and vertically dipping (dip angle of 90^∘) fracture sets with constant fracture sizes for simplicity. Slight variability in dip angle and direction is introduced by a Fisher dispersion parameter of 20. By this, we ensured that the primary orientation of the formed intersections is oriented parallel to the z direction in the model to provoke the possibly maximal effect of IFL on the network scale. For each test DFN, we vary the following structural parameters during the generation process:

• the cubic overall system side length L by 1, 10, 100, and 1000 m;

• the constant size l of all fractures in the system by 0.25, 0.5, 1, and 2 times the system side length L; and

• the total number of fractures for each set by 10, 50, and 100 times $L/\mathrm{2}l$. The latter rescaling factor is arbitrarily chosen to increase the total number of fractures for systems with lower fracture sizes.

The following parameters are varied in the prescription of the hydraulics of each test DFN:

• the scaling parameter β in a sublinear aperture-length correlation model (e.g., Olson, 2003; Klimczak et al., 2010) by 10⁻⁴, 10⁻³, and 10⁻² – this infers the fractures' hydraulic apertures a_h from their sizes according to a_h=βl^0.5; and

• the isotropic and constant permeability of the rock matrix by 10⁻¹⁷, 10⁻¹⁵, and 10⁻¹³ m².

This results in 432 test DFNs, which are discretized to an ECM with the two different methods already described in the previous section to analytically derive local cell permeability tensors (i.e., the aforementioned fracture-and-pipe and fracture-only method as described in Sect. 4). For this, we start with an ECM grid resolution of $\mathrm{128}\times \mathrm{128}\times \mathrm{128}$ numerical cells to prevent artificial connectivity (e.g., Jackson et al., 2000; Svensson, 2001) for networks with high fracture densities. If this results in nonphysical fracture porosities above unity at the local cell level, we consecutively reduce the grid resolutions by powers of 2 up to $\mathrm{4}\times \mathrm{4}\times \mathrm{4}$ until all local cells have fracture porosities below unity. If this condition cannot be achieved (e.g., for small scales with high fracture densities and high apertures), the model is excluded from the analysis. Furthermore, models with hydraulic apertures above 1 cm were excluded as well, as we assume that the simplification of laminar flow might not hold anymore. This results in a total of 381 test DFNs, whose structure we quantify in a 2D nondimensional parameter system given by (1) the ratio of the maximum intersection length of the system l_i to the system size L and (2) their global P₃₃ value (i.e., global fracture volume divided by total volume according to Dershowitz and Herda, 1992, and hence referred to as the system's global fracture porosity). Figure 7 demonstrates the distribution of the generated test DFNs within this 2D parameter space and shows the DFN structure of some chosen examples. For each geometrical DFN configuration displayed, we compute two (one for each discretization method) effective permeability tensors with the continuum flow procedure described in the Appendix. We quantify the absolute difference of the vertical component in the resulting permeability tensors of the fracture-and-pipe (k_fp) to the fracture-only model (k_fo) by ΔK_zz according to

$$\begin{array}{}\text{(16)}& {\displaystyle }{\displaystyle }\mathrm{\Delta }{K}_{zz}={\displaystyle \frac{|k_{\text{fp}}-k_{\text{fo}}|}{k_{\text{fp}}}},\end{array}$$

which serves as a measure for the magnitude of IFL effects on the network scale. Figure 8 demonstrates that the effect of IFL on network-scale flow depends linearly on global fracture porosity by $\left|\right|\mathrm{\Delta }{k}_{zz}\left|\right|=\mathrm{0.285}{P}_{\mathrm{33}}$. However, the generally very low absolute differences in vertical system permeability indicate a negligible effect of IFL at network scales. Only for networks with global fracture porosities in the range of 30 %–40 % could we observe differences of about 10 %.

Figure 9 Panels (a, b) display the test DFNs with 10 000 and 1000 fractures, respectively. Both are generated with the software ADFNE (Alghalandis, 2017), whereas input parameters are given in the text. Yellow lines depict the location of the slice shown in (c, d). Black lines indicate fractures and magenta spheres the location of fracture intersections.

Figure 10 Panels (a–d) display the norm of the permeability tensor for each cell in an ECM representation of the 10 000 fracture test DFNs displayed in Fig. 9 (a) for grid resolutions of 32³, 64³, 128³, and 256³ voxels, respectively. Panels (e, f) visualize the resulting velocity distribution for an applied pressure gradient in the z direction.

Figure 11 Absolute permeability values k for the six main components of the computed effective permeability tensor (principal components in red, off-diagonal components in black) and the norm of the permeability tensor in magenta as a function of the grid resolution in cubic voxels (number of voxels in the x–y–z direction).

7 Resolution dependency of ECM methods

The resolution dependency of ECM methods to upscale the permeabilities of fracture networks is a crucial aspect that has to be considered to provide accurate upscaling results. Artificial connectivity is one of the main issues that arises if the resolution of the ECM is insufficient, as demonstrated by Jackson et al. (2000) and Svensson (2001). Their results suggest that we can expect accurate upscaling results only when the resolution of the ECM is sufficiently large to resolve the structure of the DFN (i.e., a maximum of two fracture segments and one intersection per cell). As fracture networks typically have a multi-scale character with power-law or lognormal fracture size distributions (e.g., Bonnet et al., 2001; Davy et al., 2006), fulfilling that condition may require very large grid resolutions. Predicting the effective permeability of the DFN by solving the groundwater flow equations (Darcy's law) would then require prior upscaling of the grid cell conductivities (e.g., Zhou et al., 2010; Hauge et al., 2012), depending on the chosen flow solver and the available computational resources. However, averaging or flow-based upscaling approaches may misrepresent network-scale flow characteristics, depending on the chosen coarse grid resolution. Hence, it is often unclear how the resolution dependency affects the accuracy of effective permeability computations and whether flow anisotropy is conserved. In the following, we will demonstrate that using ECMs of DFNs with sufficiently high resolutions can achieve this while avoiding initial upscaling. For this, we compare effective permeability tensors obtained from massively parallelized continuum flow simulations (see Appendix A) for different DFN scenarios with varying resolutions of their equivalent continuum counterparts. We generate two test DFNs utilizing the open-source MATLAB toolbox ADFNE (Alghalandis, 2017). For comparability reasons, we use similar input data as Hadgu et al. (2017), who separated all fractures into three orthogonal sets based on the data reported in SKB (2010). $S\mathrm{1}:\mathrm{90}|\mathrm{090}$, $S\mathrm{2}:\mathrm{90}|\mathrm{000}$, and $S\mathrm{3}:\mathrm{00}|\mathrm{360}$ give the mean dip angle and dip direction for the three fracture sets, respectively, with a constant Fisher distribution concentration value of 5 accounting for variability around the mean. Fracture sizes l are distributed as a power law according to

$$\begin{array}{}\text{(17)}& {\displaystyle }{\displaystyle }l={\left[\left(l_{\mathrm{1}}^{\mathit{\alpha }+\mathrm{1}}-l_{\mathrm{0}}^{\mathit{\alpha }+\mathrm{1}}\right)u+l_{\mathrm{0}}^{\mathit{\alpha }+\mathrm{1}}\right]}^{\mathrm{1}/\mathit{\alpha }+\mathrm{1}},\end{array}$$

where l₁ is the upper cut-off length (500 m) and l₀ the lower cut-off length (15 m), u represents a set of uniformly distributed random numbers in the interval (0,1), and α is the power-law exponent (here $\mathit{\alpha }=-\mathrm{2.5}$). All fracture centers are randomly placed in a cube with 500 m side lengths (the resulting DFNs are displayed in Fig. 9) with a background matrix permeability of 10⁻¹⁸ m². A sublinear scaling of aperture versus length (e.g., Olson, 2003; Klimczak et al., 2010) is employed to correlate the hydraulic apertures a_h of the fractures with their lengths l:

$$\begin{array}{}\text{(18)}& {\displaystyle }{\displaystyle }{a}_{\mathrm{h}}=\mathit{\beta }{l}^{\mathrm{0.5}},\end{array}$$

with a scaling factor β of 10⁻⁴. The only difference between the two test DFNs is the overall fracture number, which is 10 000 for DFN-A (plot a in Fig. 9) and 1000 for DFN-B (b in Fig. 9), such that we obtain a densely and sparsely fractured system, respectively. DFN-A thus represents the scenario of a typical REV network according to Long et al. (1982) and Oda (1985). DFN-B, on the other hand, reflects a flow scenario closer to the percolation threshold with anisotropic, non-REV behavior (Maillot et al., 2016).

After calculating all fracture intersections with ADFNE's built-in function Intersect (see b and d in Fig. 9 for intersection spots in a 2D slice), we use the method presented in Sect. 4, which incorporates the permeability parameterization concepts from Sect. 3, to generate several ECMs with varying grid resolutions. Starting from 4³ voxels and increasing by powers of 2 up to 1024³ voxels yields nine different continuum representations for each test DFN (see Fig. 10 for examples). For each representation, we compute the effective permeability tensor of the DFN by repeatedly solving the Darcy equations in three principal flow directions (see Appendix A for a detailed description). The results are displayed in Fig. 11. For both test DFNs, the norm of the resulting effective permeability tensor ranges within the same order of magnitude. For DFN-A, we obtain a difference of about 30 % from coarse (4³ voxels, $\left|\right|k_{ij}\left|\right|=\mathrm{5.24}\times {\mathrm{10}}^{-\mathrm{11}}$) to fine (1024³ voxels, $\left|\right|k_{ij}\left|\right|=\mathrm{4.03}\times {\mathrm{10}}^{-\mathrm{11}}$) grid resolution, whereas DFN-B shows a larger difference of about 129 % (coarse $\left|\right|k_{ij}\left|\right|=\mathrm{5.07}\times {\mathrm{10}}^{-\mathrm{12}}$, fine $\left|\right|k_{ij}\left|\right|=\mathrm{2.21}\times {\mathrm{10}}^{-\mathrm{12}}$). Thus, the resolution dependence of the absolute permeability is small for fracture networks with an expected REV behavior (DFN-A) and more pronounced if fracture networks with non-REV behavior (DFN-B) are considered. Interestingly, the individual components of the permeability tensor converge to constant values above resolutions of 128³ voxels for both test cases, indicating that anisotropy magnitude depends on the level of detail of the ECM grid.

8 Discussion

Including a pipe flow model into the ECM generation process improves the representation of permeability anisotropy therein and can have impacts on overall permeabilities as well. For example, at the scale of the intersection itself, it significantly modifies the shape and absolute values of the permeability tensor (Fig. 4). However, the errors of the intersection benchmark (2.7 % and 7.9 % for the fracture-and-pipe and fracture-only model, respectively) indicate that, from a statistical perspective, the effect of IFL on overall permeability is usually a second-order effect. This is because the fracture-only ECM discretization approach by default accounts for the increased permeability of intersection cells in an isotropic manner simply by the summation of the individual fracture contributions per cell. The parameter study from Sect. 6 shows that incorporating the effects of IFL into an ECM-based upscaling approach is only necessary for fracture systems that (1) produce a significant amount of fracture intersection (i.e., for fracture systems with two fracture sets of strongly different directions) and (2) have sufficiently high global fracture porosities (above 30 %), regardless of scale. Achieving the latter in our test scenarios is strongly coupled to large fracture apertures (in the order of 1 cm; see Fig. 8) and a ratio between mean intersection and system size close to or above unity (i.e., intersections fully penetrating the system). So, if two dominant orthogonal fractures with large apertures form an intersection that penetrates the whole system along the direction of flow, the effect of IFL influences the system's effective permeability. Due to the nonlinear radius–permeability relation, this may become more important for fractures with apertures above 1 cm. However, for apertures above 1 cm, Reynolds numbers can easily exceed the critical threshold of unity (e.g., Zimmerman and Bodvarsson, 1996), which would require nonlinear concepts to relate fracture and intersection permeability (e.g., Forchheimer's law), as well as Navier–Stokes rather than Stokes simulations.

If small DFNs with sizes closer to the mean hydraulic radius of the intersections (e.g., micro-fracture networks of 1 m size that cannot be naturally resolved with current imaging techniques) are considered for permeability prediction, IFL could play an important role. Then, however, additional factors have to be considered as well. For example, de Dreuzy et al. (2012) have shown that fracture-scale heterogeneity affects network-scale connectivity due to flow channeling caused by closure in the aperture field. This may appear if the ECM cell size is similar to the internal correlation length of the fractures (e.g., Méheust and Schmittbuhl, 2003; Kottwitz et al., 2020), which would ultimately require new concepts to account for deviations from the average flow behavior instead of using fracture permeability parameterizations. A possible solution would be to introduce fracture permeability fluctuations if the ECM cell size is smaller than the individual fracture's correlation length. Unfortunately, the scaling of the correlation length in fractures is poorly understood, so further research is required before integrating these effects. Additionally, the pipe parameterization we use as a first-order approximation for intersection permeability requires refinement to account for irregular shapes, tortuosity, or closure, representing another interesting question to solve in future studies.

For flow simulations at reservoir scales (similar to the test cases considered in Sect. 7), the only computationally feasible solution is to use parameterization concepts (e.g., Sect. 3). For that, we were able to demonstrate that the presented fracture-and-pipe ECM method is capable of providing converged effective permeability tensors if the ECM resolution, i.e., the ratio of system size to discretization step size, is sufficiently large. This resolution dependency for 3D ECMs has not been reported at this level of detail so far but was expected based on previous works of Jackson et al. (2000) and Svensson (2001). The main problem is identified as artificially increased connectivity at lower resolutions, which occurs if the resolution is larger than either the average spacing of the fracture network or the minimal fracture length of the DFN, leading to overestimated permeabilities and misinterpreted anisotropy. Here, we use the average minimal distance of each fracture center to all other fracture centers in the network as a first-order approximation for fracture spacing. With an average spacing of 13.1±4.5 m, continuum grid resolutions above roughly 38 cubic voxels should theoretically start preventing artificial connectivity for DFN-A. For DFN-B with an approximated average spacing of 28.9±10.9 m, the required resolution to damp that effect is even lower (about 17 cubic voxels). Both test DFNs have the same lower cut-off fracture size of 15 m, so artificial connectivity should start decreasing above resolutions of about 33 cubic voxels. Looking at Fig. 11, we observe ongoing permeability convergence at these three mentioned resolutions. We attribute this to the fact that fractures are spaced randomly in space but sampled with a regular grid. Thus, the distance between fracture tips and continuum cell edges might be larger for low resolutions, again causing permeability overestimations. Only above a resolution of 128 cubic voxels do all these effects dampen out, allowing us to declare the solution sufficiently converged with quantitative errors below 10 % for tensor norm and individual components. Hence, we suggest a general upper boundary of a third of the minimal fracture length l₀ as the cell size for an ECM discretization of a DFN to provide constant results.

Based on analytical solutions of flow in fracture networks with constant apertures, Svensson (2001) proposed that the ratio of ECM cell size to hydraulic aperture should not exceed 2 to provide small flow errors. So far, the ratio of cell size to the minimal hydraulic aperture in the system was much larger (about 1260) due to the low scaling factor β of the sublinear aperture to length correlation (Eq. 18). To achieve similar discretization ratios of Svensson (2001) while maintaining a power-law size scaling, we would have to increase β to 10⁻¹, resulting in minimal and maximal apertures of 0.39 and 2.14 m, respectively. As this would violate the assumption of laminar flow conditions within the fractures, we cannot test their hypothesis and rather recommend staying above the maximum hydraulic aperture a_h1 of the system, as otherwise the volume-fraction-based permeability scaling factor in Eqs. (13) and (14) exceeds unity. In that case, parameterization assumptions might not hold anymore, preventing the use of continuum flow methods. However, as demonstrated here, sticking to $l_{\mathrm{0}}/\mathrm{3}>\mathit{\delta }x>a_{\text{h1}}$ as a condition for ECM discretization delivers constant effective permeabilities and conserves flow anisotropy for the upscaling. Within that discretization range, mapping a DFN onto an equivalent continuum grid can be used as a geometric upscaling procedure for further effective permeability analysis. Notably, this range strongly depends on the structural characteristics of the considered DFN, especially on the fracture size distribution and corresponding aperture correlation functions. For some DFNs this may require cropping the fracture size distributions from below to a few multiples of the cell size and compensating for the hydraulic contribution of lower size fractures with a background permeability.

9 Conclusions

This study analyzed the complexity of fracture intersection flow by conducting Stokes flow simulations in simple fracture crossings. Intersections that are aligned with the pressure gradient initiating the flow cause an increase in permeability, as they act similarly to a pipe. This results in intersection flow localization (IFL); i.e., intersections represent preferred pathways for the fluids compared to the connected fractures. We thus extended the state-of-the-art methodology to generate equivalent continuum models (ECMs) for effective permeability computations of discrete fracture networks (DFNs) to incorporate IFL effects. Those are integrated using a directional pipe flow parameterization with a hydraulic radius half the hypotenuse size in a right-angled triangle with side lengths of both intersecting hydraulic apertures. By assessing the permeabilities of fracture intersections numerically, we could demonstrate that for system sizes close to the approximated pipe radius (millimeters to centimeters), the effect of IFL on permeability can be almost 1 order of magnitude. At network scales (meters to kilometers), the impact of IFL on the system's effective permeability is generally minor. Only for fracture systems with high global fracture porosities (above 30 %) do IFL effects become noticeable. Analyzing the effects of IFL on mass transport through fracture networks poses an interesting question for a follow-up study. For example, Makedonska et al. (2016) have shown that early breakthrough times of solute transport through kilometer-scale DFNs are sensitive to local permeability fluctuations. Thus, local permeability increases induced by IFL could potentially affect transport behavior as well.

Next to the effects of IFL, we investigated the resolution dependency of current ECM-based upscaling approaches, as the cell size with which the ECM is discretized represents the most crucial aspect for the accuracy of ECM-based effective permeability predictions. Based on a resolution test with two different DFN scenarios, we suggest that the ECM cell size should be smaller than a third of the minimal fracture size and larger than the maximal hydraulic aperture of the system to conserve constant permeabilities and full anisotropy of flow. Within that range, we conclude that ECM methods equivalently serve as geometric upscaling procedures for fluid flow problems. It is important to note that the accuracy of ECM methods to predict flow is always linked to the quality of the input DFN. Improving the DFN method to better characterize natural fracture systems, especially in terms of fracture termination rules and spatial clustering, is still an ongoing topic of research.

Appendix A: ECM-based effective permeability prediction workflow

In the following, we will explain our method to obtain the effective permeability tensor of continuum cell representations for fractured porous media. The governing equations for steady-state single-phase flow equations for an incompressible, isothermal, and isoviscous fluid without sources and sinks are given in compact form by the following system of mass (Eq. A1) and momentum (Eq. A2) conservation equations:

$$\begin{array}{}\text{(A1)}& {\displaystyle }& {\displaystyle }\mathrm{\nabla }\cdot q=\mathrm{0},\text{(A2)}& {\displaystyle }& {\displaystyle }q=-\mathbf{K}\mathrm{\nabla }P,\end{array}$$

where ∇ and ∇⋅ denote the gradient and divergence operator for global 3D Cartesian coordinates, respectively. The specific discharge (flux) is given by q, pressure by P, and the positive definite and symmetric hydraulic conductivity tensor by K according to

$$\begin{array}{}\text{(A3)}& {\displaystyle }{\displaystyle }\mathbf{K}=\left[\begin{array}{ccc}{k}_{xx}& k_{yx}& k_{zx}\\ k_{yx}& k_{yy}& k_{zy}\\ k_{zx}& k_{zy}& k_{zz}\end{array}\right]{\displaystyle \frac{\mathit{\rho }g}{\mathit{\mu }}},\end{array}$$

with the principal permeability tensor components k_xx, k_yy, and k_zz, the off-diagonal components k_yx, k_zx, and k_zy, fluid density ρ, gravitational acceleration g, and fluid dynamic viscosity μ. We employ a 3D finite-element discretization scheme (e.g., Hughes, 1987; Zienkiewicz and Taylor, 2000; Belytschko et al., 2000; Lin et al., 2014) for Eqs. (A2) and (A1) to simulate boundary-driven pressure diffusion through any input grid consisting of unique permeability tensors. Using the Galerkin method (e.g., Belytschko et al., 2000; Lin et al., 2014), we transform Eq. (A1) into an expression for the nodal residual R according to

$$\begin{array}{}\text{(A4)}& {\displaystyle }{\displaystyle }R=\underset{V}{\int }\mathrm{\nabla }{\mathbf{N}}^{\mathrm{T}}\mathbf{K}\mathrm{\nabla }\mathbf{N}\mathrm{d}VP=\mathrm{0}.\end{array}$$

V denotes the domain volume, N the nodal shape function matrix, and P the nodal pressure. We use eight-node rectangular elements (voxels) with linear interpolation functions (e.g., Zienkiewicz and Taylor, 2000) for volume integral approximation, whereas element integrals are evaluated by the Gauss–Legendre quadrature rule (e.g., Belytschko et al., 2000) over eight integration points with parametric coordinates. Within each element, standard coordinate transformation is employed to compute shape function derivatives with respect to global coordinates ∇N:

$$\begin{array}{}\text{(A5)}& {\displaystyle }{\displaystyle }\mathrm{\nabla }\mathbf{N}={\mathbf{J}}^{-\mathrm{1}}{\mathrm{\nabla }}_L\mathbf{N},\mathbf{J}={\mathrm{\nabla }}_L\mathbf{N}x,\end{array}$$

where ∇_L denotes gradient operator for local 3D element coordinates, J the Jacobian matrix, and x the 3D global element coordinates. After imposing initial pressure conditions at the boundary nodes, the global residual vector R_g is assembled from elemental contributions (e.g., Hughes, 1987) according to Eq. (A4) to solve the linear system of equations,

$$\begin{array}{}\text{(A6)}& {\displaystyle }{\displaystyle }{\mathbf{C}}_{\mathbf{g}}{\mathit{P}}^{\text{new}}={\mathit{R}}_{\mathbf{g}},\end{array}$$

for the unknown pressure P^new. C_g denotes the global coefficient matrix, which is assembled from the nodal coefficient matrix C given by

$$\begin{array}{}\text{(A7)}& {\displaystyle }{\displaystyle }\mathbf{C}=\underset{V}{\int }\mathrm{\nabla }{\mathbf{N}}^{\mathrm{T}}\mathbf{K}\mathrm{\nabla }\mathbf{N}\mathrm{d}V.\end{array}$$

Following this, we evaluate the Darcy velocities at the integration points u based on the newly solved nodal pressures by

$$\begin{array}{}\text{(A8)}& {\displaystyle }{\displaystyle }u=\mathbf{K}\mathrm{\nabla }\mathbf{N}{\mathit{P}}^{\text{new}},\end{array}$$

where the velocity vectors on the nodes are averaged from the neighboring integration points.

Figure A1 Pressure boundary conditions for an applied gradient in the z direction. Here, top and bottom faces experience constant pressures of 1 and 0 Pa, respectively. A linearly interpolated pressure distribution is applied at the remaining four boundary faces, as indicated by the colored wedges next to the side faces of the model. Thus, the principal direction of flow is in the z direction, allowing for the calculation of the terms related to the z component of the permeability tensor according to Eq. (A9).

Figure A2 Benchmark case 1 from Berre et al. (2020). Panel (a) shows the benchmark geometry of an embedded fracture (aperture of 10⁻² m) in a matrix with a hydraulic conductivity of 10⁻⁶. The hydraulic conductivity in the gray band at the bottom is increased to 10⁻⁶. Constant pressures of 4 and 1 Pa are applied at the inlet band (blue) and outlet band (red), respectively. The diagonal light gray line through the model indicates the sampling line for the pressures shown in (b). The pressure distribution is plotted as a function of arc length of the gray line in (a), and the results of different resolutions are compared to the benchmark target field obtained from 17 different numerical methods. The dark gray region illustrates the area between the 10th and 90th percentiles for the highest refinement level of the benchmarked methods, whereas the light gray region illustrates the same area for their lowest refinement level.

Three principal directions of the applied pressure gradient have to be considered to predict the full tensor of permeability. Thus, the flow simulation procedure has to be repeated three times such that each principal flow direction (x, y, and z direction in a Cartesian coordinate system) is covered. For each iteration, two constant pressure values are applied at two opposing boundary faces (e.g., lower and upper face in a cube for principal flow in the z direction), and the same linear interpolation between those two values is applied at the remaining four boundary faces (see Fig. A1 for an example). This ensures the capture of both the diagonal and off-diagonal terms of the permeability tensor properly, which are computed by substituting the volume average $\overline{u}$ of all nodal velocity vectors u_I (see Eq. 3) into Darcy's law for flow through porous media in the form of Eq. (4). Figure A1 displays the situation of a vertically aligned pressure gradient ($\mathrm{\Delta }{P}_z=\frac{\mathit{\delta }P}{\mathit{\delta }z}$). The corresponding entries in the permeability tensor are computed according to

$$\begin{array}{}\text{(A9)}& {\displaystyle }{\displaystyle }\left[\begin{array}{c}{k}_{zx}\\ k_{zy}\\ k_{zz}\end{array}\right]={\displaystyle \frac{\mathit{\mu }}{\mathrm{\Delta }{P}_z}}\left[\begin{array}{c}{\overline{u}}_x\\ {\overline{u}}_y\\ {\overline{u}}_z\end{array}\right],\end{array}$$

and vice versa for the iterations with pressure gradients in the x and y direction to obtain the permeability tensor as shown in Eq. (A3).

The single-continuum discretization scheme used might appear simplistic compared to more sophisticated mesh representations (see Berre et al., 2019). However, the merits of our approach rather lay (1) in a fully anisotropic permeability representation of the individual continuum cells and (2) massive parallelization and HPC optimization. Utilizing the parallelization framework of PETSc (Balay et al., 2018) and their multigrid preconditioned solvers significantly reduces the computational cost, allowing simulations to be run routinely with 10⁹ individual grid cells. An increase in grid resolution compensates for the benefits of using conforming meshes or multi-continuum formulations (e.g., Berre et al., 2019). To test this, we compare our modeling procedure against benchmark case 1 from Berre et al. (2020), who compare 17 different methods of simulating single-phase flow in fractured porous media. The initial setup (displayed in a in Fig. A2) consists of an inclined fracture with a hydraulic aperture of 10⁻² m embedded in a cube of 100 m length with a matrix hydraulic conductivity of 10⁻⁶ m², whereas the hydraulic conductivity of a small band of 10 m width at the bottom is increased to 10⁻⁵ m². We prescribe these two values as background permeabilities and use the methodology described in Sect. 3 to incorporate fracture permeability accordingly. The boundary conditions are given by small pressure inlet (4 Pa) and outlet (1 Pa) bands as indicated in Fig. A2a. The comparison of the pressure distribution (Fig. A2b) highlights the fact that with a resolution of 32³ voxels, we can already obtain a good fit with the benchmark target field. This thus suggests that our modeling procedure is sufficiently correct for effective permeability predictions.