Effect of Shear Displacement on the Directivity of Permeability in 3 D Self-Affine Fractal Fractures

State Key Laboratory for Geomechanics and Deep Underground Engineering, China University of Mining and Technology, Xuzhou 221116, China State Key Laboratory of Mining Disaster Prevention and Control Co-Founded by Shandong Province and the Ministry of Science and Technology, Shandong University of Science and Technology, Qingdao 266510, China School of Engineering, Nagasaki University, 1-14 Bunkyo-machi, Nagasaki 852-8521, Japan Center of Rock Mechanics and Geohazards, Shaoxing University, Shaoxing 312000, China Hubei Subsurface Multi-Scale Imaging Key Laboratory, Institute of Geophysics and Geomatics, China University of Geosciences, Wuhan 430074, China


Introduction
Understanding the hydraulic response of fractured rock masses with respect to shear stress is of great importance in various fields of rock engineering, such as the construction and maintenance of dam foundations and radioactive waste repositories, unconventional oil and gas exploration, and hydrocarbon recovery [1][2][3][4].In tight rock masses, the rock matrix has much less conductivity when compared to fractures, which means that the fracture network forms the main flow channels [5][6][7].In a natural rock mass, the fracture network consists of numerous single fractures [8][9][10][11].Therefore, in order to improve the understanding of the hydraulic properties of fractured rock masses, the hydraulic response of a single fracture, which is the fundamental element of a fracture network, should be initially studied [12][13][14].
Many previous studies have focused on the hydromechanical properties of fractures through theoretical modeling, laboratory experiments, and numerical simulations.Zimmerman et al. [15] compared the flow through rough fractures with parallel plates to examine the tortuosity induced by the geometry of void spaces.They found that in a rough fracture, fluid flows through connected channels that bypass the closed zones.This phenomenon is the channeling effect that has been observed in other studies [16][17][18].Rasouli and Hosseinian [19] investigated the effect of profile roughness on the hydro-mechanical behavior of fractures and proposed correlations to relate the hydraulic aperture with the geometrical properties.The results suggest that the fracture surface roughness is one of the most important factors that influence the permeability of a fracture.
The permeability of a fracture is dominated not only by the geometrical characteristics of the fracture surface but also by the states and evolutions of the stress fields [20].Li et al. [21] conducted a series of laboratory tests to examine the shear effect on the evolution of the permeability of natural rough rock fractures.The test results showed that the aperture distribution and the permeability of the fractures were altered according to the shear displacement and normal constraint conditions.In previous studies, the flow behavior has been mostly analyzed only for flow directions parallel and/or perpendicular to the shear direction [22][23][24][25].Whereas, in nature, the fracture orientation and position are arbitrarily distributed with respect to the in situ stresses, and the flow direction, which is controlled by the hydrogeological conditions, is not always exactly parallel or perpendicular to the shear direction [26,27].As indicated in Figure 1, the flow direction is consistent with the macroscopic pressure drop direction, and the shear direction varies accordingly with the applied principal stresses.The different inclinations between the flow direction and shear direction lead to a diversity of flow patterns in rough fractures.
In this study, to characterize the hydraulic properties of natural rough fractures, a self-affine fractal model is proposed to generate fracture surfaces that have different joint roughness coefficients (JRCs) using the modified successive random addition (SRA) method.Then, the evolutions of the aperture distribution and fluid flow during shear are simulated based on the topographical data of the generated fracture surfaces.The influence of the angle between the macroscopic flow direction and the shear direction on the flow paths and permeability are estimated.The hydraulic anisotropy that results from the different inclinations between the shear direction and the macroscopic flow direction is systematically investigated.

