Spatial Characterization of Channeling in Sheared Rough‐Walled Fractures in the Transition to Nonlinear Fluid Flow

Accurate quantification of spatially resolved fluid flow within fractures is crucial for successful reservoir development, such as Enhanced Geothermal Systems. This study presents an innovative workflow designed to model and characterize preferential flow paths (channels) within rough‐walled shear fractures. A set of 30 rough‐walled self‐affine fractures, all possessing identical roughness characteristics, is stochastically generated. By solving the nonlinear Navier‐Stokes equations in 420 individual realizations, the transition from linear to nonlinear flow regimes and the two extreme flow directions perpendicular and parallel to the shearing are numerically captured. A distinguishing feature of this approach is its comprehensive statistical analysis, which encompasses both the geometric and transport properties of flow paths in the non‐simplified three‐dimensional fractured void space under typical geothermal flow conditions. In a perpendicular orientation of flow and shearing, fluid flow exhibits pronounced localization, with more than one‐third of the volumetric flow concentrated within 15% of the fracture volume. In contrast, parallel to the shearing, a complex pattern of individual tortuous channels emerges, with flow occurring in 22% of the void space. Nonlinear effects primarily manifest outside these channels, suggesting that complex flow phenomena may dominate irregular fracture structures, such as contact zones or asperities. In the parallel case, increased flow rates lead to an amplification of channeling processes resulting in less affected volume and diminished tortuosity of the main flow path, while in perpendicular orientation nonlinear effects are only of minor importance. The small‐scale flow regime of both extreme cases tends to converge with increasing flow rates.

Besides fracture geometry, the fluid flow at reservoir scales is often represented by highly simplified assumptions.The most common approach, cubic law (CL), assumes fractures as smooth parallel plates (Witherspoon et al., 1980).The local cubic law (LCL) accounts for local void space changes but depends strongly on the aperture calculation (Cheng et al., 2020;L. Wang et al., 2015;Zimmerman & Bodvarsson, 1996).All LCL-based models introduce two types of uncertainties.First, the fracture geometry is projected onto a 2D plane, causing an underestimation of geometrically induced flow effects.Second, the Darcian-type LCL is only valid for laminar flow and low Reynolds numbers (Re « 1) (Brush & Thomson, 2003).Both limitations pose problems for the analysis of spatially resolved fluid flow in fractured reservoirs, for example, during geothermal energy production and hydraulic stimulation (Chen et al., 2015;Kohl et al., 1997;Louis, 1967).The solution of the nonlinear Navier-Stokes equations (NSE) on rough fracture surfaces allows to circumvent these limitations, but they are rarely used due to high computational complexity (Javadi et al., 2014;Liu, He et al., 2020) and are mostly limited to small Re (Aghajannezhad & Sellier, 2022;He et al., 2022).Egert et al. (2021) were able to demonstrate the need for solving the NSE to consider both the three-dimensional fracture void space and growing nonlinear effects in the small-scale fracture-internal flow field by comparing NSE and LCL for Reynolds numbers up to 220.The natural small-scale flow field in a fracture is not homogeneously distributed.Preferential flow paths are formed locally, leading to a spatially varying agreement between the two flow laws with increasing flow velocity (Egert, 2021) and turbulence already at moderate flow rates (Nikuadse, 1930).
These preferential fluid pathways are the results of tectonic stresses in the earth's crust leading to the formation of new fractures as well as the shearing of existing opposing surfaces.This creates more or less connected void spaces both within fracture networks and along with individual fractures.In these void spaces, fluids tend to flow in the direction of least resistance (Moreno & Tsang, 1994).The effect is generally known as channeling and is defined as the spatial concentration of fluid flow along preferential hydraulic pathways through geologic systems with heterogeneous internal structures (Tsang & Tsang, 1987).It has been observed in small-scale laboratory experiments on individual surfaces (Auradou et al., 2006;C. Wang et al., 2022) and in field-scale fracture networks (Guihéneuf et al., 2017;Tsang & Neretnieks, 1998).As a result, both geometric and transport properties are no longer generally valid for the entire fracture, but anomalous transport properties are highly affected by the geometric nature (e.g., contact areas and local permeability) of these individual channels (Talon et al., 2012;Z. Wang et al., 2021).This concerns not only the channel length/width but also the fracture area/volume covered by actively flowing water (flow-wetted surface, flow-through volume) and solute mixing processes (Figueiredo et al., 2016).Different approaches have been developed to characterize the extent of channeling experimentally (Klepikova et al., 2020) or numerically (Watanabe et al., 2015).
Numerically, these approaches focus either on the geometric and statistical description of the channeling processes or the effects on the flow-related transport in various spatial scales (Renard & Allard, 2013).Talon et al. (2010) showed that the transport properties depend decisively on the worst permeability within the channels, the so-called bottleneck.Ronayne et al. (2008) inverted the statistical properties of an aquifer system considering channeling.Maillot et al. (2016) evaluated the connectivity of channels in a discrete fracture network.Le Goc et al. (2010) use effective permeabilities and momentum analysis for the stochastic description of channel distributions on different mathematically generated aperture fields and 2D discrete fracture networks.Knudby and Carrera (2005) compared several statistical-, flow-and transport-related indicators for their applicability to identify the connectivity of channels on a single surface.Bruderer-Weng et al. (2004) analyzed the multifractal spectrum of the flow field to quantify the effect of channeling on the dispersion of tracer.Western et al. (2001) analyzed the connectivity of flow channels and concluded that classical indicator geostatistics like variograms cannot distinguish between random and connected soil moisture patterns like channels.Marchand et al. (2020) evaluated the geometrical properties of channels concerning shear displacements of the opposing fracture surfaces.Tang et al. (2019) evaluated the channeling characteristics and transport-related properties in pore-scale heterogeneous media.Javanmard et al. (2022) compared different channel connectivity indicators regarding their applicability under normal stress conditions.In this way, several different indicators and approaches have been investigated, with all studies being limited to simple and planar 2D geometries, small flow velocities, and/or linear flow laws.Siena et al. (2019) demonstrated for 3D porous media that increasing Reynolds numbers leads to more focused flow in individual channels.Liu, Wang et al. (2020) could reveal a correlation between stress-induced local aperture, pressure gradient, and the onset of channeling for a single fracture.The importance of nonlinear and 3D geometric effects concerning channeling has been addressed but has not been further quantified.
One outcome of these preferential fluid pathways is the emergence of flow anisotropy.Flow anisotropy refers to the difference in the volume of fluid flowing through a fracture when it flows in two extreme directions: parallel or perpendicular to the shearing direction (Auradou et al., 2001).Besides laboratory experiments (Gentier et al., 1997;Rasmuson & Neretnieks, 1986), flow anisotropy was also investigated numerically.The direct link between the effect of channeling and flow anisotropy could be confirmed in small-scale experiments (Silliman, 1989) and numerically at the reservoir scale concerning the local stress field (Egert et al., 2021;Lang et al., 2018).However, it has not yet been investigated how the geometric properties (e.g., volume, area, connectivity) of these channels can be quantified in the three-dimensional void geometry and how these channel characteristics change during the transition from linear to nonlinear flow conditions and with respect to the extreme cases described as flow anisotropy.
This study focuses on two main objectives.First, a novel and practically applicable workflow concerning the quantification of small-scale fluid pathways in the unprojected 3D fracture void space is developed.Then, different indicators are employed to quantify and discuss these identified channels with respect to changes in the extent of channeling under different flow and geometric conditions, such as increasingly nonlinear and directional (anisotropic) flow.In the first step, a variety of unique rough surfaces are created and subsequently sheared to generate individual fractures.Then, the local flow field within the three-dimensional fractures is determined by directly solving the NSE considering different flow boundary conditions (BC) and directions.Using a novel multi-stage workflow, the individual channels are identified, and suitable indicators are introduced to portray their geometric features as well as flow and transport properties in the transition from linear to nonlinear flow regimes.Finally, this workflow, applied across a substantial quantity of unique small-scale realizations allows a transfer of the knowledge on the effective properties, such as heat and mass transfer, to fractured geothermal reservoirs.

