Stochastic 3 D Navier-Stokes Flow in Self-Affine Fracture Geometries Controlled by Anisotropy and Channeling

The investigation of fluid flow regimes in naturally fractured rock is crucial for assessing the performance of underground operations. While disposal of nuclear waste or tunnel construction aims for a tight subsurface (Tsang et al., 2015), highly permeable structures (e.g., fractures and faults) are essential for reservoir operations, such as enhanced oil recovery (Berkowitz, 2002) and geothermal energy production (Murphy et al., 1981). High flow velocities with Reynolds numbers (Re) >10 can occur in these reservoirs, especially in fractures and near-wellbore areas (Kohl et al., 1997; Liu et al., 2016). The flow regime is controlled by the complexity of the internal fracture morphology (i.e., roughness) and the interaction and intersection of different fractures (structural complexity). Therefore, hydraulic quantification is typically resolved using simplified numerical approaches, e.g., assuming smooth parallel plates and linear flow laws (Egert et al., 2020).


Plain Language Summary
The movement of fluid and heat through fractured rocks is massively affected by the void space of the contributing fractures and the velocity of the fluid. The opposing offset of the two fracture surfaces creates complex and variable void geometries. A calculation of the expected flow rates is very difficult. In most studies, simplifications are adopted neglecting these complex features by using a plane surface and assuming that the flow is taking place smoothly. We therefore simulate the real flow using the three-dimensional fracture void space and compare the results with simpler two-dimensional models. It can be demonstrated that the difference between the two types of simulations depends not only on the local simplification but also on the flow direction within the fractures and that the difference increases with higher flow velocities. If the fractures are not considered as a whole continuum but on the local scale, preferential fluid pathways or channels are formed in which most of the flow takes place. Depending on the flow direction, channels are more or less pronounced. We can show that in well-developed channels the differences between the calculation methods are much smaller than in the other parts of the domain. (Brown, 1989), while the effective aperture takes into account local surface tortuosity (Ge, 1997). Other authors proposed further modifications to account for different geometrical effects (Konzuk & Kueper, 2004;Wang et al., 2015). However, all LCL-based approaches still have two major limitations: (a) the fracture is simplified to a 2D geometry (Marchand et al., 2020), and (b) the Darcian-type LCL is assumed to be valid only for sub-laminar flow (Zimmerman et al., 2004). A critical Re ( crit Re ) for considering nonlinear effects is under discussion and spans from <1 to >10 (Cunningham et al., 2020). The Navier-Stokes equations (NSE) provide a general basis to overcome these limitations and cover complex hydraulic conditions (linear and nonlinear) in the fracture, but due to their high computational effort, they are mainly for single fractures (Javadi et al., 2014;. In the pre-stressed underground, the breaking up of intact rock and the subsequent shear displacement create highly variable fracture void spaces and introduce directionality into a previously nondirectional geometric problem. Interconnected pathways of larger voids can form, which enable a higher flow and associated advective transport of solutes and heat in channels in the direction of the overall pressure gradient (Klepikova et al., 2020;Moreno and Tsang 1994). Depending on the relative orientation of pressure gradient and stress field, these channeling effects can either promote or reduce the fluid flow (Lang et al., 2018). The ratio of the flow between the two extreme cases, where the flow is oriented perpendicular or parallel to the shearing, is called flow anisotropy (FA) (Auradou et al., 2006). Besides laboratory experiments (Gentier et al., 1997;Rasmuson & Neretnieks, 1986), FA was also investigated numerically (Koyama et al., 2006;Mallikamas & Rajaram, 2005;Marchand et al., 2020). Investigations were mostly limited to linear flow conditions, while nonlinear effects were only briefly addressed (Liu et al., 2020). To the authors' knowledge, no study has investigated the role of shear-dependent and flow rate-dependent channeling processes in a statistically relevant set of fractures and concerning the validity of the LCL. A quantitative analysis of the prevailing processes can only be performed in a spatially resolved manner.
The primary goal of this stochastic study is to assess the individual flow and surface effects causing the overestimation of the LCL and FA. Particular attention is paid to the formation of channels and their dynamic changes during the transition to nonlinear flow regimes. For this purpose, 840 individual realizations (30 fractures, 7 pressure gradients, 2 flow laws, and 2 directions) are simulated, and local differences between the fully 3D Navier-Stokes and 2D LCL solutions are quantified in terms of local flow rate and velocity. The significant number of individual fractures can be used to infer the range of uncertainty in large-scale reservoir models.

Rough Surface and Sheared Fracture Generation
A fracture surface can be created either by scanning a real surface or by numerically generating surface morphologies represented by fractal geometries following self-affine properties with isotropic correlation functions (Brown, 1995). In this study, a set of 30 random unique surfaces with a granite-typical Hurst roughness exponent of H = 0.8 (Schmittbuhl et al., 1993) are created using the workflow described by Méheust and Schmittbuhl (2001). The generated surfaces have a resolution of 128 × 128 nodes corresponding to a size of 20 × 20 cm to guarantee high accuracy with reasonable computational effort. The maximum surface amplitude is limited to 30 mm. Shearing is introduced by offsetting the fracture upper surface by one row of nodes in x-direction causing a sliding-up (Auradou et al., 2006). Maintaining at least one contact point results in a locally varying inter-surface void space (Figure 1a). In order to focus on the impact of the complex void geometry on channeling processes, mechanical effects, such as fracture frictional behavior or deformation of the void space, which can alter the local flow regime (Blöcher et al., 2019) and induce further uncertainties, are neglected.
While the resulting 3D fracture void geometry will be directly transformed into a 3D-mesh for the Navier-Stokes simulations, it is reduced to a 2D mesh for the LCL model including calculated local apertures. For the modeling presented here, an aperture distribution is used which has neither directional nor direct roughness dependences. Following the approach of Ge (1997), an effective aperture eff a is calculated, which is the distance between the bottom and top surfaces as surface normal to the local middle plane.

Governing Equations
Incompressible and steady-state fluid flow is governed by the NSE (Bear & Cheng, 2010): where ρ, P, u, f, and μ denote the fluid density, the total pressure, the velocity, and body force vector (e.g., the gravitational acceleration), and the dynamic fluid viscosity, respectively. For 3D problems, the equations form a set of four fully coupled partial differential equations solved for local pressure and velocity components. A simplified solution is to assume linear flow conditions and neglect the nonlinear inertia forces leading to the LCL (Nicholl et al., 1999): where i a is the local aperture

Numerical Modeling
The numerical simulations (LCL and NS) are carried out with the finite element open-source application TIGER (THMC simulator for Geoscience Research) (Gholami Korzani et al., 2020). The code is based on the MOOSE framework (Permann et al., 2020) and designed for solving thermo-hydraulic-mechanical-solute transport problems in geothermal reservoirs at different scales and dimensions. For the LCL simulations, an XY-extended FE mesh is used, which consists of four-node quadrilateral elements. The local effective aperture is represented as spatial-dependent material property. Nonlinear flow is realized via the Navier-Stokes module (Peterson et al., 2018). The simulations are performed using a structured 3D FE Mesh consisting of 127 × 127 × 10 hexahedrons created using Gmsh (Geuzaine & Remacle, 2017). Pressure-stabilized Petrov-Galerkin stabilization is applied to the equal-order elements for directly solving the PDEs. Top and EGERT ET AL.
10.1029/2020GL092138 3 of 10 bottom surfaces of the 30 unique fractures are transferred to a base mesh in such a way that it reflects their natural rough-walled and sheared fractures. A mesh sensitivity analysis was performed to avoid uncertainties regarding the discretization (see supporting information S2).
The further parameterization is done similarly for both kinds of simulations. A constant pressure of 0.01, 0.1, 1, 5, 10, 20, and 50 Pa is applied to the inlet boundary in the direction of flow, while the outlet boundary is fixed to zero pressure. Elsewhere no-slip boundary conditions are applied. Fluid density and viscosity are set constant at 10 3 kg.m −3 and 10 −3 Pa.s, respectively. All simulations are carried out for the two flow laws (LCL and NS) and the extreme cases of an overall pressure gradient parallel and perpendicular to the shearing, leading to a total of 840 realizations.

Fracture Flow
The fluid flow in the NS simulations is evaluated at each fractures' outlet to be able to consider both small-scale local disturbances as well as fracture-wide flow effects. The outlet flow rates NS Q are obtained from the mean velocities multiplied by the cross-sectional areas and are shown as function of the mean vertical fracture apertures a in Figure 2. Since the surface generation algorithm imposes a Gaussian roughness distribution the mean apertures a range between 1.2 and 2.5 mm with an accumulation of around 1.7 mm. For flow parallel to the shearing (red colors),

The Transition From Linear to Nonlinear Flow
The nonlinear effects along the fracture plane are analyzed by comparing the normalized flow rate NS normalized Q (Figure 3a): where ΔP is the global pressure gradient. Nonlinear effects arise at Re > 1, while at Re = 100, they already lead to a decrease of flow by 5%-10%, correlating well with the findings of Brush and Thomson (2003). The dispersion of the results increases strongly with increasing nonlinearity, especially for the parallel case.
In addition to the nonlinear effects, LCL and NS solutions already differ initially in laminar conditions (Figure 3b). For flow parallel to the shearing, the initial bulk overestimation of the LCL is about 9%-19%, which is consistent with previous observations (Al-Yaarubi et al., 2005;Nicholl et al., 1999). For flow perpendicular to the shearing, the mean overestimation is around 1%-9%. The significant differences between the NS simulations are the consequence of the extended travel distance between inlet and outlet (of more than 10%), which leads to reduced local pressure gradients and enhances the anisotropy compared to pure LCL investigations ( In Figures 2 and 3, the results are evaluated as the sum of the total flow through the fracture geometry, while Figure 1b depicts a fluid flow that is concentrated on small-scale fluid pathways. For further evaluation of the local validity of the LCL, we identify and separately analyze two flow domains, the channel sweeping domain, where the flow is concentrated in channels, and the surface sweeping domain, where a more dispersed flow is prevailing in off-channel areas.

Channel Identification
Channels are defined as continuous spatial features where large portions of the entire domain flux are concentrated. Various approaches have been developed to compare the extent of channeling using the LCLbased flow and effective permeability either deterministically (Watanabe et al., 2015) or statistically (Le Goc et al., 2010;Moreno & Tsang, 1994). However, given the linear relationship between Q and ΔP, these LCLbased studies cannot account for locally increasing nonlinear effects and resulting changes in the extent and localization of channeling within an individual fracture geometry. Within the scope of the present work, we tested several variables (e.g., aperture and velocity) and thresholds (e.g., 50% and 75% of the total flow) for their applicability to different fractures and pressure gradients. Inspired by the study of Knudby and Carrera (2005), flow paths are identified as channels where the upper quartile of the flow rate is localized. Subsequently, the effects of channeling are compared over a large number of realizations.
Based on the previously defined frameworks, a multi-stage post-processing concept is developed to identify and later quantify the individual channels within the studied realizations using MATLAB and the results of the NS simulations (Figure 1). Channel identification is an iterative process. All nodes are sorted in descending order in terms of their local flow rate. Starting from the highest flow rate, channel nodes are marked and their nodal flow rates are added up until 25% of the mesh's entire flow is localized within channels. Subsequently, the skeletonized channel network is extracted and further evaluated for its connectivity. The connectivity is checked by setting up a 3 × 5 nodes sampling window in the direction of the pressure gradient, which assigns nodes with two or more connections to a unique channel. Since channels are defined as a coherent and continuous structure, background noise, and poorly connected nodes (node row of less than 10% of the mesh side-length) are removed (Figure 4).
EGERT ET AL.    (Figure 4b) to shearing. A perpendicular orientation of shearing and flow leads to strong localization and formation of well-connected and continuous channels. In contrast, flow parallel to the shearing causes less localization and poorer connectivity characterized by high velocities in these regions and significantly increased flow impedance (Auradou et al., 2006). The resulting FA is more pronounced for smaller a, as the probability of well-connected channels in the parallel extreme case increases with larger a.

Spatial Assessment of the LCL Accuracy
The LCL with a local effective aperture is one of the most common flow laws adopted for computing fractured rock hydraulics. Nevertheless, its applicability is often determined at the fractures' outlet, leading to the conservative conclusion that the LCL is only valid for Re ≪ 1. However, the presented analyses show that the local Re can vary within one order of magnitude for a single realization depending on whether the local Re is calculated within or off a channel. Therefore, the global error is broken down into a local error, which is mathematically described by: x y are the volumetric flow rate calculated with the NSE and LCL, respectively. The volumetric flow rate is obtained from the mean velocity multiplied by the crossectional area.
In Figures 5a and 5b, the relative errors are shown as a function of the local Re for one individual realization considering flow parallel and perpendicular to the shearing. The mean error increases nonlinearly toward low flow rates and small apertures, indicating that the aperture correction according to Ge (1997) has a higher validity in channeled structures, while surface effects cause significantly increased local errors. In the identified channels (high Q and Re) the local error is negligibly small so that LCL and NS predict the same local flow rate. The small-scale results show a slight directional dependence on the overall trends. The mean error for flow perpendicular to the shearing is reduced and single nodes even show an underestimation by the LCL.  the differences between NS and LCL are quantified in terms of their uncertainties and variability with respect to the tectonic setting. While the differences in the identified channels remain constant, it increases significantly for the surface sweeping areas with increasing pressure gradients, and the stronger local nonlinearity results in an increased overall error (Figures 5c and 5d). Channeling thus becomes more pronounced with increasing flow rate. Since inertial effects are prevalent for surface sweeping, less void space is required to localize 25% of the entire domain flow. For higher Q and Re, the overall results tend to show higher variability and standard deviation.
The described differences, both purely laminar and incipient nonlinearity, can be traced back to the internal velocity field. In the channel sweeping, a parabolic shape is maintained, i.e., the basic assumption of the LCL, while it deviates significantly from this ideal shape for the surface sweeping domains. Here, complex internal void geometries such as contact spots, asperities, and local roughness cause transverse flows, eddy currents, or stagnant water zones that simply cannot be considered in 2D fracture models (Zou et al., 2017). Cunningham et al. (2020) showed experimentally on a mated fracture that the crit Re to consider nonlinear effects decreases with increasing local complexity of the fracture void geometry. Thus, both decreasing apertures and an increase in relative aperture fluctuations cause enhanced inertial effects. Our realizations confirm the presented relationship between crit Re and relative surface roughness. Transferred to sheared fractures and resulting anisotropic flow, a parallel orientation of shearing and flow favors the formation of irregular and highly tortuous channels with increased viscous friction. This process causes both, a decrease in the observed flow rate and an increased sensitivity to nonlinear effects (reduced EGERT ET AL.  Orange indicates channels while gray represents the error located in the surface regions. The identified channels generally tend to have less difference between local cubic law (LCL) and Navier-Stokes (NS) and fewer nonlinear effects (values closer to 0). crit Re ). In contrast, with a perpendicular orientation of shearing and flow, well-developed and smooth channels result in reduced geometry dependence with significantly higher LCL accuracy and reduced nonlinear effects at the same flow rate. In such a complex environment (in terms of individual geometry and internal flow field), it is therefore difficult to introduce a generally valid correction factor for the LCL and growing nonlinearities, as described for example by Wang et al. (2020) or used in the well-known Forchheimer equation (Chen et al., 2015). However, with respect to geothermal reservoirs, this also implies that a typical normal stress regime, with vertically oriented and displaced fracture surfaces, favors the formation of channel sweeping for horizontal flow between two wells (Lang et al., 2018). Strong channeling increases the range of validity of the LCL at higher Re (<10% deviation for Re ∼ 100) but could cause reduced effective heat exchange within the reservoir (Guo et al., 2016). An additional strike-slip component would however lead to more dispersed surface sweeping and thus reduces the reliability of LCL predictions.

Conclusions
Originally focused on the impact of higher Re on flow through self-affine fracture geometries under various shear displacements, our stochastic analysis of 3D Navier-Stokes flow reveals the importance of local complexity on the spatial validity of the LCL and explains contradicting observations of different studies. In fracture void geometries, the definition of a single critical Re (e.g., crit Re = 1) is not sufficient to describe the flow regime and associated processes; an additional spatial analysis is crucial. The initial Gaussian aperture distribution is not reflected in the local flow regimes; instead, the character of the flow could vary between the two extreme cases, a rather focused channel sweeping or a rather dispersed surface sweeping, resulting from the orientation of flow and shearing. Channel sweeping, typically favored by flow perpendicular to the shearing, is characterized by highly localized channels with an ideal parabolic velocity field. A weak dependence on the relative roughness and tortuosity of the fracture results in accurate LCL prediction (<10% overestimation) and mostly linear flow occurring even at Re ∼ 100. In contrast, local geometrical complexity in the dispersed surface sweeping entails increased frictional and inertial effects resulting in increasing differences between NSE and LCL (>15%). The small-scale stochastic nature of our study allows us to infer the range of uncertainty in large-scale fluid flow phenomena.