Equivalent continuum-based upscaling of flow in discrete fracture networks: The fracture-and-pipe model

Abstract. Predicting effective permeabilities of fractured rock masses is a key component of reservoir modelling. This is often realized with the discrete fracture network (DFN) method, where single-phase incompressible fluid flow is modelled in discrete representations of individual fractures in a network. Depending on the overall number of fractures, this can result in significant computational costs. Equivalent continuum models (ECM) 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 this has the advantage of lower computational costs and the possibility to include matrix properties, choosing the right cell size for discretizing 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 termed intersection flow localization (IFL). We propose a new method to generate ECM's that includes this effect with a directional pipe flow parametrization: the fracture-and-pipe model. Our approach is tested by conducting resolution tests with a massively parallelized Darcy-flow solver, capable of representing the full permeability anisotropy for individual grid cells. The results suggest 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, ECM's are applicable to upscale flow in fracture networks, which reduces computational expenses for numerical permeability predictions of fractured rock masses. Furthermore, incorporating the off-diagonal terms of the individual permeability tensors into numerical simulations results in an improved representation of anisotropy in ECM's that was previously reserved for the DFN method.


flow simulations. So-called stair-case 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.
This study focuses on an often ignored but important aspect in fracture network modelling, that is given by the complexity of fracture intersection flow. To our knowledge, only few studies have presented 3D flow simulations within fracture intersections 70 (Zou et al., 2017;Li et al., 2020), revealing the fact that flow velocities will increase within the fracture intersections compared to the fractures 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 effective permeability should increase by a certain amount within the intersection. To demonstrate that, we systematically conduct 3D numerical simulations of Stokes flow within differently oriented, planar fracture crossings to analyse the permeability increase 75 caused by intersection flow localization (IFL). Using these results, we extend the current state-of-the-art methodology for equivalent continuum representations of DFN's to account for IFL in a quantitative manner and analyse its impact on effective permeability computations. There, it is still unclear, at which level of detail the ECM has to be discretized in order to conserve the structural complexity of the DFN, as aforementioned stair-case patterns and artificial connectivity cause resolution dependencies. Subsequently, resolution tests are performed on two DFN test-cases with a newly developed, massively paral-80 lelized and high-performance-computing (HPC) optimized finite element Darcy-flow solver, that is 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.

Fracture intersection flow modelling
Fluid flow in porous and fractured media is described by the well-known Navier-Stokes equations (Bear, 1972). It is commonly 85 assumed that sub-surface 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, isoviscous and the impact of gravity to be negligible, steady-state flow at constant temperature is defined by Stokes momentum balance (eq. 1) and continuity (eq. 2) equations (Bear, 1972): 90 1 and 2, utilizing PETSc (Balay et al., 2018) for HPC optimisation. 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 figure 1) while setting the other boundaries to no-slip (velocity component normal to the boundary is zero) enables the prediction of the mediums 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: with domain volume V. Using Darcy's law for flow through porous media (Darcy, 1856), that relates the specific discharge Q according to: with intrinsic permeability k and cross-sectional area A in combination with the fact that Q =vA, the directional permeability 105 k z is calculated by: As demonstrated by Eichheimer et al. (2019); Kottwitz et al. (2020); 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 (see comparison of errors to the result at largest 110 resolution in plot b, figure 5). 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 (plot 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 (plot 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.

Permeability parametrization concepts
As the two main structural features (fractures and intersections) composing a fracture network differ significantly in terms of their hydraulics (figure 1), they require independent concepts to parametrize their permeabilities for formulating their effective  grid block permeability tensor. For fractures, it is usual practice to use the cubic-law parametrization (e.g., Snow, 1969;Long et al., 1982), relating the specific discharge Q through a void system between two parallel plates according to: with the fractures 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: Natural fractures deviate from the assumptions of parallel plates, which is why a m in eq. 7 is commonly replaced with a 125 hydraulic aperture a h that corrects the parametrization 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 parametrization concept tailored for fracture intersections. The simulations shown on figure 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 t 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, cross-sectional 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 135 r h to compensate for the structural difference. As a first order approximation, we use half the size of the hypotenuse in a rightangled triangle whose legs are given by the two intersecting apertures (called half-hypotenuse assumption in the following, see figure 2 for details). This delivers sufficiently good results, as will be demonstrated later ( figure 6).

Equivalent continuum representation of DFN's
The use of the ECM approach instead of the DFN method to predict the effective permeabilities of fractured media crucially 140 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 145 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 (figure 1) into account.
Thus, we will explain our new approach to generate continuum representations based on DFN structures -the fracture-andpipe model.