Sheared Fracture Generation
A large number of statistically generated and unique fracture geometries are required for the probabilistic determination of the different channel characteristics.An individual surface can be either created from scanning data of a real fracture surface or by numerically generating surfaces following mathematical functions (Brown, 1995).In this study, it is assumed that the fracture surfaces can be represented by fractal geometries that follow self-affine properties with isotropic correlation functions (Méheust & Schmittbuhl, 2003).The probability density function p d is defined as (Equation 1): where λ is an arbitrary scaling factor, H is the Hurst roughness exponent, ∆h and ∆r are the difference in height and distance, respectively.Every generated surface is unique, as each is based on a distribution extracted from a white noise generator.A roughness exponent of H = 0.8 is used (Figure 1c), which is often observed for natural fractures in granitic rock (Schmittbuhl et al., 1993).The generated rough surfaces have a resolution (m × n) of 128 × 128 nodes resulting in a size of 20 × 20 cm and a total node number N of 16'384.The maximum surface amplitude is limited to 3 cm.A total of 30 different surfaces are generated and analyzed to obtain statistically reliable results.
In the geothermal context, freshly fractured surfaces with well-connected void space and good hydraulic conductivity are of interest.These voids are created when initially well-mated fracture surfaces are sheared against each other under natural or artificially induced tectonical stresses (Egert et al., 2018).To introduce this shearing into the numerical models, the rough surface is doubled.The upper surface is raised and laterally displaced by one row of nodes.The upper surface is lowered until at least one contact point between the two opposing surfaces is attained, leading to a locally varying void space.The shear displacement is considered small compared to the entire fracture length (<1%).No deformation or rotation is assumed, as the study focuses on the nonlinear effects on the channel formation within a developing fracture in the space of rigid approximation (Figure 1a).
The resulting three-dimensional rough-walled fracture surfaces serve as the foundation for adjusting the numerical finite element (FE) mesh to solve the NS equations.However, to translate the resulting flow velocities (the NS equations are solved for the pressure and velocity vector) into a volumetric flow rate (represented as Q), we require the cross-sectional area perpendicular to the pressure gradient.This is achieved by multiplying the constant grid spacing by the local aperture, as illustrated in Figure 1b.A vertical or mechanical aperture a vert is not useful for rough fractures and leads to an overestimation of the volumetric flow rates (Cheng et al., 2020;L. Wang et al., 2015).We therefore use the effective aperture a eff according to Ge (1997).This method involves determining a "true aperture" perpendicular to the local orientation of the centerline between the two rough fracture surfaces.This approach has been shown in previous studies (Marchand et al., 2020) to be a good approximation for the local volumetric flow rate in rough fractures.To gain a eff , the local median m(x,y) is established and the median curve has to be calculated as follows: where h bot (x,y) and h top (x,y) denote the height of the local bottom and top surface (Figure 1b).Subsequently, the slope of the local median plane grad (m) is calculated (Equation 3) and the normal vector n (m) of the center line can be obtained (Equation 4): (3) 10.1029/2022WR034362 5 of 17 The effective aperture a eff is defined as the distance between h bot and h top intersecting the point m(x,y) and using the normal vector n(m): where B and T denote the intersection of the bottom and top surfaces with n(m).a eff ranges between 0 and 3 mm and is typically slightly reduced compared to the vertical aperture (Figure 1d).The mean aperture  eff is calculated to compare across the different meshes and to compute the variability in the local aperture distribution.Figure 2 shows that  eff ranges between 1.2 and 2 mm, with some outliers to higher mean apertures.