Fracture Surface Generation and Aperture Distribution
Due to the difficulty associated with the repeatability of fracture specimen analysis and sealing treatments in laboratory tests during shearing, numerical computation is an appropriate method for analyzing the flow behaviors in various flow directions, and the precise modeling of a fracture surface is a critical issue.Previous research has shown that the complex geometry of an object in nature can be modeled by fractal theory [28][29][30][31].The most important property of fractal geometry is self-similarity or self-affinity.In case of self-similarity, the object is scaled by the same amount in all directions to ensure that the enlarged profile looks similar to the original, but in self-affinity scaling, it is not necessary to be identical in all directions.The geometrical features of a rough fracture surface are fractal in nature and typically follow the self-affinity distribution, which can be modeled using fractional Brownian motion (fBm) [32,33].The SRA algorithm, which is recognized as an efficient and fast algorithm, is used to employ fBm to generate the fracture surface.As shown in Figure 2, four self-affine rough fracture surfaces of 256 mm in size with different Hurst exponents (H) are generated using the modified SRA algorithm described in Liu et al. [34].The roughness of the generated fracture surfaces can be represented by a scalar parameter denoted as H that is related to the fractal dimension (D) by D = 3 − H [35].The JRC of each fracture surface is estimated according to the following equation [36]: where Z 2 is the root mean square of the first derivative of the profile that can be expressed in the following discrete form: where x i and z i are the coordinates of the fracture surface profile and N t is the total number of sampling points along the length of fracture.A sampling point interval of 0.5 mm is selected for the JRC calculation.The mean value of JRC for all the selected profiles is calculated to characterize the roughness of the fracture surface.A larger value of H results in a smoother surface with a smaller JRC (Figure 2(a)).For a fracture surface with a small H (Figure 2(d)), the calculated JRC is close to 17, and plenty of precipitous asperities exist within the fracture plane. 2 Geofluids The hydraulic response of a rough fracture is highly sensitive to the fracture aperture, which is defined as the difference between the two opposing walls of a fracture.In this study, the fracture surfaces are assumed to be well matched at the initial state.After implementing a shear displacement, the two fracture walls shift horizontally due to the shear stress and separate vertically due to shear dilation, as shown in Figure 3.For each shear displacement interval, the opposing points on the two fracture walls are changed, and the local aperture b is recalculated according to the following equation: where u denotes the shear displacement, Z x, y represents the aperture asperity height of the fracture surface, and the value of b x, y = 0 indicates the contact point between the two fracture surfaces.u v is the normal displacement generated at u, which can be estimated by where γ is the dilation rate, which can be calculated based on the analytical model proposed by Indraratna et al. [37].
A total shear displacement of 20 mm is applied to each fracture.Because of the gradually decreasing nominal contact area between the upper and lower fracture surfaces during shear, the aperture field of 200 mm × 200 mm is extracted from the central part of the original model to ensure that the analyzed area is constant.For each shear displacement, the aperture fields in several radial directions are extracted with inclination angles of θ = 0 °, 30 °, 60 °, 90 °, 120 °, 150 °, and 180 °to the shear direction, as shown in Figure 3. Since the flow behaviors in the extracted models are symmetrical along the center of the original model, the permeability of the model having an inclination angle = θ equals to that having an inclination angle = (θ + 180 °).Therefore, the inclination angles ranging from 180 °to 360 °are not presented.