150
Generally, the DFN approach offers a straightforward way to characterize structurally complex fracture networks. Most commonly, every fracture is modelled as a geometric primitive (here a disc) with a prescribed length l, center coordinate p 0 and unit normal vectorn 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ī between the two points defines its orientation. The goal of the ECM method is to generate a 3D regular grid with constant x-y-z 155 spacing δx, whereas every grid cell contains a symmetric, positive definite permeability tensor that is based on the fractures and their intersections.
To map each individual fracture to its corresponding grid cells, we first assume a horizontal disc (normal vectorḡ = [0, 0, 1]) at center point p g = (0, 0, 0) with corresponding fracture radius r (r = l/2) and represent it with an equally spaced set of points in the x-y plane P g , with the condition ||P g − p g )|| ≤ r. By that, we obtain a constantly spaced grid of points representing the 160 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 planē g = [0, 0, 1] with the actual normal vector of the fracturen. Utilizing Rodrigues's rotation formula (Rodrigues, 1840) around the rotation axis w = (ḡ ×n)/||(ḡ ×n)|| (unit vector orthogonal toḡ andn) yields the rotation matrix R f according to: 165 with ×, ·, and ||x|| denoting the cross-product, dot-product and vector norm of x, respectively. I represents the 3-by-3 identity matrix and C the cross-product matrix of the rotation axis 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 170 fracture in its actual 3D position: where * denotes matrix-matrix multiplication. By ensuring that the lower left corner coordinate of the rectangular grids bounding box is initially located at (0, 0, 0) (this may require a translation of all center points to incorporate all fractures), we obtain the grid-indices (i,j and k in x,y and z-direction, respectively) of the fracture by dividing P r with the cell size δx and round-175 ing the results. Finally, we compute the individual anisotropic permeability tensor K ijk for the cells by using a parametrized 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  (0.8003|1.000|0.0505) and its orientation is given by the unit vector (0.1270, 0.5839, −0.8018). The hydraulic pipe radius resulting from the half-hypotenuse assumption is 7.0711 · 10 −4 . b) visualizes the shape of the permeability tensor for an ECM model that considers only fracture permeability (grey, 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. c) presents the norm of the permeability tensor K ijk as a function of the ratio between initial point spacing δp and ECM grid spacing δx (see text for explanations). The dashed black line denotes the condition δp/δx ≥ 16, that is used to provide an correct approximation of the fracture and intersection volume per cell.
In case 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. To map all previously found intersections to the grid cells and formulate their permeability tensors, we follow the same workflow as presented for 185 individual fractures. A horizontal line of the same length as the intersection (||i 1 − i 0 ||), 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), whereas the mean point of the line is again located at (0, 0, 0). We then calculate the rotation matrix R i (eq. 10) by usingḡ = [1, 0, 0] andn = (i After identifying the corresponding grid i,j and k indices as described above, their permeability tensors are increased by using a parametrized interesection 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 πr 2 h , 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. plot 195 a) in figure 4), 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, that could cause an overall permeability increase by almost one 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 at the network needs to 200 be evaluated. To test the half-hypotenuse assumption (see figure 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 are calculated from high-resolutions Stokes-flow simulations (e.g., section 2) and compared to their analytically represents a sufficiently converged solution. This is done by calculating the L2-error-norm ||δ k || according to: whereas 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., accounting for variability around the mean. Fracture sizes l are distributed as a power law according to: 265 with a scaling factor β of 10 −4 . The only difference between the two test DFN's is the overall fracture number, which is 10000 for the DFN-A (plot a in figure 7) and 1000 for the DFN-B (b in figure 7), 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); Oda (1985). ||k ij || k xx k yy k zz k yx k zx k zy Figure 9. Absolute permeability values k for the 6 main components of the computed effective permeability tensor (principial 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 x-y-z direction).
behaviour (Maillot et al., 2016). if the ECM cell size is smaller than the individual fractures 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 parametrization we use as 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 here), the only computationally feasible solution is to use parametrization concepts (e.g. section 3). For that, we were able to demonstrate that the presented fracture-and-pipe ECM method is capable to provide 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 ECM's 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). There, the main problem 310 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 fractures 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, 315 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 DFN's 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 figure 9, 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 perme-320 ability overestimations. Only above a resolution of 128 cubic voxels all these effects seem to dampen out, allowing to declare the solution as 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 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 325 ECM cell size to hydraulic aperture should not exceed two to provide small flow errors. So far, the ratio of cell size to minimal hydraulic aperture in the system was much larger (about 1260) due to the low scaling factor β of the sub-linear aperture to length correlation (eq. 17). To achieve similar discretization ratios of Svensson (2001) while maintaining a power-law 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 330 rather recommend to stay above the maximum hydraulic aperture a h1 of the system, as otherwise the volume-fraction based permeability scaling factor in equations 13 and 14 exceed unity. In that case, parametrization assumptions might not hold any more, preventing the use of continuum flow methods. However, as demonstrated here, sticking to l 0 /3 > δx > a h1 as 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 an geometric upscaling procedure for 335 further effective permeability analysis. Notably, this range strongly depends on the structural character of the considered DFN, especially on the fracture size distribution and corresponding aperture correlation functions. For some DFN's this may require to crop the fracture size distributions from below to a few multiples of the cell size and compensate the hydraulic contribution of lower sized fractures with a background permeability.