Flow Governing Equations
The approach assumes an equivalent grid block for the fracture void space, where fluid flow and transport processes can take place.Incompressible and steady-state fluid flow is governed by the NSE (Equation 6) and the mass conservation equation (Equation 7) (Bear & Cheng, 2010): where ρ is the fluid density, P is total pressure, u denotes the velocity vector (u, v, w) and μ is dynamic fluid viscosity.
For three-dimensional problems, the equations form a set of four fully coupled partial differential equations governing the velocity components and the local pressure.Simpler solutions of the NSE such as neglecting the inertia term (e.g., the Stokes equation or the LCL) are not permissible for this study due to high Re > 1 (Brush & Thomson, 2003).

Numerical Modeling
The generated fracture surfaces are directly used for the FE modeling (FEM) by projecting them onto a three-dimensional FE mesh.For this purpose, we create a cuboid structured mesh with the dimensions of 128 × 128 × 11 comprising equally sized hexahedra, using the software Gmsh (Geuzaine & Remacle, 2017).Subsequently, the vertically superimposed nodes of the original mesh are deformed in height so that they align with the upper and lower fracture surfaces created earlier.The mesh nodes located between these upper and lower surfaces are evenly spaced in the x-y direction, mirroring the characteristics of the rough surfaces.Their vertical positions remain constant and are determined by the local height, as shown in Figure 1b.
The numerical simulations are carried out with a FE open-source application called TIGER (TMHC sImulator for Geoscience Research) (Gholami Korzani et al., 2020).The code is based on the MOOSE (Multiphysics Object-Oriented Simulation Environment) framework (Lindsay et al., 2022) and is designed for solving thermo-hydro-mechanical-solute transport problems in geothermal applications on different spatial and temporal scales.Nonlinear flow is realized via the Navier-Stokes module (Peterson et al., 2018).Each 3D realization is solved in parallel on an HPC system using Newton's method combined with the parallel direct solver SuperLU_ dist at each linear iteration (X. S. Li & Demmel, 2003).Pressure-Stabilized Petrov Galerkin (PSPG) method is introduced to avoid instabilities arising from unstable element pairs and allow the usage of equal-order discretization in incompressible Navier-Stokes modeling (Hughes et al., 1986).The maximum mesh size is determined by a mesh sensitivity analysis to achieve highly accurate convergence with reasonable computational effort.
Uniformly refining an initial mesh with 16 × 16 × 1 elements did not lead to any improvement in local volume flow beyond the third refinement step (vertical resolution of eight elements), as shown in Figure 3.
Dirichlet-type BC are applied in the direction of the overall pressure gradient (parallel and perpendicular to the shearing), with p = 0.01, 0.1, 1, 5, 10, 20, and 50 Pa at the inflow and p = 0 at the outflow boundary.The

