Analysis of Fracture Roughness Control on Permeability Using SfM and Fluid Flow Simulations : Implications for Carbonate Reservoir Characterization

Fluid flow through a single fracture is traditionally described by the cubic law, which is derived from the Navier-Stokes equation for the flow of an incompressible fluid between two smooth-parallel plates. Thus, the permeability of a single fracture depends only on the so-called hydraulic aperture which differs from the mechanical aperture (separation between the two fracture wall surfaces). This difference is mainly related to the roughness of the fracture walls, which has been evaluated in previous works by including a friction factor in the permeability equation or directly deriving the hydraulic aperture. However, these methodologies may lack adequate precision to provide valid results. This work presents a complete protocol for fracture surface mapping, roughness evaluation, fracture modeling, fluid flow simulation, and permeability estimation of individual fracture (open or sheared joint/pressure solution seam). The methodology includes laboratory-based high-resolution structure from motion (SfM) photogrammetry of fracture surfaces, power spectral density (PSD) surface evaluation, synthetic fracture modeling, and fluid flow simulation using the Lattice-Boltzmann method. This work evaluates the respective controls on permeability exerted by the fracture displacement (perpendicular and parallel to the fracture walls), surface roughness, and surface pair mismatch. The results may contribute to defining a more accurate equation of hydraulic aperture and permeability of single fractures, which represents a pillar for the modeling and upscaling of the hydraulic properties of a geofluid reservoir.