Flow Calculation
Fluid flow through a rough fracture is governed by the Navier-Stokes equations.Due to the high nonlinearity of the system, it is usually difficult to obtain the exact solutions of these equations.For a rock fracture consisting of two perfectly smooth parallel plates, the fluid yields the well-known cubic law that relates the volumetric flow Q to the macroscopic pressure gradient ∇P [38]: where μ is the dynamic viscosity.Since natural rock fractures usually have anisotropic apertures formed by two rough surfaces that are in partial contact, it is difficult to assign a unique value to b in (5).When considering the local flux as a conservative quantity, the cubic law yields the well-known Reynolds equation: 3 Geofluids In ( 6), b can be defined point by point to represent the 3D distribution of aperture fields in the rough fracture.This equation has been widely utilized to study fluid flow through rough fractures [39][40][41].This approach is accomplished by simplifying the void space between the two fracture surfaces into a series of connected parallel plates.By assigning spatially variable apertures to these parallel plates, the heterogeneities of the fracture aperture field can be incorporated into the model.Then, the total flux through the fracture is equal to the integral of local fluxes through the cross section of the fracture.The equivalent permeability K of the fracture is evaluated by where L is the fracture size and A is the cross-sectional area.
The finite element method (FEM) software COMSOL Multiphysics is used to simulate the fluid flow through fractures during shearing.The aperture field is divided into 160,000 (400 × 400) small square elements in the FEM model.Shear displacements of 3 mm, 5 mm, 10 mm, 15 mm, and 20 mm are applied to the fracture.At each shear step, the previous aperture field is redistributed, and some new contact points and void spaces are generated.Thus, the aperture distribution should be recalculated for each shear interval.Two fixed hydraulic heads of 0.2 m and 0 m along the two opposite boundaries are applied to the model.The other two boundaries are assumed to be impermeable.The influence of shear displacement on the aperture fields at other inclination angles shows similar tendencies, which therefore are not displayed.The frequency corresponding to b = 0 decreases as u increases, which indicates that the contact area between the two fracture surfaces decreases and more void spaces are generated by shearing.An increase in the JRC also has a similar impact on the contact pattern and reduces the contact areas.For a fracture with JRC = 2 45, the distribution of the aperture varies from 0 to 0.8 mm at u = 0 mm.As u increases, the shape of the frequency curve changes from a sharp curve to a flatter curve, in which the aperture varies from 0 to 2.3 mm at u = 20 mm.This indicates that a more heterogeneous aperture field is generated as u increases.A comparison of the frequency curves of the four fractures at the same u shows that as the JRC increases, more significant flattening is observed and accompanied by gradual increases in the mean aperture and deviation.All aperture fields are extracted from the models parallel to the shear direction, with θ = 0 °.At u = 5 mm, a number of contact areas exist, with a small void space distributed sporadically in the fracture.Some flow paths are generated, but the channeling flow effect is not obvious.As u increases, the contact areas decrease dramatically, which further confirms the decrease of the frequency corresponding to b = 0 in Figure 4.At the same time, the void spaces gradually concentrate in a few major channels in which the aperture is relatively large.A significant channeling flow is observed with the flow clustering in the higher aperture channels and bypassing the lower aperture regions and contact areas.This phenomenon is consistent with the results obtained by many previous studies [42,43].

Results and Analysis
Figure 6 compares the aperture and flow fields obtained at u = 10 mm for four fractures with different JRCs.All aperture fields are extracted from the models parallel to the shear direction, with θ = 0 °.Since the surface of fracture with JRC = 2 45 is flat and smooth, the generated contact areas are larger, and the void spaces are smaller compared with the widely distributed asperities on the fracture surface for JRC = 17 28.Shear displacement induces different aperture fields for the four fractures and thus results in a diversity of channeling flow patterns, depending on how the voids and barriers are oriented with respect to the macroscopic flow.
Figure 7 indicates the aperture and flow field for fractures with JRC = 17 28 and u = 10 mm when the inclination (θ) between the shear direction and macroscopic flow direction is 0 °, 30 °, 60 °, 90 °, 120 °, and 150 °, respectively.As the shear direction is fixed, the hydraulic gradient is applied such that the angle between its direction and the shear direction varies from 0 to 150 °, correspondingly.The main flow paths through the fracture vary as θ changes.Comparisons between the aperture distributions at θ = 0 °and θ = 90 °show that the shear process is more prone to generate striped contacting asperities and void space channels in  6 Geofluids the direction perpendicular to the shear displacement.These striped contacting asperities tend to change the flow paths of the model at θ = 0 °, and, at the same time, the void space channels tend to promote the flow of the model when θ = 90 °.When θ is closer to 0 °(i.e., θ = 30 °and 150 °), the flow path becomes more tortuous.However, when θ is closer to 90 °(i.e., θ = 60 °and 120 °), the macroscopic flow becomes more parallel to the void space channels, which results in a higher permeability.4.3.Directivity of Permeability during Shear.Figure 8 shows the calculated equivalent permeability (K) for fractures with different JRC, u, and θ values.For all cases, as u increases from 3 mm to 20 mm, K increases by approximately 1~2 orders of magnitude.As the JRC increases from 2.45 to 17.28, K increases by approximately 1~3 orders of magnitude.The variation of K is consistent with the change in the fracture apertures induced by u and the JRC.
The directional permeability contours that are plotted in a polar coordinate system for fracture surfaces with different JRCs are shown in Figure 9.The curve is symmetrical about its center; thus, only magnitudes of K for models with θ = 0 °-180 °are calculated.The largest K is mainly oriented in the direction perpendicular to the shear direction, with θ = 90 °, and the smallest is oriented in the direction parallel to the shear direction, with θ = 0 °.The ratios (K/K 0 ) relate the directional permeability to the permeability at θ = 0 °(K 0 ) for fracture surfaces with different JRCs.These ratios are calculated and plotted in Figure 10.K/K 0 generally increases as θ approaches 90 °.The largest K/K 0 for each fracture varies from 2 to 3 when θ = 90 °, indicating that the permeability in the direction perpendicular to the shear direction can be two or three times larger than that in the direction parallel to the shear direction.For each fracture, K/K 0 decreases as u increases, indicating that the directivity of the permeability becomes less obvious with increasing shear.Comparisons of the four fractures at the same u show that K/K 0 decreases as JRC increases, which indicates that the directivity of permeability is less variable in the rougher fracture.This result is  The results of this paper can provide some fundamental understanding of the hydraulic properties for the fracture with an arbitrary direction with respect to in situ stress.
However, many investigations have shown that the hydromechanical behavior of rock fractures is also scale dependent, and the hydraulic properties determined by small laboratory samples cannot be directly extended to the fracture at field scale.Therefore, when the results of this paper applied to the engineering field, the scale effect should be considered.In future works, this issue will be systematically studied.