Channel and Connectivity Identification Workflow
In general, channeling is defined as the spatial concentration of volume flow.However, no universal criterion or threshold exists neither for its identification nor the quantification of the associated processes.Consequently, commencing from the 3D Navier-Stokes flow simulations, a multi-step image sampling workflow (Figure 4a) is developed.This workflow allows us to pinpoint individual channels within a specific flow-through medium, and the selection of input variables is customized to the nature of the problem at hand.The developed workflow was inspired by the 3D work of Ebrahimi et al. (2014) and the 2D work of Le Goc et al. (2010) and was first presented in Egert (2021) and Egert et al. (2021).Related modifications were later presented in Javanmard et al. (2022) with respect to percolation theory.The individual steps of the workflow, from the fracture generation (Figure 4b) to the identified channels, are vividly illustrated using a base case, a 3D NS realization, and the volume flow as a critical variable (Figure 4).
In the first step, the variable of interest (e.g., localized volumetric flow Q, Figure 4c) is identified, and the first guess for a threshold, here a local volumetric flow rate threshold Q t , is made.For all nodes N, each node i with a volumetric flow rate (Q i ) greater than or equal to the threshold value (   ≥  ) is categorized as a channel node.Subsequently, we assess whether these marked nodes satisfy the primary criterion, which, for example, may involve verifying if 35% of the entire mesh's volumetric flow rate is concentrated within these marked nodes (as depicted in Figure 4d).If this is not the case, the threshold value Q t is continuously reduced, and the workflow is iterated again over the entire mesh.If true, the skeletonized channel network is extracted, and the connectivity between the identified channel nodes is further evaluated.This is done by setting up a 3 × 5 nodes sampling window, which assigns nodes with two or more connections to the same unique channel (Figure 4e).Given that a channel is characterized by its coherent and uninterrupted structure, the subsequent phase involves the elimination of background noise in the form of isolated or weakly connected channels, which occupy less than 10% of the total fracture length (l frac , Figure 4f).As suggested for example, in Western et al. (2001) it is currently not possible to distinguish between randomly short and patterned structures in the skeletonized channel network.Therefore, in the final step, the longest ongoing channel is identified, and the length of the geometrically tortuous flow path is evaluated using Dijkstra (1959) shortest path algorithm.This path length is the shortest connection along the identified channel and considers all geometrical effects, including both the local roughness and the geometric tortuosity caused by the undulating fracture surfaces (Figure 4g).
Once an identification workflow is established, a critical variable must be identified to detect and quantify channeling processes in terms of flow geometry and regime.While this criterion should be rather simple and general, it must allow the identification of the flow backbone.Although it is common to use the variable of main interest, the volumetric flow rate (or the recalculated effective permeability), as a critical variable (Le Goc et al., 2010), several other variables (aperture, velocity, flux, pressure gradient, participation ratio) have been tested for their applicability to different flow constraints.Figure 5a shows the identified channels exemplified for a single realization with a boundary pressure gradient of 50 Pa, a threshold ratio of 25%, and different tested variables.
The local volumetric flow rate Q i , as a multiplication of flux and aperture, provides the most promising results for the determination of channel characteristics and is therefore retained.While a minimum volumetric flow threshold Q t for fracture networks is obvious (is the flow percolating along a fracture?, Maillot et al., 2016), a deterministic absolute value for the local volumetric flow rate threshold Q t in combination with single fractures or porous media is not useful.We are seeking an approach that is as universal as possible in all types of flow-through media and that allows comparability between all 420 different realizations and pressure gradients.The latter can only be achieved by normalizing the individual realization-dependent volumetric flow rates and keeping the ratio of normalized flow localized in the channels (R chan ) constant (Equation 8): where N is the total number of elements, Q i is the volumetric flow rate vector at the ith element of the fracture, Q chan the amount of volume flow localized in channels, Q tot the sum of the entire domain volume flow.