Introduction
Fractures exert a vital contribution on determining the migration and storage for geofluids, such as groundwater, and hydrocarbons. Thus, the analysis and modeling of fractures are imperative for characterizing reservoirs and simulating their behavior during the production stage. Fluid flow through fractures is traditionally described by the cubic law, derived from the Navier-Stokes equation for the flow of an incompressible fluid between two smooth-parallel plates [1]. Thus, the permeability (intrinsic permeability) of a single fracture may be represented by the equation where e corresponds to the hydraulic aperture. Since fractures are normally rough, the hydraulic aperture differs from the mechanical aperture (separation between the two fracture wall surfaces). The hydraulic aperture is in fact an equivalent value which can be derived from field analysis like tracer tests and hydraulic tests [2] and laboratory fluid flow experiments (e.g., [3]). Several authors have studied the effect of roughness of the walls on fracture permeability working with various materials, such as glass [4], rocks [5], and concrete [6]. These authors included a correction term, a so-called friction factor (f ), on the permeability of rough fractures: where the friction factor is defined by the common base equation where r a is the difference between the highest peak and the lowest valley of the physical roughness, and both terms c 1 and c 2 are constants with slightly different values depending on the author. Thus, the friction factor depends only on the relative hydraulic roughness r a /2e ignoring the frequency (or wavelength) of the asperities. Another widely used methodology derives the hydraulic aperture from the mechanical aperture, E, and the joint roughness coefficient JRC as proposed by Barton et al. [7]: e = E 2 JRC 2 5 , 4 where JRC is derived by comparing the fracture profile obtained with the Barton Comb with the standard tables provided by Barton and Choubey [8]. This methodology is perhaps the simplest and cheapest way to obtain fracture surface roughness values and has been widely used in outcrop studies [9][10][11][12][13]. The disadvantages of this method are related to the moderate resolution (about 1 mm) and the inaccuracy of equation (4) at relatively wide apertures (with respect to the JRC value). For instance, considering a fracture with 100 μm mechanical aperture and a JRC equal to 2, equation (4) gives a hydraulic aperture equal to 1767 μm. Considering the previous arguments, the main objective of this work is to find empirical equations that describe the effect of fracture roughness on permeability at different apertures. In order to reach this goal, some problems which should be overcome are (i) to develop a protocol for mapping the fracture surface, (ii) to evaluate the fracture roughness as a function of the wavelength of the asperities, and (iii) to validate the relationships using a significant number of samples, roughness values, and aperture scenarios.
Various approaches have been reported in the literature for mapping the surface of fractures and faults in the field or laboratory involving the use of Lidar [14,15], laboratory profilometers [16][17][18], and SfM photogrammetry (e.g., [12,19]). Corradetti et al. [19] applied SfM photogrammetry for mapping fracture surfaces obtaining 3D reconstructions with point-cloud densities of equal quality to Lidar-derived data. Among these methods, SfM photogrammetry shows great promise as it is inexpensive (photo-camera and processing software) and extremely flexible regarding the scales and conditions (applicable in the field and laboratory). For example, SfM photogrammetry has been successfully used as an analytical tool to gather geologic data from outcrop studies [20][21][22][23], as well as at smaller scales on fault surfaces [19] and fossilized human footprints [24].
Evaluation of fracture roughness is achieved by implementing the power spectral density (PSD), which provides a more objective description based on the frequency distribution of the asperities in the Fourier domain. This approach has been successfully applied by previous authors for describing the roughness of fractures (e.g., [25]) and fault surfaces [14,15,18,19]. To increase the statistical significance of the results, approximately 2000 fractures were modeled using the software SYNFRAC [25][26][27] creating computer-generated synthetic fractures, using input parameters of the fractal dimension, fracture roughness (an output of the PSD analysis of real fractures), and the standard deviation of the asperities height.
A key benefit of incorporating computer-generated synthetic fractures is the capability to work with a large amount of fracture data to perform direct fluid flow simulations, such as (i) the finite difference method (e.g., [28]), (ii) the finite element method (e.g., [29,30]), (iii) the finite volume method (e.g., [31]), and the lattice-Boltzmann method (e.g., [32]). The lattice-Boltzmann method (LBM) describes the flow of many particles interacting with a medium and among themselves following the Navier-Stokes equation at the macroscopic scale (e.g., [33]). The LBM has been implemented to compute permeability using 3D images of rocks and soft sediments obtained by micro-CT imaging techniques [34][35][36][37][38] and from reconstructed models [39][40][41] generating results consistent with laboratory measurements [40,42]. The simplest LBM is based on the Bhatnagar-Gross-Krook (BGK) collision operator, which consists of a single relaxation time approximation [43]. Despite its widespread use, the BGK-LBM brings some issues, for example, the computed permeability values may be viscosity-dependent [44]. An alternative approach involves the use of multiple relaxation times (MRT) methods, which solve the drawbacks of the BGK method and are characterized by more stability [45][46][47].
In this study, the permeability values of single isolated fractures (synthetic and natural) were calculated via LBM, using the PALABOS open source library [48]. The method and the code have been previously implemented by Zambrano et al. [38] for quantifying the permeability in deformed porous carbonates using X-ray microtomography synchrotron-based images [49]. Following these authors, rather than the BGK method (e.g., [34]), the MRT approach has been adopted in the present study to assure that permeability values are viscosity-independent.
The selected study area, the Roman Valley Quarry (Figure 1(a)), is an inactive quarry located at the northern termination of the Majella Mountain (Abruzzo region, Italy). The Majella Mountain is the orogenic expression of a thrust- 2 Geofluids related anticline, with internal deformation characterized by high-angle normal, strike-slip and oblique-slip normal faults, small folds, multiple sets of opening-mode fractures, pressure solution seams, and deformation bands [50][51][52][53][54][55][56][57]. The Roman Valley quarry has been heavily studied by previous authors focusing on the structural, sedimentological, and diagenetic properties [54,56], and the fluid flow behavior of the fractured carbonates at the macroscale [13].
Here, the bitumen distribution suggests that the main hydrocarbon flow occurred through the damage zones of the principle NW-SE-oriented oblique slip faults [54] ( Figure 1(b)). The distribution of major lithofacies at the Roman Valley Quarry is another element influencing the presence of bitumen which has been previously described by Rustichelli et al. [56] (see Table 1). Concerning the fractures, the most pervasive ones are represented by both   Table 1 for details). (b) Stratigraphic and structural scheme of the study area; oil distribution is arbitrarily representative considering field observations (see Volatili et al. [13] for more details).
Considering the significance of these fractures, this work focused on investigating both cases of open mode and sheared fractures with a small sliding/tearing mode displacement, in the order of millimeters, allowing the assumtion of a negligible wall wearing. For this last case, the mismatch between opposite walls was also computed due to its importance as a mechanism for maintaining fracture openings even at reservoir depths.

