Research article 05 Oct 2021
Research article  05 Oct 2021
Investigating the effects of intersection flow localization in equivalentcontinuumbased upscaling of flow in discrete fracture networks
 ^{1}Johannes Gutenberg University, Institute of Geosciences, JohannJoachimBecherWeg 21, 55128 Mainz, Germany
 ^{2}Igem, Institute for Geothermal Resource Management, Berlinstr. 107a, 55411 Bingen, Germany
 ^{3}Johannes Gutenberg University, M3ODEL – Mainz Institute of Multiscale Modeling, Staudingerweg 7, 55128 Mainz, Germany
 ^{1}Johannes Gutenberg University, Institute of Geosciences, JohannJoachimBecherWeg 21, 55128 Mainz, Germany
 ^{2}Igem, Institute for Geothermal Resource Management, Berlinstr. 107a, 55411 Bingen, Germany
 ^{3}Johannes Gutenberg University, M3ODEL – Mainz Institute of Multiscale Modeling, Staudingerweg 7, 55128 Mainz, Germany
Correspondence: Maxmilian O. Kottwitz (mkottwi@unimainz.de)
Hide author detailsCorrespondence: Maxmilian O. Kottwitz (mkottwi@unimainz.de)
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 singlephase 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 fractureandpipe model. Our approach is compared against an ECM method that does not take IFL into account by performing ECMbased 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 networkscale 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 ECMbased 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 resolutionindependent. Within that range, ECMs are applicable to upscale flow in fracture networks.
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_{2} 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 longstanding 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 socalled direct flow modeling approaches crucially rely on a digital representation of a rock that separates pore space from the matrix, which results from highresolution Xray computed tomographies (Andrä et al., 2013a; Cnudde and Boone, 2013). However, they are limited in maximum scannable size and respective tradeoff 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), socalled 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 socalled 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 multiscale 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 reducedorder 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 parallelplate assumption (Brown, 1995; Oron and Berkowitz, 1998; Méheust and Schmittbuhl, 2000). Semiempirical functions derived from numerical simulations in roughwalled 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) cubiclaw assumptions. Improved discretization techniques with individual fracture treatment (DFN method), inclusion of matrix properties in multidimensional meshes (discrete fracture and matrix – DFM – method), and multicontinuum methods come at the cost of high computational expenses. Discretizing the fractured media as equivalent singlecontinuum 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 multiscale 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 multiscale 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 socalled 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, steadystate, singlephase 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. Socalled 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 highresolution simulations. However, assessing permeabilities in a MonteCarlolike 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 stateoftheart 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 highperformance computing (HPC) optimized finiteelement 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.
Fluid flow in porous and fractured media is described by the wellknown 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, steadystate flow at constant temperature is defined by the Stokes momentum balance (Eq. 1) and continuity (Eq. 2) equations (Bear, 1972):
with the fluid's dynamic viscosity μ, pressure P, and velocity vector $\mathit{v}=({v}_{x},{v}_{y},{v}_{z})$. ∇, ∇⋅, and ∇^{2} denote the gradient, divergence, and Laplace operator for 3D Cartesian coordinates, respectively.
Here, the 3D staggeredgrid, finitedifference 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 noslip (velocity component normal to the boundary is zero), enables the prediction of the medium's directional permeability. After obtaining the steadystate solution, the volume integral of the pressuregradientaligned velocity component v_{z} (e.g., Osorno et al., 2015) is computed according to
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
with intrinsic permeability k and crosssectional area A in combination with the fact that $Q=\overline{v}A$, the directional permeability k_{z} is calculated by
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^{3}, 256^{3}, 512^{3}, and 1024^{3} 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 nopressure 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.
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 cubiclaw 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
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, crosssectional area A=wa_{m}) leaves the intrinsic permeability of a fracture k_{f} defined by
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 readytouse parameterization concept tailored for fracture intersections. The simulations shown in Fig. 1 suggest that the flow in the intersection is approximately pipelike. 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
Again, combining this equation with Darcy's law (Eq. 4, crosssectional area A=πr^{2}) results in the following expression for the intrinsic permeability of a pipe k_{i}:
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 firstorder approximation, we use half the size of the hypotenuse in a rightangled triangle whose legs are given by the two intersecting apertures (called halfhypotenuse assumption in the following; see Fig. 2 for details). This delivers sufficiently good results, as will be demonstrated later (Fig. 6).
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 DFNmapped 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 fractureandpipe 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_{0}, 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_{0} and i_{1}, 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
with ×, ⋅, and $\left\right\mathit{x}\left\right$ denoting the crossproduct, dot product, and vector norm of x, respectively. I represents the 3 × 3 identity matrix and C the crossproduct matrix of the rotation axis $\mathit{w}=[{w}_{x},{w}_{y},{w}_{z}]$:
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_{0} to produce a rotated set of points P_{r} representing the fracture in its actual 3D position:
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
V_{c} denotes the cell volume (δx^{3}) 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):
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^{3} 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.
To test the halfhypotenuse assumption (see Fig. 2 for details) as a firstorder 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 highresolutions Stokes flow simulations (e.g., Sect. 2) are compared to their analytically derived ECM singlecell counterparts (i.e., δx is equal to the full system size L) using (1) the halfhypotenuse parameterization for intersection flow (fractureandpipe model) and (2) omitting this intersection flow parameterization (fractureonly 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 nonboundary 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 parallelplate 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 noslip yields the vertical permeability component of the permeability tensor, k_{z}). We systematically increased the numerical resolutions of the Stokes flow simulations (128^{3}, 256^{3}, 512^{3}, and 1024^{3} 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
where k_{1024} represents the directional permeability obtained at the highest resolution (i.e., 1024^{3} voxels) and k_{x} the directional permeability from simulations with lower resolution (i.e., 128^{3}, 256^{3}, 512^{3} 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^{3} voxels compared to the simulations at 1024^{3} voxels, we assume that the solution at 1024^{3} 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 singlecell 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 (intersectiontoflowdirection 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 fractureandpipe model and heavily underestimated by the fractureonly approach (e.g., Hadgu et al., 2017). Using the halfhypotenuse assumption sufficiently integrates the effect of IFL at the scale of a continuum cell. If intersections that connect nopressure 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 fractureandpipe model and a median error of 7.9 % for the fractureonly model). Thus, the systematic error observed for $\mathit{\gamma}>\mathrm{40}{}^{\circ}$ appears to be negligible.
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 socalled crossjoint 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 aperturelength correlation model (e.g., Olson, 2003; Klimczak et al., 2010) by 10^{−4}, 10^{−3}, and 10^{−2} – 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^{−17}, 10^{−15}, and 10^{−13} m^{2}.
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 fractureandpipe and fractureonly 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_{33} 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 fractureandpipe (k_{fp}) to the fractureonly model (k_{fo}) by ΔK_{zz} according to
which serves as a measure for the magnitude of IFL effects on the network scale. Figure 8 demonstrates that the effect of IFL on networkscale 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 %.
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 multiscale character with powerlaw 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 flowbased upscaling approaches may misrepresent networkscale 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 opensource 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
where l_{1} is the upper cutoff length (500 m) and l_{0} the lower cutoff length (15 m), u represents a set of uniformly distributed random numbers in the interval (0,1), and α is the powerlaw 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^{−18} m^{2}. 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:
with a scaling factor β of 10^{−4}. The only difference between the two test DFNs is the overall fracture number, which is 10 000 for DFNA (plot a in Fig. 9) and 1000 for DFNB (b in Fig. 9), such that we obtain a densely and sparsely fractured system, respectively. DFNA thus represents the scenario of a typical REV network according to Long et al. (1982) and Oda (1985). DFNB, on the other hand, reflects a flow scenario closer to the percolation threshold with anisotropic, nonREV behavior (Maillot et al., 2016).
After calculating all fracture intersections with ADFNE's builtin 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^{3} voxels and increasing by powers of 2 up to 1024^{3} 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 DFNA, we obtain a difference of about 30 % from coarse (4^{3} voxels, $\left\right{k}_{ij}\left\right=\mathrm{5.24}\times {\mathrm{10}}^{\mathrm{11}}$) to fine (1024^{3} voxels, $\left\right{k}_{ij}\left\right=\mathrm{4.03}\times {\mathrm{10}}^{\mathrm{11}}$) grid resolution, whereas DFNB 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 (DFNA) and more pronounced if fracture networks with nonREV behavior (DFNB) are considered. Interestingly, the individual components of the permeability tensor converge to constant values above resolutions of 128^{3} voxels for both test cases, indicating that anisotropy magnitude depends on the level of detail of the ECM grid.
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 fractureandpipe and fractureonly model, respectively) indicate that, from a statistical perspective, the effect of IFL on overall permeability is usually a secondorder effect. This is because the fractureonly 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 ECMbased 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., microfracture 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 fracturescale heterogeneity affects networkscale 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 firstorder 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 fractureandpipe 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 firstorder 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 DFNA. For DFNB 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 cutoff 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_{0} 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 powerlaw size scaling, we would have to increase β to 10^{−1}, 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 volumefractionbased 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.
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 stateoftheart 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 rightangled 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 followup study. For example, Makedonska et al. (2016) have shown that early breakthrough times of solute transport through kilometerscale 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 ECMbased upscaling approaches, as the cell size with which the ECM is discretized represents the most crucial aspect for the accuracy of ECMbased 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.
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 steadystate singlephase 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:
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
with the principal permeability tensor components k_{xx}, k_{yy}, and k_{zz}, the offdiagonal components k_{yx}, k_{zx}, and k_{zy}, fluid density ρ, gravitational acceleration g, and fluid dynamic viscosity μ. We employ a 3D finiteelement 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 boundarydriven 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
V denotes the domain volume, N the nodal shape function matrix, and P the nodal pressure. We use eightnode 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:
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,
for the unknown pressure P^{new}. C_{g} denotes the global coefficient matrix, which is assembled from the nodal coefficient matrix C given by
Following this, we evaluate the Darcy velocities at the integration points u based on the newly solved nodal pressures by
where the velocity vectors on the nodes are averaged from the neighboring integration points.
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 offdiagonal 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
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 singlecontinuum 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^{9} individual grid cells. An increase in grid resolution compensates for the benefits of using conforming meshes or multicontinuum 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 singlephase 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^{−2} m embedded in a cube of 100 m length with a matrix hydraulic conductivity of 10^{−6} m^{2}, whereas the hydraulic conductivity of a small band of 10 m width at the bottom is increased to 10^{−5} m^{2}. 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^{3} 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.
All numerical simulation tools used in this paper are past or modified versions of the opensource software package LaMEM (https://bitbucket.org/bkaus/lamem/src, last access: 11 September 2021, Kaus et al., 2016). The version for finitedifference Stokesflow simulations is available at: https://doi.org/10.5281/zenodo.5520927 (last access: 22 September 2021, Kottwitz et al., 2021a). The version for finiteelement Darcyflow simulations is available at: https://doi.org/10.5281/zenodo.5521001 (last access: 22 September 2021, Kottwitz et al., 2021b). The scripts to reproduce the input data sets are available upon request from the authors.
MOK wrote the initial draft of the paper, performed numerical simulations, analyzed the data, and generated the figures. AAP supervised and helped design the study, helped implement the numerical methods, and edited the paper. SA helped design the study and edited the paper. BJPK edited the paper and discussed the results.
The authors declare that they have no conflict of interest.
Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
The authors gratefully acknowledge the computing time granted on the supercomputer Mogon II at Johannes Gutenberg University Mainz (hpc.unimainz.de). The authors sincerely thank two anonymous referees for reviewing this paper.
This work has been funded by the Federal Ministry of Education and Research (BMBF; program GEO:N, grant no. 03G0865A) and the M3ODEL consortium at Johannes Gutenberg University Mainz.
This paper was edited by Juliane Dannberg and reviewed by two anonymous referees.
Alghalandis, Y. F.: ADFNE: Open source software for discrete fracture network engineering, two and three dimensional applications, Comput. Geosci., 102, 1–11, https://doi.org/10.1016/j.cageo.2017.02.002, 2017. a, b, c, d
Andrä, H., Combaret, N., Dvorkin, J., Glatt, E., Han, J., Kabel, M., Keehm, Y., Krzikalla, F., Lee, M., Madonna, C., Marsh, M., Mukerji, T., Saenger, E. H., Sain, R., Saxena, N., Ricker, S., Wiegmann, A., and Zhan, X.: Digital rock physics benchmarksPart I: Imaging and segmentation, Comput. Geosci., 50, 25–32, https://doi.org/10.1016/j.cageo.2012.09.005, 2013a. a
Andrä, H., Combaret, N., Dvorkin, J., Glatt, E., Han, J., Kabel, M., Keehm, Y., Krzikalla, F., Lee, M., Madonna, C., Marsh, M., Mukerji, T., Saenger, E. H., Sain, R., Saxena, N., Ricker, S., Wiegmann, A., and Zhan, X.: Digital rock physics benchmarkspart II: Computing effective properties, Comput. Geosci., 50, 33–43, https://doi.org/10.1016/j.cageo.2012.09.008, 2013b. a
Bai, T., Maerten, L., Gross, M. R., and Aydin, A.: Orthogonal cross joints: Do they imply a regional stress rotation?, J. Struct. Geol., 24, 77–88, https://doi.org/10.1016/S01918141(01)000505, 2002. a
Balay, S., Abhyankar, S., Adams, M., Brown, J., Brune, P., Buschelman, K., Dalcin, L., Dener, A., Eijkhout, V., Gropp, W., Karpeyev, D., Kaushik, D., Knepley, M., May, D., Curfman McInnes, L., Mills, R., Munson, T., Rupp, K., Sanan, P., Smith, B., Zampini, S., Zhang, H., and Zhang, H.: PETSc Users Manual: Revision 3.10, Tech. rep., Argonne National Lab.(ANL), Argonne, IL (US), 2018. a, b
Batchelor, G. K.: An Introduction to Fluid Dynamics, Cambridge Mathematical Library, Cambridge University Press, Cambridge, United Kingdom, https://doi.org/10.1017/CBO9780511800955, 1967. a
Bear, J.: Dynamics of fluids in porous media, Elsevier, New York, 1972. a, b, c, d, e
Belytschko, T., Liu, W. K., Moran, B., and Elkhodary, K.: Nonlinear finite elements for continua and structures, John Wiley and Sons, Chichester, West Sussex, United Kingdom, 2000. a, b, c
Berre, I., Doster, F., and Keilegavlen, E.: Flow in fractured porous media: A review of conceptual models and discretization approaches, Transport Porous Med., 130, 215–236, https://doi.org/10.1007/s1124201811716, 2019. a, b, c
Berre, I., Boon, W. M., Flemisch, B., Fumagalli, A., Gläser, D., Keilegavlen, E., Scotti, A., Stefansson, I., Tatomir, A., Brenner, K., Burbulla, S., Devloo, P., Duran, O., Favino, M., Hennicker, J., Lee, I.H., Lipnikov, K., Masson, R., Mosthaf, K., Nestola, M. G. C., Ni, C.F., Nikitin, K., Schädle, P., Svyatskiy, D., Yanbarisov, R., and Zulian, P.: Verification benchmarks for singlephase flow in threedimensional fractured porous media, Adv. Water Res., 147, 103759, https://doi.org/10.1016/j.advwatres.2020.103759, 2020. a, b
Boersma, Q., Hardebol, N., Barnhoorn, A., and Bertotti, G.: Mechanical Factors Controlling the Development of Orthogonal and Nested Fracture Network Geometries, Rock Mech. Rock Eng., 51, 3455–3469, https://doi.org/10.1007/s0060301815528, 2018. a
Bogdanov, I. I., Mourzenko, V. V., Thovert, J.F., and Adler, P. M.: Effective permeability of fractured porous media in steady state flow, Water Resour. Res., 39, 1023, https://doi.org/10.1029/2001WR000756, 2003. a
Bonnet, E., Bour, O., Odling, N. E., Davy, P., Main, I., Cowie, P., and Berkowitz, B.: Scaling of fracture systems in geological media, Rev. Geophys., 39, 347–383, https://doi.org/10.1029/1999RG000074, 2001. a, b
Botros, F. E., Hassan, A. E., Reeves, D. M., and Pohll, G.: On mapping fracture networks onto continuum, Water Resour. Res., 44, W08435, https://doi.org/10.1029/2007WR006092, 2008. a, b, c
Brown, S. R.: Fluid flow through rock joints: the effect of surface roughness, J. Geophys. Res.Sol. Ea., 92, 1337–1347, https://doi.org/10.1029/JB092iB02p01337, 1987. a, b
Brown, S. R.: Simple mathematical model of a rough fracture, J. Geophys. Res.Sol. Ea., 100, 5941–5952, https://doi.org/10.1029/94JB03262, 1995. a
Cacas, M. C., Ledoux, E., de Marsily, G., Tillie, B., Barbreau, A., Durand, E., Feuga, B., and Peaudecerf, P.: Modeling fracture flow with a stochastic discrete fracture network: calibration and validation: 1. The flow model, Water Resour. Res., 26, 479–489, https://doi.org/10.1029/WR026i003p00479, 1990. a, b
Cartwright, J. and Huuse, M.: 3D seismic technology: the geological “Hubble', Basin Res., 17, 1–20, https://doi.org/10.1111/j.13652117.2005.00252.x, 2005. a
Chen, M., Bai, M., and Roegiers, J.C.: Permeability tensors of anisotropic fracture networks, Math. Geol., 31, 335–373, https://doi.org/10.1023/A:1007534523363, 1999. a, b, c
Chen, T., Clauser, C., Marquart, G., Willbrand, K., and Mottaghy, D.: A new upscaling method for fractured porous media, Adv. Water Resour., 80, 60–68, https://doi.org/10.1016/j.advwatres.2015.03.009, 2015. a
Cnudde, V. and Boone, M. N.: Highresolution Xray computed tomography in geosciences: A review of the current technology and applications, EarthSci. Rev., 1–17, https://doi.org/10.1016/j.earscirev.2013.04.003, 2013. a
Darcel, C., Bour, O., Davy, P., and de Dreuzy, J. R.: Connectivity properties of twodimensional fracture networks with stochastic fractal correlation, Water Resour. Res., 39, 1272, https://doi.org/10.1029/2002WR001628, 2003. a
Darcy, H. P. G.: Les Fontaines publiques de la ville de Dijon. Exposition et application des principes à suivre et des formules à employer dans les questions de distribution d'eau, etc, V. Dalamont, Paris, 1856. a, b, c
Davy, P., Bour, O., De Dreuzy, J.R., and Darcel, C.: Flow in multiscale fractal fracture networks, Geol. Soc. Sp., 261, 31–45, https://doi.org/10.1144/GSL.SP.2006.261.01.03, 2006. a, b
Davy, P., Le Goc, R., and Darcel, C.: A model of fracture nucleation, growth and arrest, and consequences for fracture density and scaling, J. Geophys. Res.Sol. Ea., 118, 1393–1407, https://doi.org/10.1002/jgrb.50120, 2013. a
de Dreuzy, J.R., Méheust, Y., and Pichot, G.: Influence of fracture scale heterogeneity on the flow properties of threedimensional discrete fracture networks (DFN), J. Geophys. Res.Sol. Ea., 117, B11207, https://doi.org/10.1029/2012JB009461, 2012. a, b, c
Dershowitz, W. S. and Einstein, H. H.: Characterizing rock joint geometry with joint system models, Rock Mech. Rock Eng., 21, 21–51, https://doi.org/10.1007/BF01019674, 1988. a, b
Dershowitz, W. S. and Herda, H. H.: Interpretation of fracture spacing and intensity, 757–766, 1992. a, b
Eichheimer, P., Thielmann, M., Popov, A., Golabek, G. J., Fujita, W., Kottwitz, M. O., and Kaus, B. J. P.: Porescale permeability prediction for Newtonian and nonNewtonian fluids, Solid Earth, 10, 1717–1731, https://doi.org/10.5194/se1017172019, 2019. a, b
Eichheimer, P., Thielmann, M., Fujita, W., Golabek, G. J., Nakamura, M., Okumura, S., Nakatani, T., and Kottwitz, M. O.: Combined numerical and experimental study of microstructure and permeability in porous granular media, Solid Earth, 11, 1079–1095, https://doi.org/10.5194/se1110792020, 2020. a, b
Gross, M. R.: The origin and spacing of cross joints: examples from the Monterey Formation, Santa Barbara Coastline, California, J. Struct. Geol., 15, 737–751, https://doi.org/10.1016/01918141(93)90059J, 1993. a
Hadgu, T., Karra, S., Kalinina, E., Makedonska, N., Hyman, J. D., Klise, K., Viswanathan, H. S., and Wang, Y.: A comparative study of discrete fracture network and equivalent continuum models for simulating flow and transport in the far field of a hypothetical nuclear waste repository in crystalline host rock, J. Hydrol., 553, 59–70, https://doi.org/10.1016/j.jhydrol.2017.07.046, 2017. a, b, c, d, e, f, g
Hauge, V. L., Lie, K.A., and Natvig, J. R.: Flowbased coarsening for multiscale simulation of transport in porous media, Computat. Geosci., 16, 391–408, https://doi.org/10.1007/s105960119230x, 2012. a, b
Healy, D., Rizzo, R. E., Cornwell, D. G., Farrell, N. J. C., Watkins, H., Timms, N. E., GomezRivas, E., and Smith, M.: FracPaQ: A MATLAB™ toolbox for the quantification of fracture patterns, J. Struct. Geol., 95, 1–16, https://doi.org/10.1016/j.jsg.2016.12.003, 2017. a
Hughes, T. J. R.: The Finite Element Method: Linear Static and Dynamic Finite Element Analysis, PrenticeHall Inc., Englewood Cliffs, New Jersey, 1987. a, b
Hyman, J. D., Karra, S., Makedonska, N., Gable, C. W., Painter, S. L., and Viswanathan, H. S.: dfnWorks: A discrete fracture network framework for modeling subsurface flow and transport, Comput. Geosci., 84, 10–19, https://doi.org/10.1016/j.cageo.2015.08.001, 2015. a
Jackson, C. P., Hoch, A. R., and Todman, S.: Selfconsistency of a heterogeneous continuum porous medium representation of a fractured medium, Water Resour. Res., 36, 189–202, https://doi.org/10.1029/1999WR900249, 2000. a, b, c, d, e
Jing, L.: A review of techniques, advances and outstanding issues in numerical modelling for rock mechanics and rock engineering, Int. J. Rock Mech. Min., 40, 283–353, https://doi.org/10.1016/S13651609(03)000133, 2003. a
Kaus, B. J., Popov, A. A., Baumann, T., Pusok, A., Bauville, A., Fernandez, N., and Collignon, M.: Forward and inverse modelling of lithospheric deformation on geological timescales, in: Proceedings of NIC Symposium, available at: http://hdl.handle.net/2128/9842 (last access: 11 September 2021), 2016. a, b
Klimczak, C., Schultz, R. A., Parashar, R., and Reeves, D. M.: Cubic law with aperturelength correlation: implications for network scale fluid flow, Hydrogeol. J., 18, 851–862, https://doi.org/10.1007/s1004000905726, 2010. a, b
Kottwitz, M. O., Popov, A. A., Baumann, T. S., and Kaus, B. J. P.: The hydraulic efficiency of single fractures: correcting the cubic law parameterization for selfaffine surface roughness and fracture closure, Solid Earth, 11, 947–957, https://doi.org/10.5194/se119472020, 2020. a, b, c, d, e
Kottwitz, M. O., Popov, A. A., Abe, S., Kaus, B. J. P.: LaMEM version used for finitedifference Stokesflow computations in “Investigating the effects of intersection flow localization in equivalent continuumbased upscaling of flow in discrete fracture networks”, Zenodo [code], https://doi.org/10.5281/zenodo.5520927, last access: 22 September 2021. a
Kottwitz, M. O., Popov, A. A., Abe, S., Kaus, B. J. P.: LaMEM version used for finiteelement Darcyflow computations in “Investigating the effects of intersection flow localization in equivalent continuumbased upscaling of flow in discrete fracture networks”, Zenodo [code], https://doi.org/10.5281/zenodo.5521001, last access: 22 September 2021. a
La Pointe, P. R., Wallmann, P., and Follin, S.: Estimation of effective block conductivities based on discrete network analyses using data from the Äspö site (SKBTR–9515), Tech. rep., Swedish Nuclear Fuel and Waste Management Co., 1995. a
Lei, Q., Latham, J.P., and Tsang, C.F.: The use of discrete fracture networks for modelling coupled geomechanical and hydrological behaviour of fractured rocks, Comput. Geotech., 85, 151–176, https://doi.org/10.1016/j.compgeo.2016.12.024, 2017. a
Leung, C. T. O., Hoch, A. R., and Zimmerman, R. W.: Comparison of discrete fracture network and equivalent continuum simulations of fluid flow through twodimensional fracture networks for the DECOVALEX–2011 project, Mineral. Mag., 76, 3179–3190, https://doi.org/10.1180/minmag.2012.076.8.31, 2012. a, b
Li, B., Mo, Y., Zou, L., Liu, R., and Cvetkovic, V.: Influence of surface roughness on fluid flow and solute transport through 3D crossed rock fractures, J. Hydrol., 582, 124284, https://doi.org/10.1016/j.jhydrol.2019.124284, 2020. a
Li, L. and Ji, S.: A new interpretation for formation of orthogonal joints in quartz sandstone, Journal of Rock Mechanics and Geotechnical Engineering, 13, 289–299, https://doi.org/10.1016/j.jrmge.2020.08.003, 2021. a
Lichtner, P. C., Hammond, G. E., Lu, C., Karra, S., Bisht, G., Andre, B., Mills, R., and Kumar, J.: PFLOTRAN User Manual: A Massively Parallel Reactive Flow and Transport Model for Describing Surface and Subsurface Processes, https://doi.org/10.2172/1168703, 2015. a
Lie, K.A.: Upscaling Petrophysical Properties, in: An Introduction to Reservoir Simulation Using MATLAB/GNU Octave: User Guide for the MATLAB Reservoir Simulation Toolbox (MRST), edited by: Lie, K.A., pp. 558–596, Cambridge University Press, Cambridge, https://doi.org/10.1017/9781108591416.020, 2019. a
Lin, G., Liu, J., Mu, L., and Ye, X.: Weak Galerkin finite element methods for Darcy flow: Anisotropy and heterogeneity, J. Comput. Phys., 276, 422–437, https://doi.org/10.1016/j.jcp.2014.07.001, 2014. a, b
Long, J. C. S., Remer, J. S., Wilson, C. R., and Witherspoon, P. A.: Porous media equivalents for networks of discontinuous fractures, Water Resour. Res., 18, 645–658, https://doi.org/10.1029/WR018i003p00645, 1982. a, b, c, d, e
Maillot, J., Davy, P., Le Goc, R., Darcel, C., and de Dreuzy, J. R.: Connectivity, permeability, and channeling in randomly distributed and kinematically defined discrete fracture network models, Water Resour. Res., 52, 8526–8545, https://doi.org/10.1002/2016WR018973, 2016. a, b
Makedonska, N., Hyman, J. D., Karra, S., Painter, S. L., Gable, C. W., and Viswanathan, H. S.: Evaluating the effect of internal aperture variability on transport in kilometer scale discrete fracture networks, Adv. Water Resour., 94, 486–497, https://doi.org/10.1016/j.advwatres.2016.06.010, 2016. a
Malehmir, A., Bellefleur, G., Koivista, E., and Juhlin, C.: Pros and cons of 2D vs 3D seismic mineral exploration surveys, First Break, 35, 8, https://doi.org/10.3997/13652397.35.8.89805, 2017. a
McDonald, M. G. and Harbaugh, A. W.: A modular threedimensional finitedifference groundwater flow model, US Geological Survey, United States, 1988. a
Méheust, Y. and Schmittbuhl, J.: Flow enhancement of a rough fracture, Geophys. Res. Lett., 27, 2989–2992, https://doi.org/10.1029/1999GL008464, 2000. a
Méheust, Y. and Schmittbuhl, J.: Scale effects related to flow in rough fractures, Pure Appl. Geophys., 160, 1023–1050, https://doi.org/10.1007/PL00012559, 2003. a, b
Mourzenko, V. V., Thovert, J.F., and Adler, P. M.: Conductivity and Transmissivity of a Single Fracture, Transport Porous Med., 123, 235–256, https://doi.org/10.1007/s112420181037y, 2018. a
Neuman, S. P.: Trends, prospects and challenges in quantifying flow and transport through fractured rocks, Hydrogeol. J., 13, 124–147, https://doi.org/10.1007/s1004000403972, 2005. a
Oda, M.: Permeability tensor for discontinuous rock masses, Géotechnique, 35, 483–495, https://doi.org/10.1680/geot.1985.35.4.483, 1985. a, b, c
Odling, N. E., Gillespie, P., Bourgine, B., Castaing, C., Chiles, J. P., Christensen, N. P., Fillion, E., Genter, A., Olsen, C., Thrane, L., Trice, R., Aarseth, E., Walsh, J. J., and Watterson, J.: Variations in fracture system geometry and their implications for fluid flow in fractures hydrocarbon reservoirs, Petrol. Geosci., 5, 373–384, https://doi.org/10.1144/petgeo.5.4.373, 1999. a, b
Olson, J. E.: Sublinear scaling of fracture aperture versus length: An exception or the rule?, J. Geophys. Res.Sol. Ea., 108, 2413, https://doi.org/10.1029/2001JB000419, 2003. a, b
Oron, A. P. and Berkowitz, B.: Flow in rock fractures: The local cubic law assumption reexamined, Water Resour. Res., 34, 2811–2825, https://doi.org/10.1029/98WR02285, 1998. a
Ortega, O. J., Marrett, R. A., and Laubach, S. E.: A scaleindependent approach to fracture intensity and average spacing measurement, AAPG Bull., 90, 193–208, https://doi.org/10.1306/08250505059, 2006. a
Osorno, M., Uribe, D., Ruiz, O. E., and Steeb, H.: Finite difference calculations of permeability in large domains in a wide porosity range, Arch. Appl. Mech., 85, 1043–1054, https://doi.org/10.1007/s0041901510254, 2015. a, b
Patir, N. and Cheng, H. S.: An average flow model for determining effects of threedimensional roughness on partial hydrodynamic lubrication, J. Lubr. Tech., 100, 12–17, https://doi.org/10.1115/1.3453103, 1978. a, b
Reeves, D. M., Benson, D. A., and Meerschaert, M. M.: Transport of conservative solutes in simulated fracture networks: 1. Synthetic data generation, Water Resour. Res., 44, W05404, https://doi.org/10.1029/2007WR006069, 2008. a, b
Renshaw, C. E.: On the relationship between mechanical and hydraulic apertures in roughwalled fractures, J. Geophys. Res.Sol. Ea., 100, 24629–24636, https://doi.org/10.1029/95JB02159, 1995. a, b
Rodrigues, O.: Des lois géométriques qui régissent les déplacements d'un système solide dans l'espace: et de la variation des cordonnées provenant de ces déplacements considérés indépendamment des causes qui peuvent les produire, Journal de mathématiques pures et appliquées 1re série, tome 5, pp. 380–440, 1840. a
Ruf, J. C., Rust, K. A., and Engelder, T.: Investigating the effect of mechanical discontinuities on joint spacing, Tectonophysics, 295, 245–257, https://doi.org/10.1016/S00401951(98)001231, 1998. a
Rutqvist, J., Leung, C., Hoch, A., Wang, Y., and Wang, Z.: Linked multicontinuum and crack tensor approach for modeling of coupled geomechanics, fluid flow and transport in fractured rock, Journal of Rock Mechanics and Geotechnical Engineering, 5, 18–31, https://doi.org/10.1016/j.jrmge.2012.08.001, 2013. a
SKB: Data report for the safety assessment SRSite, Technical Report SKB TR1052, Technical Report TR1052, Swedish Nuclear Fuel and Waste Management Co., Stockholm, Sweden, 2010. a
Snow, D. T.: Anisotropie permeability of fractured media, Water Resour. Res., 5, 1273–1289, https://doi.org/10.1029/WR005i006p01273, 1969. a, b
Svensson, U.: A continuum representation of fracture networks. Part I: Method and basic test cases, J. Hydrol., 250, 170–186, https://doi.org/10.1016/S00221694(01)004358, 2001. a, b, c, d, e, f, g
Sweeney, M. R., Gable, C. W., Karra, S., Stauffer, P. H., Pawar, R. J., and Hyman, J. D.: Upscaled discrete fracture matrix model (UDFM): an octreerefined continuum representation of fractured porous media, Computat. Geosci., 24, 293–310, https://doi.org/10.1007/s10596019099219, 2020. a, b, c
Witherspoon, P. A., Wang, J. S. Y., Iwai, K., and Gale, J. E.: Validity of cubic law for fluid flow in a deformable rock fracture, Water Resour. Res., 16, 1016–1024, https://doi.org/10.1029/WR016i006p01016, 1980. a
Xu, C. and Dowd, P.: A new computer code for discrete fracture network modelling, Comput. Geosci., 36, 292–301, https://doi.org/10.1016/j.cageo.2009.05.012, 2010. a
Zhou, H., Li, L., and Jaime GómezHernández, J.: Threedimensional hydraulic conductivity upscaling in groundwater modeling, Comput. Geosci., 36, 1224–1235, https://doi.org/10.1016/j.cageo.2010.03.008, 2010. a, b
Zienkiewicz, O. C. and Taylor, R. L.: The Finite Element Method, 5th edn., ButterworthHeinemann, Oxford, 2000. a, b
Zimmerman, R. W. and Bodvarsson, G. S.: Hydraulic conductivity of rock fractures, Transport Porous Med., 23, 1–30, https://doi.org/10.1007/BF00145263, 1996. a, b, c, d
Zou, L., Jing, L., and Cvetkovic, V.: Modeling of flow and mixing in 3D roughwalled rock fracture intersections, Adv. Water Resour., 107, 1–9, https://doi.org/10.1016/j.advwatres.2017.06.003, 2017. a
 Abstract
 Introduction
 Fracture intersection flow modeling
 Permeability parameterization concepts
 Equivalent continuum representation of DFNs
 Fracturescale intersection flow benchmark
 Intersection flow effects at network scales
 Resolution dependency of ECM methods
 Discussion
 Conclusions
 Appendix A: ECMbased effective permeability prediction workflow
 Code availability
 Author contributions
 Competing interests
 Disclaimer
 Acknowledgements
 Financial support
 Review statement
 References
 Abstract
 Introduction
 Fracture intersection flow modeling
 Permeability parameterization concepts
 Equivalent continuum representation of DFNs
 Fracturescale intersection flow benchmark
 Intersection flow effects at network scales
 Resolution dependency of ECM methods
 Discussion
 Conclusions
 Appendix A: ECMbased effective permeability prediction workflow
 Code availability
 Author contributions
 Competing interests
 Disclaimer
 Acknowledgements
 Financial support
 Review statement
 References