Geometric Threshold Identification
Once the critical variable has been determined, an appropriate threshold must be identified.As introduced in the previous chapter, there is no common   crit chan to serve as a threshold for the identification of channeling.One approach describes the use of a specific percentile of the nodal volumetric flow rate as a critical limit for channel identification (Knudby & Carrera, 2005;Marchand et al., 2020).Such a criterion is quite simple and obvious but depends on the arbitrary choice of the threshold.We seek a statistical characteristic that captures the maximum extent of channeling and, at the same time, can be used to compare channel properties with the onset of nonlinear flow conditions and the resulting changes in flow characteristics.Inspired by the rather statistical indicators of Le Goc et al. (2010) and the application-related indicators of Marchand et al. (2020), we discuss the two geometric characteristics of interconnected transport volume and heat exchange area with increasing channeling threshold R chan .Subsequently, we define a characteristic   crit chan to quantify channeling.Addressing this question, both properties must be defined first.The effective transport volume V1 of the identified channels is the sum of the volume of all previously identified channel elements compared to the entire fracture void space.V1 as a function of the channeled flow R chan is determined by Equation 11: where V1 is the normalized channel volume, V tot is the entire fracture volume, V i is the element volume.Similarly, the heat exchange area or more generally, the flow-wetted surface A1 can be defined as the sum of the surface area (top and bottom surface) of all identified channel elements compared to the entire fracture surface area.A1 is therefore affected by both local aperture and height differences.A1 as a function of the channeled flow R chan is determined by Equation 12and considers, that a fracture is confined by two opposing surfaces: where A1 is the normalized channel half area with respect to Q t .A tot ,   top  and   bot  are the entire 3D fracture surface area, the element top area, and the element bottom area, respectively.
In the case of homogeneous flow through the entire fracture (e.g., the parallel plate approach and valid CL), a linear relationship exists between the flow ratio R chan and the channeled volume/area with always V1 = V pp and A1 = A pp .1 and A1 are biased toward smaller values when the volumetric flow rate contrast within the fracture increases.With increasing channeling, the distance between the identified flow paths increases.Thus, not only the absolute values of V1 and A1 are of interest, but also the distance to the parallel plate solution is an indicator of the amount of channeling expressed as V pp − V1 and A pp − A1. Figure 6a shows V pp − V1 as a function of R chan for all realizations with a constant pressure gradient of 50 Pa and for both two extreme cases of the orientation of flow and shearing oriented parallel (red) and perpendicular (blue).
All 30 realizations, regardless of the flow direction, show a strongly localized flow pattern, where the effective flow-through volume represents only a small portion of the total fractures' void space (V1).However, notable distinctions exist between the two extreme cases.In a perpendicular orientation of flow and shearing, substantial portions of the fluid become concentrated within a few extensively developed channels.The results show strong sensitivity to the initial fracture geometry.This extreme case, characterized by highly localized flow and minimal volume involvement, is therefore defined as channel sweeping.A parallel orientation of flow and shearing leads to a significantly lower degree of channeling, and large portions of the fluid are distributed over large fracture volumes.The outcomes are uniform across all fractures and show no geometric dependence.This case is therefore defined as surface sweeping.
Similarly, Figure 6b shows A pp − A1 as a function of R chan for all realizations with a constant pressure gradient of 50 Pa and the two extreme cases of the orientation of flow and shearing oriented parallel (red) and perpendicular (blue).Localized flow tends to adhere to the regular geometries that offer the least resistance.Tortuous regions within the entire fracture with large surface areas predominantly fall within the surface sweeping domain.Contrasts between flow directions are less pronounced than for V1.As a result, when channeling intensifies, there is a substantial reduction in the relative flow-wetted surface area.
We now return to the question of a suitable threshold value for   crit chan for our analysis of nonlinear channeling effects.Since R chan is not uniform but varies between 0 and 1 (see Figure 6), we seek a statistically critical   crit chan .This value should be sufficiently high to effectively capture channeling for our analysis while allowing for the evaluation of discrete channels and alterations in their connectivity.In this regard, several authors (e.g., Knudby & Carrera, 2005;Marchand et al., 2020or Western et al., 2001) have discussed the optimal choice of volumetric flow rate to study locally resolved channeling processes targeting mainly 2D fractures and the LCL.However, our study (Figures 5 and 6a) reveals that the third quartile of elemental volumetric flow rate (   crit chan = 0.25 ) proposed in these studies is not optimally suited for our analysis of 3D fractures to fully capture changes in channeling due to nonlinear flow effects.Instead, the analysis (e.g., Figure 6) shows that for 13) maximum channeling and the onset of percolation are best captured across all geometries and for a parallel orientation of flow and shearing.The derived threshold is therefore equally applied in the same way to all realizations to ensure comparability in the subsequent investigation of nonlinear channeling characteristics.

Nonlinear Channel Characteristics
The channeling characteristics of all realizations are evaluated based on the previously discussed critical flow ratio   crit chan = 0.33 .The use of seven pressure gradients between 0.05 and 250 Pa•m −1 ensures a range of mean Re « 1 up to 220 (Egert et al., 2021), which allows the quantification of nonlinear effects on channeling processes in the transition from linear to nonlinear flow regimes.In the first step, the changes in the geometrical indicators V1 and A1 are presented.Then two indicators describing the channel continuity and interconnection (C1 and C2) are further investigated and a final summary is given (Table 1).not alter the mean value.In contrast, surface sweeping dominates with a parallel orientation of flow and shearing.

Geometry Indicators
The results exhibit an initial volume between 20% and 23%, which is rather independent of the particular fracture geometry.Upon increasing pressure gradient, nonlinear effects in the surface sweeping areas cause increased localization of flow within channels resulting in a reduction in channel volume in the nonlinear flow region (−1% for a pressure gradient of 250 Pa•m −1 ).Siena et al. (2019) showed comparable trends with shrinking channel size in porous media and increasing Re in the transitional regime (Re between 1 and 10).Egert et al. (2021) were able to show that the nonlinear effects are significantly more pronounced in the surface sweeping regions outside the identified flow channels, while in the channels an ideal flow field and a higher LCL validity tend to be maintained even at high flow rates.Indicator A1 has similar trends as V1 (Figure 7b).A perpendicular orientation of flow and shearing results in a flow-wetted surface between 12% and 18%, which is only slightly affected by increasing flow rates.The mean area remains constant, but the variance of the results increases.Once flow and shearing are parallel, values for A1 between 17% and 21% are observed.Nonlinear effects cause a reduction of A1 from 19% to 18% and the more pronounced channeling leads to less dispersion in the results.