Methodology
In this work, we present a multiphase integrated methodology for characterizing fracture surfaces and their effect on permeability. This approach combines fracture surface scanning using structure from motion photogrammetry, a statistical and spectral description of individual natural fracture surfaces, modeling of synthetic fractures, and computational calculation of permeability by fluid flow simulation.

Sample Collection.
During the summer of 2018, a suite of oriented hand samples was collected from the study site comprising three (i.e., Au, B, and C) of the four major lithofacies of the Bolognano Formation present in the quarry ( Figure 2, Table 1). The field sampling procedure involved manually removing blocks containing fracture surfaces that showed minimal signs of physical and chemical weathering. Sampling targeted these specific lithofacies based on accessibility, quality of well-developed fracture surfaces, and the fact that they have been well documented in previous studies focused on fracture distribution and mechanical properties [13,54,56]. After removal from the outcrop but prior to analysis, the rock samples were cleaned using a soft-bristled brush to remove debris and other obstructions but without abrading the surface.

2.2.
Mapping Surface Topography. The workflow for mapping surface topography involves the following key stages: (1) fracture surface image set acquisition and (2) image alignment and three-dimensional digital rock model creation using SfM.
2.2.1. Image Set Acquisition. SfM photo scanning was performed at the University of Camerino photogrammetry laboratory (Figure 3(a)). We used a tripod-mounted Canon EOS 100D with the standard kit lens fixed at 55 mm. Images of samples were recorded inside a photo-lightbox to maintain full control of camera positions and ambient lighting and to reduce shadows and glare ensuring high image quality. To achieve the desired >70% inter-photo overlap, fracture samples were placed on a 360-degree rotating stage and manually rotated at 10-degree angular increments between photos. After completing a full orbit of the object, the camera was reset at a new vertical position, and the next orbit was conducted adhering to the same 10-degree increments but offset one position from the initial starting point. This procedure was repeated along three horizontal rotations from different vertical positions followed by 2-3 photos directly above the object oriented normal to the fracture surface.