Conclusions
In this study, the fluid flow in 3D self-affine fractal fractures under shear was simulated, and the effects of shear displacement on the anisotropy of the aperture field and the directivity of the permeability were analyzed.The modified successive random addition (SRA) method was proposed to generate a series of rough fracture surfaces with self-affine fractal characteristics.The variation in the aperture distribution during shear was estimated based on the analytical model that describes the complete shear behavior of a rough fracture.The calculated aperture field was imported into COMSOL Multiphysics to simulate the fluid flow using finite element method (FEM) code.The evolutions of the flow patterns and the permeability that result from shear displacement with different inclinations with respect to flow direction were investigated.The results show that as the shear displacement increases, the contact areas decrease and the void spaces increase.The mean aperture and its deviation increase, and the frequency curve of the aperture distribution changes from sharp to flat.This change is mainly due to the shear-induced dilation and roughness of the fracture surface.Under the same shear displacement, a smaller contact area and a larger mean aperture are generated for a fracture surface with a larger JRC.Channeling flow becomes gradually obvious through the 10 Geofluids fracture as shear displacement increases, and a diversity of channel flow patterns exists in fractures with different JRCs.The equivalent permeability of the fracture increases by almost two orders of magnitude as the shear displacement increases from 3 mm to 20 mm.As the JRC increases from 2.45 to 17.28, the equivalent permeability of the fracture increases by almost three orders of magnitude when the shear displacement is 20 mm.The inclinations between the shear direction and macroscopic flow direction have significant effects on the permeability of fractures.The largest equivalent permeability usually exists in the direction perpendicular to the shear direction, and the smallest usually exists in the direction parallel to the shear direction.The equivalent permeability tends to become larger as the macroscopic flow direction approaches the direction perpendicular to the shear.The ratio of directional permeability to the permeability in the direction parallel to the shear direction varies between 1.03 and 2.71.This ratio tends to decrease as the shear displacement and JRC increase, which indicates that the directivity of the permeability is more obvious for fractures with small JRCs under small shear displacement.
In the future, the investigation will be extended to the hydraulic properties of 3D fracture networks with geometries similar to those of natural fractured rock masses.From a general perspective, we will study how different factors, such as the fracture heterogeneity and shear stress, impact the permeability tensors of 3D fracture networks.

1 FFigure 1 :
Figure1: A schematic view of shear direction and flow direction in a fracture[24].

Figure 3 :
Figure 3: Generation of the original aperture distribution with the upper fracture surface displaced from the lower fracture surface and extraction of fracture model with different rotations.

4. 1 .
Aperture Distribution during Shear.Figure4shows the statistics of aperture distribution at different u for fracture surfaces with different JRCs.All aperture fields are extracted from the models parallel to the shear direction with θ = 0 °.