Continuity Indicators
The two previously introduced indicators V1 and A1 describe the channeling process as a whole, but no statement can be made about the distribution of the channels within the fracture domain and the connectivity and continuity of the individual channels concerning the flow regime The interconnection and connectivity of the channel network are further evaluated by identifying the main flow path within each realization (Figure 4g).Using the shortest path theory, the two flow-based indicators C1 and C2 are developed to describe the flow path in the unprojected three-dimensional space.The combined interconnection and continuity indicator C1 is the normalized channel length for the entire fracture length (Equation 14): where C1 is the normalized distance dist between the start point p start and endpoint p end of the main channel and l frac is the fracture length in the direction of the global pressure gradient.C1 ≥ 1 indicates channels that can be traced in the direction of the pressure gradient throughout the fracture and percolation is achieved.A maximum value of  C1 = √ 2 is reached when the main flow path is diagonally connected through the fracture.
The continuity indicator C2 is not a measure of the absolute channel length within the fracture but describes the relative three-dimensional continuity of the main flow path within each realization following the fracture void space (M.Wang et al., 2016).Therefore, 3D channel tortuosity C2 is determined as the shortest connection within the extreme points and following the identified main channel (red line in Figure 4g) and calculated using Equation 15: where C2 is the normalized and path(p start , p end ) the absolute length of the 3D tortuous flow path along the local median plane m(x,y) and in the direction of the overall pressure gradient.Indicator C2 can take values greater than 1, where C2 = 1 means that the channel suffered no topological height changes.A high C2 is a measure of the reduced pressure gradient within the fracture.
Figure 8a depicts the resulting continuity indicator C1 concerning the global pressure gradient and for flow and shearing oriented perpendicular (blue) and parallel (red).The averages and standard deviations of the statistical indicator C1 are shown in Table 1.A perpendicular orientation of flow and shearing results in well-formed and straight channels that exhibit continuity throughout the fracture with  C1 =1 , regardless of the initial sheared geometry.This implies that the channels are close to percolation and straight through the fracture or that the channels are oriented diagonally and slightly shortened.The transition to nonlinear flow regimes results in increased channeling and an associated increase in mean continuity for the given threshold.In the case of parallel orientation of flow and shearing, C1 exhibits a strong dependence on the fracture geometry, being a function of the occurrence of interconnected permeable structures.C1 ranges between 50% and 120%, indicating that both, poorly connected and short channels (C1 < 1), or connected but diagonal channels (C1 > 1) exist.With increasing pressure gradient, C1 approaches 1.This leads to enhanced continuity in shorter channels, while longer pre-existing channels adopt a straighter path due to improved interconnection.A similar directional dependence of continuity was also demonstrated by Talon et al. (2010) who determined a significantly higher effective permeability perpendicular to shearing, although a two-dimensional analysis using effective apertures leads to more homogeneous results (Marchand et al., 2020).
Figure 8b depicts the continuity indicator C2 as a function of the 3D fracture flow path, the global gradient, and with respect to the flow direction.The averages and standard deviations of the statistical indicator C2 are summed in Table 1.A perpendicular orientation of flow and shearing results in a uniform channel tortuosity between 105% and 110%.An increase in the pressure leads to a decrease in channel tortuosity within the already straight channels.A different situation occurs with a parallel orientation of flow and shearing.An initial high channel tortuosity exceeding 115% indicates a substantially prolonged 3D path length within the fracture void geometry.This elevated path length is caused by the inadequate interconnection of large void spaces and higher local pressure gradients.The path of least resistance is often not the direct route and is highly affected by the geometrical tortuosity of the undulating 3D fracture void geometry.As the pressure gradient increases, the channeling becomes more pronounced and independent of the local fracture geometry.Enhanced channel sweeping results in less tortuous and more direct connections and decreases tortuosity and C2.Similar values (∼110%) for the tortuosity of flow paths (C2) when flow and shearing are oriented in parallel were also obtained by M. Wang et al. (2016) in their study on surface roughness for comparable Hurst coefficients (0.8).In addition, increasing roughness and decreasing Hurst coefficient were shown to favor the development of more tortuous main flow paths and an increase in C2 (M.Wang et al., 2016).Crandall et al. (2010) were able to quantify that an increase in roughness results in increased formation of preferential fluid pathways, consequently reducing A1/V1, and causing an increase in channel tortuosity, denoted as C2.
The determined indicators can be compared not only against the pressure gradient but also with each other.A comparison of the geometric indicators V1 and A1 reveals the two dominant regimes of channel and surface sweeping, respectively, in the form of two separated clusters (Figure 9a).Channel sweeping (blue), has an identical A1/V1 ratio whose absolute value depends solely on the underlying fracture geometry.In the case of dominant surface sweeping (red), the A1/V1 ratio is shifted toward the volume and shows significant dispersion effects.As the flow rates increase, the A1/V1 ratio changes toward unity, both clusters tend to converge, and the influence of channel sweeping increases.The combination of both indicators is thus a measure of the proportion between the dominant flow processes within the fracture and the transition from surface to channel sweeping.
The link between the two continuity indicators C1 and C2 is depicted in Figure 9b.In the case of parallel orientation of flow and shearing (red), surface sweeping dominates, and no dependence between C1 and C2 can be identified.In the case of a perpendicular orientation of flow and shearing, C1 dominates but no quantitative statement about the connection of the indicators C1 and C2 can be made.As pressure gradients increase, the ratio for all fracture geometries, independent of flow direction and regime, tends toward unity, representing well-formed and continuous channels without any topographic height changes.