Conclusions
In this study, we have analysed 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 similar to a pipe. This results in intersection flow localisation (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 (ECM) for effective permeability computations of discrete fracture networks (DFN) to incorporate IFL effects. Those are integrated by using a directional pipe-flow parametrization with a hydraulic radius of half the hypotenuse size in a right-angled triangle with side lengths of both intersecting hydraulic apertures. By assessing numerically the permeabilities of fracture intersections, we could demonstrate, that for system sizes close to the approximated pipe-radius (typically mm to cm), the effect of IFL on permeability can be almost one order of magnitude. At larger scales (system size of several hundred m) on the other hand, the impact of IFL on overall flow is minor. There, 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 to be lower than a third of the minimal fracture size and larger then 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 355 problems. Whether this holds for transport problems as well, needs to be determined in future studies.
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 im compact form by the following system of mass (eq. A1) and 360 momentum (eq. A2) conservation equations: whereas ∇ and ∇· denote the gradient and divergence operator for global 3D Cartesian coordinates, respectively. The specific 365 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 off-diagonal components k yx , k zx and k zy as well as fluid density ρ, gravitational acceleration g and fluid dynamic viscosity µ. We employ a 3D finite-element discretization 370 scheme (e.g., Hughes, 1987;Zienkiewicz and Taylor, 2000;Belytschko et al., 2000;Lin et al., 2014) for equations 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 equation 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 8-node rectangular elements (voxels) with linear interpolation functions (e.g., Zienkiewicz and Taylor, 2000) for volume integral approximation, whereas element integrals are evaluated by Gauss-Legendre quadrature rule (e.g., Belytschko et al., 2000) over 8 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:

385
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 nodal Darcy velocities u based on the newly solved nodal pressures by:

390
whereas the velocity vectors on the nodes are averaged from the neighbouring 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 395 boundary faces (e.g., lower and upper face in a cube for principal flow in z-direction) and the same linear interpolation between those two values is applied at the remaining four boundary faces (see figure A1 for an example). This ensures to capture both, the diagonal and off-diagonal terms of the permeability tensor properly, which are computed by substituting the volume averageū 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 (∆P z = δP δz ). The corresponding entries in the permeability Figure A1. Pressure boundary conditions for an applied gradient in z-direction. Here, top and bottom faces experience constant pressures of 1 and 0 P a, respectively. A linearly interpolated pressure distribution is applied at the remaining four boundary faces, as indicated by the coloured wedges next to the side-faces of the model. Thus, principal direction of flow is in z-direction, allowing to calculate the z-component related terms of the permeability tensor according to eq. A9 tensor are computed according to: and vice versa for the iterations with pressure gradients in x-and y-direction to obtain the permeability tensor as shown in eq.

A3.
The used single-continuum discretization scheme might appear simplistic compared to more sophisticated mesh-representations 405 (see Berre et al., 2019). However, the merits of our approach rather lay (1) on a fully anisotropic permeability representation of the individual continuum cells and (2) massive parallelisation and HPC optimization. Utilizing the parallelisation framework of PETSc (Balay et al., 2018) and their multigrid preconditioned solvers significantly reduces computational cost, allowing routinely simulations with 10 9 individual grid cells. An increase in grid resolution compensates the benefits of using conforming meshes or multi-continuum formulations (e.g., Berre et al., 2019). To test this, we compare our modelling procedure against in a matrix with a hydraulic conductivity of 10 −6 . The hydraulic conductivity in the grey band at the bottom is increased to 10 −6 . Constant pressures of 4 P a and 1 P a are applied at the inlet band (blue) and outlet band (red), respectively. The diagonal light grey line through the model indicates the sampling line for the pressures shown plot in b). There, 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 grey region illustrates the area between the 10th and 90th percentiles for the highest refinement level of the benchmarked methods, whereas the light grey region illustrates the same area for their lowest refinement level.
permeabilities and use the methodology described in section 3 to incorporate fracture permeability accordingly. The boundary 415 conditions are given by small pressure inlet (4 P a) and outlet (1 P a) bands as indicated in plot a in figure A2. The comparison of the pressure distribution (plot b in figure A2) highlights, that already with a resolution of 32 3 voxels, we obtain a good fit with the benchmark target field. This thus suggest that our modelling procedure is sufficiently correct for effective permeability predictions.
FE Darcy-Flow: https://bitbucket.org/mkottwitz/anisotropicdarcyupscaling/src/master/ ; commit: feaed524bcef636725ee30b2d4d3136b6525d83e Author contributions. MOK wrote the initial draft of the manuscript, performed numerical simulations, analysed the data and generated the figures. AAP supervised and helped designing the study, helped implementing the numerical methods and edited the manuscript. SA helped designing the study and edited the manuscript. BJPK edited the manuscript and discussed the results.