Image Alignment and Three-Dimensional Model
Creation. Fracture surface models were aligned and processed using Agisoft PhotoScan Pro (http://www.agisoft.com). For each fracture surface, approximately 63 photos were used as input images to create the digital point cloud model. We follow the procedure described by Carrivick et al. [58] and Zimmer et al. [24] for camera settings and photo procedure, and Pitts et al. [22] for image alignment and point cloud generation.
Agisoft-generated coded targets were placed inside the scene to aid in the imagery processing; these coded targets are automatically recognized by the software and help build connection points between the image sets ( Figure 3(b)). Additionally, a 5-centimeter scale was placed on the sample surface to calibrate the spatial reference.
As a measure to define the error of the model, we follow the methodology established by Corradetti et al. [19]. This calls for modeling a piece of graph paper under the same condition as fracture imaging (same light, relative distance, and number of photos). Under the assumption that the graph paper is perfectly planar, the standard deviation of the height of the scanned asperities is considered as the error of the model [19]. In our case, this value is approximately 20 μm, whereas the point density (points/area) is 34 points/mm. Table 1: Characteristics of lithofacies exposed in the Roman Valley Quarry.

Lithofacies
Thickness Φ m (%) k m (mD) Bitumen distribution Au: Alternation of medium-to coarse-grained bioclastic grainstones (Au1) and medium-grained bioclastic grainstones (Au2 Alternations of two echinoid plates and spines rich facies: fine-grained bioclastic grainstones (C1) and fineto very fine-grained bioclastic packstones (C2). Argillaceous to marly beds (<3 cm thick) are common. Model. Trimmed 3D point clouds were exported from Agisoft as ".xyz" text files. Then, a rectangular subregion of each fracture surface of interest was extracted from the point cloud and processed to remove undesirable trend and eventual noise ( Figure 4). This technique, previously documented by Corradetti et al. [19], consists of (i) removing the artificial trend of the surface, (ii) surface interpolation, and (ii) sampling in a regular grid.

Fracture Roughness Assessment.
A complete description of the fracture roughness is given by the specification of two functions: the probability density function (depending on the media and standard deviation) for heights and the PSD [59] ( Figure 5).
The Fourier power spectrum, P k , is defined as the square of the modulus of the Fourier transform [60]. Considering a cross section of the rock, this profile can be represented as a summation of sinusoidal components, each with its own amplitude, wavelength, and phase. The squared amplitude of each sinusoid component is referred to as its "power" (Figure 5(b)). The power spectrum regulated in an appropriate manner is referred to as the PSD, G k , and it provides a valuable definition of surface roughness. The PSD as a function of k in a bi-logarithmic scale graph of a self-affine function exhibits an apparent linear slope, which is defined from the following power law equation: where the exponent of the power law α is related with the fractal dimension, D [59], as follows: From a physical aspect, the fractal dimension shows the proportion of high-frequency to low-frequency sinusoid components (roughness). High D values are related to greater surface roughness. By stacking and normalizing the power spectra, it is possible to reduce the noise associated with a single profile and create a single spectrum which represents the entire rough surface in each direction (e.g., [15,19]). The MATLAB script used to perform the procedures described above is partially modified from Corradetti et al. [19]. Since a limited number of natural fracture surfaces were available, additional synthetic fracture surfaces were used to strengthen the statistical significance of the results. Following the procedure described by Ogilvie et al. [25], more than 2000 computer-generated synthetic fractures were created using the software SYNFRAC [25][26][27]. SYNFRAC is based on a mathematical model of a rough surface reported by Brown [59]. The software can model open fractures by introducing mismatch values with the spatial and spectral roughness parameters. For the scope of the present study, the mismatch values were not considered for surface modeling, but were measured after the fracture modeling (see Section 2.5.1). The individual fracture surfaces (natural and synthetic) were used to model dilation associated with ( Figure 6): (i) opening mode displacement (joint and/or opened pressure solution seam), various ranges of aperture, and (ii) sliding/tearing mode displacement (sheared joint and/or sheared pressure solution seam). In both cases, it is assumed that both walls are identical. In the second case, because the shear process is minimum, and the displacement is in the order of mm, it is assumed that no physical wearing of the fracture surfaces has occurred. To illustrate different scenarios, a wide range of displacements (opening and sliding/tearing mode) were considered (the PYTHON code for fracture modeling is available at https://github.com/ superrostom/synthetic-fracture).

Lattice-Boltzmann Method and Permeability Computation.
Lattice-Boltzmann simulations were performed using the open-source computational fluid dynamics software PALABOS [48] following the methodology described by Zambrano et al. [38].
This procedure consists of imposing a single-phase fluid flow through a 3D porous media maintaining a fixed pressure gradient between the inlet and outlet opposing faces of the model; the rest of the faces were padded (Figure 7). A bounce-back boundary condition was assigned to the fracture surfaces. An MRT collisional operator [45,46], with a D3Q19 lattice, is used instead of the popular BGK [43] as in Degruyter et al. [34]. Moreover, the geometry of the media is obtained by the SfM photogrammetry outputs and modeling differently than Degruyter et al. [34] and Zambrano et al. [38] who used X-ray micro CT images.
The simulation ended once the imposed steady-state condition was reached (standard deviation of the average 10 where δP/δx is the pressure gradient, μ is the fluid kinematic viscosity, and U is the average fluid velocity per unit of area. The permeability was calculated, using the same procedure, in two orthogonal directions: along strike and along dip (k x and k y , respectively). In the case of sliding/tearing mode displacement, the x-direction corresponds to the slip direction. All the variables are handled in lattice units prior to permeability calculation. For convenience, permeability values were converted to millidarcy which is the most commonly used permeability unit in the oil industry. All the provided values of permeability correspond to a volume of 1 25 × 10 −4 m 3 (50 × 50 × 50 mm 3 ).

Mismatch Evaluation.
The mismatch between the opposite fracture walls is of extreme importance since this factor may keep fractures open even at reservoir depths. Since the mismatch was not imposed during fracture modeling, it was measured after the generation of the synthetic fractures. The mismatch was evaluated only for the sliding/tearing mode fractures, whereas it was unnecessary in the case of opening mode fractures since the aperture is constant. For the evaluation of the mismatch value, the methodology of the power spectral density ratio (PSDr), introduced by Ogilvie et al. [25], was followed. The methodology consists of obtaining a relationship between the PSD of the aperture and both surfaces of the fracture, as follows: PSDr = PSD aperture PSD upper wall + PSD lower wall 8 The results of this calculation can be represented in a graph where the parameters associated with the mismatch and the degree of mismatch between the surfaces at different wavelengths can be obtained (Figure 8).

Geofluids
Following the definition of Ogilvie et al. [25], these parameters are the following: The calculation of these parameters was made using a MATLAB code. In this case, (ML_min) is considered as the only reliable indicator of the mismatch since (ML_max) often falls outside the scale of the study (Figure 8).

Results
The results of this work consist of an analysis of surface topography performed on fracture samples from three lithofacies (Au, B, and C), and the computed permeability in function of the fracture properties, including fractal dimension, opening and sliding/tearing displacement, and minimum mismatch length. Table 2, the values correspond to the fractal dimension (D), the average height of asperities, and their standard deviation. The traditional roughness measurement, JRC, was added to compare both techniques and the results with previous works in the same outcrop. For the fracture description, we followed Agosta et al. [54]. The set PS2a corresponds to pressure solution seams (generated during background deformation) often sheared, with normal or left lateral kinematics, often impregnated by tar. The set PS2b corresponds to pressure solution seams (generated during background deformation) which are generally un-sheared.

Permeability Results.
The results of the work are presented in Figures 9-14. In each graph, the different surface roughness values (expressed in fractal dimension) are illustrated. When pertinent, error bars are added to show the variability of the results.
In the case of opening mode displacement, the results indicate that the permeability increased proportionally to the mechanical aperture following a positive power-law relationship (Figure 9(a)). Similarly, the hydraulic aperture (derived from equation (1)) is related with the mechanical aperture by means of a positive power law with exponents varying from 1.6 to 1.8 for fractal dimension (D) values between 1.6 and 2.5, respectively (Figure 9(b)).

Geofluids
With respect to the case of sliding/tearing displacement, the results indicate that the permeability component parallel to the displacement (k x ) increases proportionally to the sliding/tearing displacement (S x ) following a positive power-law relationship (Figure 10(a)). The permeability component perpendicular to the displacement (k y ) is also related by a power law with the sliding/tearing displacement (S x ). The anisotropy ratio, k y /k x , is generally higher for low values of fractal dimension (Figure 10(b)). The anisotropy ratio value tends to decrease as a function of the sliding/ tearing displacement (S x ) following a negative power-law relationship. The highest value of anisotropy was near 2.6 for fractal dimension (D) equal to 2.0 and sliding/tearing displacement (S x ) equal to 0.5 mm.
The fracture roughness, expressed in terms of fractal dimension (D), showed a different influence in the permeability on single fractures depending on their kinematic: opening or sliding/tearing mode. For the opening mode displacement case, the permeability is inversely proportional to the fractal dimension (Figure 11(a)). This relationship   the fractal dimension following a positive power law (Figure 11(b)). In this case, an increment of the fractal dimension (roughness) from 2.0 to 2.5 represents an enhancement of the permeability of approximately one order of magnitude. It is expected that this positive relationship between displacement and permeability should stabilize at a certain 11 Geofluids point, as the permeability and porosity cannot increase indefinitely. This behavior is observed when the porosity is evaluated at higher displacement values ( Figure 11). In this case, thanks to the simplified calculation of porosity in comparison with permeability, a large volume of data was considered (more than 2000 fractures). The porosity is proportional to the sliding/tearing displacement following a nonlinear relationship with variable slope. The most significant change of slope occurs near 10 mm and after approximately 30 mm of sliding/tearing displacement, where the porosity seems to stabilize at values between 3% and 4%. The fractures with higher fractal dimension (roughness) tend to have higher porosity, which agrees with the permeability results. As previously mentioned, these results do not consider the possible 12 Geofluids physical wearing of the surfaces due to shearing particularly at high sliding/tearing displacement. The fracture permeability is related to the porosity following a power law (Figure 13), which differs from the theoretical relationship based on the smooth parallel-plate equation; here, we also assumed smooth parallel plates for the porosity. The power-law relationship seems to be unaffected by the roughness (fractal dimension). Similar behaviors were obtained for both permeability components k x and k y .
The minimum mismatch length, ML_min, was evaluated as a function of the displacement (Figure 14(a)). Results indicate that within the evaluated range (a maximum displacement of 5 mm), the minimum mismatch is linearly proportional to the displacement. The permeability is proportional to the minimum mismatch length following a power-law relationship (depending on the fractal dimension). In other words, higher values of fractal dimension imply a higher permeability for similar mismatch values.

Discussion
The present work evaluates the effect of fracture surface features such as roughness, aperture, and mismatch on permeability using fracture surface scanning by SfM photogrammetry, numerical modeling, and lattice-Boltzmann fluid flow simulation.

SfM Photogrammetry Surface
Scanning. The results of this study demonstrate the versatility of the SfM procedure as an analytical tool which can be applied at a wide range of scales including millimeter-scale features such as fracture surfaces. The controlled conditions in the photogrammetry laboratory allowed a highly detailed scan and extraction of the micro surface topography of samples sized 30 × 30 × 30 cm producing a point cloud with a density of 34 points/mm and an estimated error of 20 μm. This method produces more realistic and applicable results than the traditional Barton Comb, with results comparable to those reported by Candela et al. [14], Renard et al. [18], and Corradetti et al. [19] using both Lidar or laser profilometers. However, the SfM methodology is several orders of magnitude more cost-effective and is readily accessible. A future implementation of this study could include developing a workflow suitable for in situ field studies. However, more variables need to be controlled and results yielding lower accuracy are expected.

Fracture Roughness
Characterization. This methodology proved to be highly efficient in expressing the fracture roughness allowing a more accurate and representative measure with respect to the relative hydraulic roughness [4][5][6] or the JRC [7,8]. However, more data is necessary for evaluating control of the lithofacies in the fracture roughness, as it is observed for other properties such as distribution and spacing (e.g., [61]).
Another important aspect of these results is the reproducibility of synthetic fractures with similar characteristics. This step was important to increase the data volumes leading to a greater statistical significance of the results and validity of the inferred relationships. The lattice-Boltzmann procedure also played a key role in this study as it allows the estimation of permeability values for controlled scenarios with different imposed properties (i.e., roughness, opening mode, and sliding/tearing displacement). This permits evaluation of the relationship between permeability, porosity, mismatch, and other imposed properties. The computed permeability may present some inaccuracy in low-resolution models as previously reported by Zambrano et al. [38]. Therefore, permeability results may be considered not as the real values but only as approximations.

Permeability in Function of Fracture Properties.
Two situations were considered to explain the presence of open fractures: (i) dilation due to opening mode displacement (joint or opened pressure solution seam) and (ii) dilation due to mismatch caused by shearing and sliding/tearing displacement.
In the first case, the results followed the expectation and confirmed previous interpretations: (i) permeability tends to increase with opening following a nonlinear relationship, (ii) a higher fractal dimension (greater roughness) correlates to lower permeability, and (iii) the effect of roughness is less significant at greater opening values. It is expected that higher roughness (higher frequencies of asperities) may expose a wider area in contact with the migrating fluid, diminishing its velocity due to friction. Evidently, at higher opening values, this effect should be less evident because it 14 Geofluids is the specific area (area/volume) which has a control on the permeability, as has been previously reported by Zambrano et al. [38] for porous media. The second case (sliding/tearing displacement) creates an aperture due to the mismatch between the opposite walls of the fracture. We found significant differences between these two cases concerning the effect of the roughness of the permeability. In fact, the effect of roughness on permeability is inverse. Given the same displacement, fractures with higher roughness values permit the creation of larger voids and therefore enhance the fluid flow. So, the effect of friction exerted by the roughness on fluid flow has a secondary role in the case of mismatch due to sliding/tearing displacement. The continuous dilatancy of the fracture due to sliding/ tearing displacement should cease at a certain value depending on the asperity frequencies present in the fracture. Nevertheless, it is difficult to verify this behavior for real fractures, where a 20 mm sliding/tearing displacement likely leads to  fracture wall wearing and the generation of cataclastic material [54], eventually reducing permeability. However, our model has greater applicability to small displacements where the damage of the fracture walls is negligible. The permeability anisotropy in fractures with sliding/ tearing displacement is significant and dependent on the roughness (fractal dimension). For low displacement (0.5 mm), the anisotropy can reach values up to 2.6 for fractures with a fractal dimension of 2.0. For the same displacement, fractures with high roughness (D = 2 5) showed lower values of permeability anisotropy near 1.5. For higher sliding/tearing displacement, the permeability anisotropy decreases approaching the value 1.
The mismatch itself has a positive control on the permeability. The importance of this result is that the mismatch could also be produced by diagenetic processes (e.g., cementation, dissolution) and shearing. Zambrano et al. [38], using X-ray microtomographic images in shear compactive bands hosted in porous carbonates, showed the existence of complex channelized porous networks along the shear surface within these structures. Permeability is directly proportional to the minimum mismatch length, which is related to the shear displacement. 16 Geofluids values of aperture and the most important bitumen impregnation [13,54,56]. The relationship between permeability and porosity for rough fractures clearly deviates from the ideal smooth parallel plate case (for the studied scenarios). Fracture permeability is lower for the same porosity range (<0.2%) in comparison to the theoretical values. Instead, the powerlaw slope is higher, indicating a more important control of porosity as it was expected. The equation itself may be useful to estimate the permeability of fractures if the fracture porosity is known.
After their formation, both closing and opening mode fractures are often subjected to a shear process, and even with a small imperceptible sliding/tearing displacement, they cannot be modeled as simple opening mode fracture. At reservoir depth, preexisting fractures (joints and pressure solution seams) that are favorably oriented to be sheared (accordingly to the orientations of the stress field which affected the area) may be characterized by a mismatch between the fracture walls enhancing the fracture opening. Therefore, the findings of this work may have a significant impact on fracture modeling workflows for subresolution faults (e.g., [62]). In fractures with small shear displacement, contrary to the open mode case, the roughness may influence positively the permeability. Sliding/tearing displacement also may enhance the fracture permeability eventually decelerating when the cataclasis process starts. If the shear displacement is small (S x < 1 mm), the permeability anisotropy is significant enough to be considered in the fracture modeling workflow.

Conclusions
We presented a new multifaceted approach to characterize surface fracture roughness by SfM photogrammetry, numerical modeling, and computational fluid dynamics simulation. This methodology provides a better quantification of surface parameters that are not possible to obtain using former surface roughness measurement and analysis tools.
In addition, this study illustrates the crucial relationships between permeability and other fracture properties, such as roughness, porosity, opening mode-sliding/tearing displacement, and mismatch. The obtained relationships pointed out the following statements: (i) In joints (opening mode fractures) and/or opened pressure solution seams, the roughness tends to reduce the permeability. Thus, the permeability is inversely proportional to the fractal dimension (ii) In sheared joints and/or pressure solution seams (assuming an insignificant surface wearing), the sliding/tearing mode displacement may cause mismatch and therefore enhance the porosity and permeability. The validity of this behavior may depend on the point that displacement starts to produce cataclastic material. Small shear displacements and mismatch may be extremely important to guarantee storage and migration of geofluid at depth thanks to asperities-supported aperture. Permeability anisotropy is very significant for small shear displacements, characterized by higher values of permeability component perpendicular to the shear displacement (iii) Porosity exerts a more important control on permeability in rough fractures (higher power-law slope).
The empiric relationship may result in greater utility for estimating the fracture permeability if the fracture porosity is known

Data Availability
The data used to support the findings of this study are available from the corresponding author upon request.

Conflicts of Interest
The authors declare that they have no conflicts of interest.