Implications on Large-Scale Reservoirs
The significance of preferential fluid pathway formation extends beyond laboratory-scale applications.In combination with observations from underground research laboratories, such as Grimsel and GeoLaB, direct implications for fractured (geothermal) reservoirs can be derived (Kohl et al., 2023;Schill et al., 2016).Given an open fractured system formed by the displacement of opposing fracture surfaces, the nature of flow can vary between the two extreme cases.Channel sweeping, characterized by highly localized volumetric flow rates, typically occurs when the flow is perpendicular to shearing.On the other hand, surface sweeping predominates when flow aligns parallel to the shearing direction.
Therefore, the extent of channeling in a fracture will also be affected by the geomechanical process both during its genesis and during later tectonic activity.At reservoir conditions, for a typical normal stress setting with vertical shear fractures close to σ 1 -σ 2 and vertical displacement, a horizontal flow between two wells tends to form channel sweeping and improved hydraulic conductivity (Egert et al., 2021;Lang et al., 2018).Having a strike-slip component with horizontal displacement along σ 1 , a surface sweeping with less expectable flow rates and increased nonlinear effects could result.In the geometrical complexity of a fracture system in a reservoir, the co-existence of improved, and impeded flow regimes would result in a strong dependency upon the orientation of the most conductive features and could explain contradicting observations regarding the nonlinear behavior of geothermal systems (Kohl et al., 1997).A dominating flow regime could even change under mechanical stimulations employed to enhance existing fractures (Evans, 2005).Not only the stress-induced shearing affects the channeling processes.The imposed normal stress and resulting deformation of the fracture void space leads to an increased contact area between the two opposing fracture surfaces.Thus, M. Li et al. (2022) showed that an increase in contact area leads to strongly enhanced channeling.Z. Wang et al. (2021), in turn, were able to quantify that a 45% increase in contact area results in a doubling of flow resistance with a single large channel dictating the flow field within the entire fracture.However, flow anisotropy is also evident here.Liu, Wang et al. (2020) show that a parallel orientation of flow and shearing under normal stress conditions leads to single bottlenecks and thus channeling is significantly enhanced, resulting in a reduction in active volume and surface area (V1 and A1).In contrast, a perpendicular orientation of flow and shearing causes only minor changes in channeling.Blöcher et al. (2019) could show experimentally and numerically, that applied normal stresses cause elastic and plastic deformation mainly around asperities, causing an aperture and permeability reduction, enhanced contact ratio, and therefore more localized flow field along the pressure gradient.In addition to stress conditions, other effects such as changes in rock mechanical properties, mineralization, and fracture filling affect the available flow cross-section within the fracture and thus the local flow regime and channeling (Vidal et al., 2015).
The formation of channels affects the transport and exchange of heat with the surrounding host rock, which can have contradictory effects on reservoir performance (Chen & Zhao, 2020).Severe channeling leads to a decrease in V1 and A1 compared to the common parallel plate approach, resulting in highly focused flow and a significantly reduced effective exchange area between the flowing fluid and the matrix as well as reduced heat and solute transfer coefficients.A large proportion of stagnant water zones, which are primarily present in low-aperture regions, asperities, and outside the identified channels, contribute little to the overall heat exchange.In contrast, channeled flow can also increase the thermal reservoir performance by extending the tortuous pathlength C2 of the flowing fluid within the fracture geometry (Guo et al., 2016).As a result, both the heat exchange area and exchange time increase, as a reduced pressure gradient leads to reduced flow velocities.Applied to the two extreme cases of channel and surface sweeping, this implies that a perpendicular orientation of flow and shearing tends to be detrimental to reservoir performance.It tends to give rise to percolated and highly localized channels with low exchange area and channel tortuosity.In contrast, a parallel orientation of flow and shearing tends to form less channeling and a prolonged flow distance, which can be advantageous for heat exchange.In this extreme case, an increase in flow rates and rising nonlinear effects outside the identified channels lead to a decrease in both, the tortuous flow distance as well as the heat exchange area.
Channeling effects are also identified via solute breakthrough curves, as matrix diffusion can be neglected for conservative tracers in typical EGS applications with low permeable reservoir rock (Egert et al., 2020), making the arrival time a direct representation of the subsurface flow regime (Klepikova et al., 2016).In this context, the relationship between channel and surface sweeping is of interest, because localized flow paths in channel sweeping lead to the formation of a dominant peak, while surface sweeping causes strong mixing and long tailing (Moreno et al., 1988).Sorptive tracers are affected by a reduced flow-wetted surface and changing sorption characteristics.

Conclusions
In complex natural systems such as fractured (geothermal) reservoirs, accurate representation of spatially resolved and highly channeled fluid flows is crucial for ensuring sustainable and dependable predictions of reservoir performance.Unfortunately, there are only a limited number of numerical and analytical models available for characterizing the natural flow field and the formation of preferential flow paths within individual fractures without necessitating prior simplifications regarding geometry and/or flow law.As a remedy, we developed a novel workflow to identify individual channels and defined four indicators that can geometrically and statistically quantify the extent of channeling processes.The derived indicators are used to characterize transport-related properties during altered flow regimes and therefore allow predictions of large-scale flow within the fractured reservoir.The two geometric indicators V1 and A1 capture the flow-through volume and heat exchanger area, respectively.The two connectivity indicators C1 and C2 describe the continuity of the main flow path relative to the fracture length and, for the first time, in relation to the unprojected three-dimensional flow path.
Given the four indicators developed, two distinctive flow regimes are identified across all fractures: A focused channel sweeping and a dispersed surface sweeping.The former occurs with a perpendicular orientation of flow and shearing.Short, well-connected, and low-tortuosity fluid pathways lead to severe localization of volumetric flow.Only small portions of the fracture void space actively flowed through, leaving large off-channel areas as quasi-stagnant water zones.In contrast to the widely assumed parallel plate approach, our results reveal that over a third of the total volume flow interacts with only 15% of the fracture volume/area.This flow regime favors earlier thermal and solute breakthroughs.For a flow orientation parallel to the shearing direction, the volume flow becomes less focused, resulting in enhanced mixing.Under laminar conditions, there is a 22% active flow volume and a 19% exchange surface area.With increased flow velocity, nonlinear effects beyond the channels drive a transition toward channel sweeping.The participating volume/area decreased by approximately 2%, accompanied by a reduction of about 3% in tortuous flow path length.Although we can quantify flow and channeling within an individual fracture, extending these findings to fractured reservoirs necessitates additional research.Moreover, numerical findings must be validated through experiments conducted on natural rocks within underground research laboratories.can be reproduced using the properties provided in the methodology and channeling section.Exemplary results of a Navier-Stokes realization (as shown in Figure 4c), with the corresponding MOOSE input file can be accessed through: http://doi.org/10.5281/zenodo.4621333.

Figure 1 .
Figure 1.(a) Example of a generated fracture after shearing and lowering.The surface roughness exponent is H = 0.8 and the upper surface is sheared to the right; (b) Simplified sketch of the grid cells in the xy-and xz-planes and center coordinate/ line; (c) Power spectral density of the surface topography; (d) Boxplot of the effective aperture for the shown fracture.

Figure 2 .
Figure 2. Mean effective aperture distribution of the 30 rough-walled fractures.

Figure 3 .
Figure 3. Results of the mesh sensitivity analysis by uniformly splitting a three-dimensional fracture mesh.Step 0 corresponds to an initial mesh consisting of 16 × 16 × 1 elements, achieving a resolution of eight elements for the third refinement step.

Figure 4 .
Figure 4. (a) Workflow from aperture calculation to channel identification; (b-g) The different steps of the multi-stage filter concept for channel identification and the resulting figures for R chan = 0.35.(b) 3D morphology of the fracture and the local aperture.(c) The calculated volumetric flow rates of the 3D Navier-Stokes solution.(d) Identified and skeletonized channel nodes.(e) Assigned channels with their ID.(f) Channels after background noise removal.(g) The red line indicates the shortest connection along the direction of flow in the identified longest percolated channel.

Figure 5 .
Figure 5. (a) Different tested critical variables for channel criterion with a constant channeling ratio of 25%.(b) Volume flow ratio (R chan ) in the identified channels in proportion to the total volume flow within one fracture.Shown are the identified channels for four different proportions of volumetric flow (R chan = 0.25; 0.33; 0.50; 0.75).All variables shown in this figure are tested with the same realization for flow and shearing oriented parallel and a pressure gradient of 50 Pa.

Figure 6 .
Figure 6.(a) Difference between the localized channel volume V1 and the volume of the non-channeled parallel plate approach V pp as function of the amount of localized flow ratio R chan for all 30 fracture geometries.(b) Difference between the localized channel half-area A1 and the area of the non-channeled parallel plate approach A pp as function of the amount of localized flow ratio R chan .

Figure
Figure 7a depicts indicator V1 as a function of the global pressure gradient and with respect to the flow direction.Mean values and standard deviation are reported in Table1.A perpendicular orientation of flow and shearing results in channel sweeping.Depending on the initial fracture geometry, V1 ranges between 12% and 18%.Increasing flow rates only lead to an increasing total variance, but do

Figure 7 .
Figure 7. Normalized (a) channel volume (V1) and (b) channel area (A1) for flow parallel (red) and perpendicular (blue) to the direction of shearing as a function of the pressure gradient.

Figure 8 .
Figure 8. Normalized (a) channel length (C1) and (b) channel tortuosity (C2) for flow parallel (left) and perpendicular (right) to the direction of shearing as a function of the pressure gradient.

Figure 9 .
Figure 9. (a) V1 versus A1 for all realizations.Channel sweeping tends to show a linear increase of A1 with V1 whereas surface sweeping results in dispersed results.(b) C1 versus C2 for all realizations.Channel sweeping tends to have a uniform continuity, while surface sweeping results in highly scattered data.The red data tends to have higher C1 and less C2 for increasing pressure gradient.

Table 1
Mean and Standard Deviation for the Identified Indicators Concerning the Flow Direction, Minimum and Maximum Pressure Gradient, and Critical Channeling Threshold