Abstract

This study explores the effects of normal loading and shearing on hydraulic properties in roughness-walled rock fractures. The geometries of five fractures were measured by the 3D scanning technology. The flow simulation was performed for rough rock fractures with large displacements during normal loading and shearing by finite volume method (FVM). The results demonstrate that the deformation of fracture with increasing normal stress and shear causes nonuniform changes in void space geometry and further influences fracture permeability. Associated with normal displacement are an increase in contact area and a decrease in mechanical aperture. The transmissivity is decreasing by 3 orders of magnitude response to applied normal displacement values of 0.0 mm to 1.8 mm. In contrast, an increase in mechanical aperture and contact ratio that occurs with increasing shear displacement values of 0.0 mm to 4.0 mm is associated with decreasing distinctly transmissivity by 1.5–2 orders of magnitude. Based on the numerical results, an empirical equation is proposed to evaluate the effects of contact area and roughness of fracture on the hydraulic aperture. The good agreement between numerical results and the predicted results by the new model indicates that the proposed model is capable of estimating the hydraulic aperture of rock fractures through parametric analyses, compared with other published models from available literature. In addition, the new model succeeds in predicting the transmissivity in Develi and Babadagli (2014) water flow experiments.

1. Introduction

Fractured rocks are complicated geological media that contain numerous fractures of various attitudes and scales [1, 2]. The in-depth exploration of hydraulic characteristics of rock fractures is an important issue for safety assessment of rock engineering, such as hydraulic engineering, migration of contaminant control, geothermal exploitation, and hazardous wastes isolation [36]. In particular, the evolution of permeability of rock mass with complex disturbed stress conditions has been prompted by an interest in identifying the development of excavation damage zones around opening [3, 7].

Redistribution of in situ stresses changes the void space geometry of fractures and its fluid flow behaviors, when being active in fractured rock mass [8, 9]. In situ stresses are divided into normal and shear components. The effect of normal stress and shear on aperture variations and asperity geometry is different for fractures. Gentier et al. [9] found that the value of closure was not uniform in the fracture and it increased with increasing normal stress (1 to 25 MPa). The nonlinearly decreasing relations of mean aperture and normal stress were obtained by Kulatilake et al. [10]. It was showed that the aperture fractal dimension decreased rapidly at low normal stresses and kept constant value with increasing stress. Compared with normal stress effects, there is more complex variation in geometry of fractures with shearing. Yeo [11] analyzed the aperture changed during shear and showed both mean aperture and its standard deviation increased, with increasing shear displacements. The same results were reported by Koyama et al. [12] and Matsuki et al. [13]. In addition, Watanabe et al. [14] found that, due to damage of asperities in fracture, the flattener fracture surface was generated during shear, and this would provoke the decreasing surface roughness.

For steady laminar flow in smooth fracture, the cubic law is commonly used to describe the permeability which increases in proportion to the square of the mechanical aperture width. However, the natural rock fracture is characterized as complex geometrical characteristics, such as rough surface, contact areas, and uniform aperture distribution, leading to flow rate being lower than predicted by cubic law. A few modified cubic laws [4, 13, 15] were therefore proposed to quantify these effects on flow behaviors. To further take into account the geometrical characteristics of fracture, fluid flow behaviors were mainly studied by means of the laboratory test and numerical simulation.

The laboratory flow and tracer tests were conducted in Gentier et al. [9], considering the effect of normal stress. The results showed that the transmissivity of fractures decreases with increasing normal stress. The effect of contact areas on flow channeling had been fully understood in flow visualization experiments performed by Develi and Babadagli [16] and Zhang et al. [17]. Chen et al. (2015) reported the flow tests results through fracture under confining pressure (5–40 MPa). One notable difficulty in testing fluid flows through rock fractures with shear was sealing of fluid. A new shear-flow testing apparatus with specially designed fluid sealing techniques for rock fractures was developed by Koyama et al. [18], under constant normal load (CNL) or constant normal stiffness (CNS) constraint. Yeo et al. [19], Watanabe et al. [14], and Xiong et al. [15] also carried out fluid flow experiment in fractures with shear displacements. These researches make a significant progress in deepening understanding of the hydromechanical and transport behavior of rough-walled rock fractures.

Although the laboratory tests are more perceptive and actual presentation of flow behavior, there are some obvious limitations, such as the difficulties of quantitative measurements of changing surface roughness and observations of fluid behaviors inside fractures [18]. Numerical methods therefore provide a good treatment. The tortuous flow fields and channeling effects in fracture were investigated by Matsuki et al. [13], Xiong et al. [15], and Koyama et al. [18], using numerical simulation. Zou et al. [20] developed a 2D finite volume method (FVM) and observed the dynamic evolution of eddy and back flow in fractures. The perpendicular velocity profile at the cross section of fracture did not obey parabolic function, when increasing flow rate. The obstacles effect of contacts for flow was evaluated by Yeo [11], with help of 3D flow simulation methods. Likewise, Auradou et al. [21] found the phenomenon of flow anisotropy induced by shear with different directions.

Traditionally, the Reynolds equations and Stokes equations, which serve as governing equation, have been used to model flow behaviors through rough-walled rock fractures [3, 22] (Thompson and Brown 1991). The Stokes equations and Reynolds equations, obtained by neglecting the inertia terms in Navier-Stokes equations, were verified to work in fractures at low flow rates [19, 23]. However, the experimental transmissivities were about 10–60% less than those predicted by the Stokes equations and Reynolds equation, as found in [23, 24]. The above-mentioned investigation results show the necessity of using the Navier-Stokes equations to model flow through rock fractures.

The influences of fracture geometry and stress on fluid flow through rock fractures have been well studied; however, there is still lack of mathematical models to fully describe their interactions. Numerical simulations by solving the 3D Navier-Stokes equations were employed to investigate the fluid flow behavior through the fractures with complex geometries in this study. A series of flow rates were imposed on the five groups of fracture with normal and shear deformations. The response of transmissivity and hydraulic aperture was thus evaluated. The evolution of geometrical characteristics of fractures during closure and shear and its effects on the fluid flow behavior were discussed.

2. Measurement of Fracture Specimen

2.1. Fracture Specimens Preparation

The granite is widely distributed in the continental crust. Granite fractures in the bam areas control the safety of the bam and hydraulic structure. The granite is originated from quarry of Dagangshan hydropower project in Dadu Reviver, China (see in Figure 1). The granite was biotite monzonitic granite in Jinning-Chengjiang period and experienced main squeezing and luffing tectonic and hydrothermal alteration activity. Several blocks with dimensions of 150 mm × 150 mm × 150 mm were obtained by cutting the granite body. The blocks were then split by Brazilian tensile test (Figure 2(a)) into two halves. As shown in Figure 2(a), the granite block was placed between two steel needles. Compressive loads from the top platen were successfully converted to tensile loads by means of these needles. Increasing loads resulted in crack development and propagation from one needle to the other. Once the indirect tensile strength of the rock was exceeded, the block was fractured along a plane and separated into two halves. According to surface texture and topography, five represented fractures were selected and marked as the Gr1, Gr2, Gr3, Gr4, and Gr5, shown in Figure 2(b).

2.2. Measurement of Fracture Geometry

In order to obtain the fracture roughness, the 15 cm × 15 cm area of rock fractures was scanned by an advanced stereotopometric scanning system with an accuracy of ±25 μm introduced by Chen et al. [25]. 151 × 151 data points representing the surface topography were obtained in the direction of flow. The resolution to represent the roughness was illustrated by Develi et al. [26] and found that a 1 mm grid size was capable of capturing the impact of surface roughness on the hydraulic behavior.

It is necessary to determine the initial void space geometry of the fractures. Tatone and Grasselli [27] provided a novel noncontact measurement technique based on the scanning system. Hence, the void space geometry of closed fractures was obtained using the technique. The parameters of fractures are listed in Table 1.

2.3. Fractal Dimension of Fracture Surface

After quantifying the fracture geometry, the fracture surface parameters were measured using the mathematical procedures. The variogram method for determining fractal dimension was employed in 2D profiles extracted from the surfaces. The variogram function is given aswhere is the semivariogram, and are height of the 2D profile from baseline, and is the number of pairs of at a lag distance between them. e can simplify as a power law function in self-affine profile, as :where is a proportionality constant and is the Hurst exponent which is related to the fractal dimension by . For each fracture surface, 31 sectional profiles in the lower part of the fracture surface, parallel to the flow direction at intervals of 5 mm, were extracted to calculate the fractal dimension . The values of are given in Table 1. As seen in Table 1, the lowest value of for the specimen Gr2 is 1.520, while the highest is 1.424 for the specimen Gr1.

Selection of the fractal dimension to be correlated to the fracture hydraulic characteristics is on the basis of previous study. The natural fracture surface is characterized by self-affine fractal profiles. Consequently the fracture roughness can be quantified by fractal dimension [10]. Brown [3] pointed out that the fractal dimension had an impact on fluid flow. Babadagli et al. [28] demonstrated that the exhibited a correlation to the fracture transmissivity in the water flow experiments.

3. Numerical Methodology

3.1. Mechanical Apertures during Normal Loading and Shear

The detailed distribution of aperture inside the fracture cannot be directly measured during normal loading or shear tests but can be accurately simulated by numerical simulations if the geometry data of the fractures are available. In addition, the complex flow can be captured by the simulation. The mechanical aperture is the mean of all local apertures inside fracture induced by stress. Under the complex stress conditions, mechanical aperture is evaluated by the following equation [15]:where is the initial mechanical aperture of fracture, is the increment of mechanical aperture due to normal stress, and is the increment of mechanical aperture induced by shear. The initial mechanical aperture directly uses the measured mean value of the vertical distance between the two surfaces before stress, which follows in Table 1.

The measured upper and lower fracture surface information is composed of the point cloud data. The points were imported to the preprocessor program of ANSYS 12.0 and built the 3D fracture meshes. For instance, Figure 3 shows the Gr3 mesh. The digitized fractures with deformations were generated by applying the normal stress and shear on upper surface. The relationship of normal stress and normal displacement is according to Bandis et al. [29]:where and are constants which depend on initial stiffness and maximum aperture. Kulhaway (1975) proposed the following formula representing shear and shear displacement relations:where and are constants which depend on initial shear stiffness and shear strength. During the applying stress, aperture distribution of fracture would be change and some asperities were overlapped. Overlapping asperities were assumed to be contacting asperities, and fracture surface deformations were not considered. Figure 4 shows the void space geometry at normal displacements (Figure 4(b)) and shear displacements (Figure 4(c)), respectively. The locations where apertures are equal to zero represented contact areas. The contact characteristics can be presented by the contact ratio , which is defined as the percentage of area fraction occupied by the contact areas to total areas: where [M2] is contact areas that include dead voids (i.e., Figure 4(a)). The area is no contribution to the flow [11], and the [M2] is total areas in fracture.

3.2. Flow Simulation

The flow simulation for rough rock fracture is conducted by solving the Navier-Stokes equations (NSE) using finite volume method (FVM). For a stable, incompressible, and low velocity fluid, the NSE can be written aswhere denoting the density of a fluid, the pressure, the viscosity coefficient, and velocity vector, respectively. In this study, the density and viscosity coefficient of water are used as 997.295 Kg/m3 and 0.001088 Pa·s in the room temperature of 25°C. The discrete equation is obtained and solved by integrating the governing equation in a series of finite volume, into which the calculate zones are divided. And the finite volume is represented by a node.

Before solving (7) using the FVM, these equations need to be adopted as the integral forms:where , , and are the normal vector, area, and volume, respectively. The discretized FVM formulation of the above integral equations in the nodes becomeswhere is the length increment of nodes and . To solve (9)~(12), the commercial program FLUENT (Ansys Inc., Canonsburg PA) is used to simulate flow processes in this study. The solved method of momentum and pressure is adopted using 2-order upwind scheme and 2-order scheme, respectively; and the pressure-velocity coupling is calculated by SIMPLE algorithm. The solution gradients are solved by direct interpolation with a least square method at the center of each cell [30].

Three types of boundary conditions are taken into account, which is the same as the one used in Xiong et al. [15] and Babadagli et al. [28] laboratory flow tests. The inlet boundary of fracture is set as a constant velocity inlet, and a pressure outlet is modeled on the opposite end of the fracture with constant pressure 0 (see in Figure 3). Both the fracture walls and two sides of fracture are assumed to be external no-flux boundaries. And the treatment of fracture contacts is as internal no-flux boundaries, with no flow into or out of the contact areas.

To accurately simulate flow behaviors, each fracture specimen is divided into approximately 1.5 million elements with an edge length of 0.3 mm, and most of the elements are hexahedrons. When the two opposing fracture surfaces are in touch, the contact zone is assigned with void element (Figure 4). Hence, the fracture asperity damage at the contact zones is considered by removing the overlapping parts during the normal and shear deformation.

4. Numerical Results

4.1. Mechanical and Hydraulic Behaviors during Normal Loading

Figure 5(a) shows the - curves, and Figure 5(b) shows the variation of contact ratio and mechanical aperture for fractures Gr1–Gr5 with normal displacements during normal loading, respectively. The contact ratio is shown using solid points, and yet the open points represent mechanical aperture. It is indicated that the contact ratio nonlinearly increases during normal closure in two stages. In initial stage, the contact ratio keeps slow growth and increases abruptly in the second stage. The contact areas range between approximately 0.0% and 38% at a normal displacement of 0 and 1.8 mm. On the contrary, the mechanical aperture has negative linear relationships with normal displacement, ranging from 0.6 mm to 3.0 mm.

The aperture distribution of fractures during normal closure is also plotted in Figure 6. The shape of frequency histograms remains unchanged, but the peak of frequency decreases, indicating that the mean aperture and its standard deviation decrease during closure. It is also observed that a special point of mechanical aperture jumps abruptly except for the Gr2 when normal displacement is more than 1.0 mm, signifying the contact ratio increases. The similar enlargement of contact areas was observed by Koyama et al. [12].

The flow simulation through fractures with normal displacement of 0–1.8 mm was performed. The transmissivity is calculated using (13) based on the numerical results.

Due to -direction flow rate being far less than -direction value, the -direction flow rate was used to evaluate the transmissivity. The transmissivity of rock fractures under 7 m3/s during normal closure is plotted, as shown in Figure 7. For each fracture, the transmissivities decrease with normal displacement increases. Two increasing stages are detected, ranging from 0.0 to 0.8 mm and beyond 0.8 mm (Figure 7), just like the changes of contact ratio with normal displacement. The transmissivities gradually decrease during the first stage and then keep decreasing more than two orders but with high gradients as normal displacement increases from 0.8 mm to 1.8 mm (second stage). The values of transmissivity in each fracture are in the range of to 0 m2/s, of which the difference is influenced by the roughness of fracture surface. A strong conformity between our results and experiment results in previous literature is noticeable [16, 17, 31].

Figure 8 shows the results of flow simulation in the Gr3 with 0.5 mm, 0.8 mm, and 1.0 mm normal displacement (corresponding to contact ratio of approximately 7.01%, 15.75%, and 24.76%). The white area indicates the area in contact, where the velocity equals 0. It is clearly identified that the predominant flow paths occur in the cases of 0.5 mm and 0.8 mm normal displacement. As increasing normal displacement (see in Figure 8(c)), the previous flow paths are broken down, and the emergence of velocity concentration and back region leads to decrease in transmissivity, due to increasing contact areas.

4.2. Mechanical and Hydraulic Behaviors during Shear

The digital upper surfaces were continuously displaced tangentially with a 1 mm interval to simulation shearing without normal loading for the specimens Gr1, Gr3, and Gr5. Their aperture distribution was then determined at each shear displacement. Figure 9 shows the aperture distribution of three samples at different shear displacements. For all samples, the shape of frequency histogram becomes flatter during the increasing shear displacements, indicating that the mechanical aperture and standard deviation increase. Similar behavior was also reported in Koyama et al. [12]. It can be seen that the frequency of special zero point keeps increasing with increasing shear displacements, indicating expansion of contact areas. This is most likely due to the effect of the fracture roughness. Some statistical mechanical parameters of those rough fractures subjected to shear and shear loading curves are shown in Figure 10. The results show that, in all the shear conditions, the mechanical aperture linearly increases, while the contact ratio nonlinearly increases with increasing shear displacements. Watanabe et al. [32] reported that the contact areas increase significantly beyond shear displacements of 1 mm. The same results are also seen in Xiong et al. [15].

The FVM method was adapted to study flow behavior during the shear. The injection flow rate is parallel to the shear displacement with  m3/s. The relations between transmissivity and shear displacement obtained by all flow analyses are shown in Figure 11. The results show that transmissivity slowly decreases with shear displacements for all fracture samples. And the decrease of transmissivity with the shear displacement of 0–4 mm varies with geometry of fracture, which is about 1.5–2 orders of magnitude. This result is very close to the result of normal deformation conditions.

To clarify the change of fluid transportation capability of fractures during shear, the images of flow velocity fields can be taken as evidence. Figure 12 shows the flow velocity field in the specimen Gr3 at shear displacement of 1.0 mm, 2.0 mm, and 3.0 mm. It is implied that no flow channel was detected in the numerical model during initial stage (Figure 12(a)), due to small contact areas in fractures. The contact areas block off fluid flow and decrease the transmissivity of fracture during shear.

4.3. Evolutions of Mechanical Aperture and Hydraulic Aperture

Based on the pressure and flow rate data, the value of was calculated by the cubic law (). In order to understand effect of the fracture roughness and the contact area, Figure 13 depicts a three-dimensional plot of the simulative data in the form of the normalized hydraulic aperture (/) against fractal dimension and contact ratio c, which indicates a clear correlation between the hydraulic aperture and two variables. Hence, a proper hydraulic aperture model needs to incorporate both and . Motivated by the characteristics of the curves plotted in Figure 13, the empirical power law expression of can be accurately established based on Zimmerman and Yeo [4] model and fractal power relation. It can be written aswhere , , and are regression coefficients. The Levenberg-Marquardt nonlinear method was used in (14). This afforded , , and , and the equation had best correlation () with that data. Therefore, the new equation for can be expressed asThis model is useful for prediction of hydraulic aperture in real fractured rock that has both rough void area and contact obstacles. It can be inferred from (15) that the hydraulic aperture has negative correlation with both contact ratio and fractal dimension, and the hydraulic aperture is less than mechanical aperture due to fracture surface roughness and contact characteristics. It explains the reduction of flowrate by interaction of contact obstacles and roughness in fracture Darcy flow.

In order to assess the effect of fracture morphology, Figure 14 plots the predicted variation of normalized hydraulic aperture (/) with fractal dimension at different values of the contact ratio (%, 16%, 24%, 32%, 40%, and 48%) using (15). The result implies that the value of decreases drastically with at the same . When decreased by two orders of magnitude, the values of (/)3 experience an decrement by 2.3 times. While the values of decrease only 1–1.5 times, the increases by two orders of magnitude. It is demonstrated that the influences of fractal dimension are smaller than contact ratio. Hence, the importance of incorporating the fracture contact ratio in the development of the hydraulic aperture models is illustrated [4, 11].

5. Discussion

5.1. Comparison of Models

The variation of hydraulic aperture by Louis (1969), Zimmerman and Yeo [4], Matsuki et al. [13], and Xiong et al. [15] and the proposed model are listed in Table 2. Louis (1969) model only considers the effect of small relative roughness and large opening in fractures. By using the effective medium theory and experimental data, Zimmerman and Yeo [4] proposed a hydraulic model related to the contact ratio. The effect of fracture geometries on hydraulic aperture has been investigated numerically by Matsuki et al. [13] and experimentally by Xiong et al. (2010). However, a larger relative roughness commonly corresponds to a mated rock fracture with two confined fracture surfaces tightly in contact with numerous contact points. These contact points result in tortuosity of the flow paths and much larger flow resistance. Both the Zimmerman and Yeo [4] and the new proposed model can reflect the contact areas effect.

5.2. Model Evaluation

To quantitatively evaluate the existing hydraulic aperture models using the experimental data, the Normalized Objective Function (NOF) and the slope of the predicted versus experimental values of the hydraulic aperture are adopted. NOF is defined as the ratio of the root mean square error to the mean of the experimental data:where is the experimental value of hydraulic aperture, is the predicted value of hydraulic aperture models, and is the total number of the data points. The slope is determined by performing a linear regression between the predicted and experimental values of hydraulic aperture. The optimum values for the model evaluation are and (Zhou et al. 2016).

The values of NOF and , which are evaluated using the numerical data, are listed in Tables 3 and 4. It shows that the NOF values of the proposed model are below 0.6 and its values are close to 1.0 for five fracture specimens, indicating a good match between the model predictions and numerical data.

The NOF and values for the specimens are also plotted in Figures 15 and 16. It seems that the model is the most reliable one for estimation of the hydraulic aperture, compared to other models. Although the aperture values predicted by Louis model and Xiong model agree well with experimental results, especially for the Gr3, the results from proposed model for the hydraulic aperture prediction are better than those from Louis and Xiong models. On contrast, Zimmerman and Matsuki models seem to deviate from the experiment data, especially the NOF > 1.2 for Fr2 and Fr4. The reason may be that the models are developed based on the smoother fractures with , compared to the of the granite fractures.

5.3. Validation Based on Other Experimental Data

Laboratory water flow experiments were performed to investigate the effect of roughness of rock fracture on flow by Develi and Babadagli [16]. The findings showed that the transmissivity was related to the fractal dimension and contact ratio. The relationships of transmissivity with fractal dimension are plotted in Figure 17. Note that the values of mechanical aperture were not provided in Develi and Babadagli [16]. The proposed model was fitted to the experiment data by taking the as variable. And according to the cubic law, the transmissivity is easily known by making use of hydraulic aperture from (14) as follows:Applying the nonlinear fitting method, the best-fitted is 0.5 mm. Figure 17 plots the curve-fitting of the experiment dada with the transmissivity based on the proposed model, and it is demonstrated that the new model can capture the variation tendency of with , despite the great scatterness.

5.4. Limitations of Proposed Model

In fact, the new model fails to predict the hydraulic aperture when the contact ratio is greater than 57%, to ensure that (14) is greater than 0 identically. Watanabe et al. [14], however, pointed out that the higher contact ratio exists inside the natural rock fractures. In addition, Xiong et al. [15] considered that the variation of hydraulic aperture is related to the inertial effect. The transmissivity decreases with increasing Reynolds number (), but the decrement is small, especially in low flow rate. The reason is that the viscous effect is predominant in low flow rate stage. Zimmerman et al. [33] believed that the existence of a weak inertia regime for Reynolds numbers was in the range of 1–10. In the study, Reynolds number is 1.6, according to . There is no effect of inertial on flow behavior in the range. So the proposed model fits well the numerical data, without incorporating Reynolds number. Permeability anisotropy is induced by the shear displacement when flow is parallel and perpendicular to shear direction [21]. In the study, we only consider the flow along the shear direction. And the proposed model can accurately characterize these flow behaviors. However, a more generalized model could be developed considering the Reynolds number and anisotropy.

6. Conclusions

The hydraulic and geometric characteristics on rough-walled fractures were studied using numerical simulation. The geometry of five artificial-splitting granite fractures was measured by 3D scanning technology. Flow behaviors through granite fractures were modeled using FVM during normal loading and shear process. On the basis of these works, the geometrical characteristics and corresponding hydraulic characteristics of fracture are analyzed:

(1) The mechanical aperture of fracture decreases from 0.6 to 3.0 mm with increasing normal displacements, while contact ratio increases from 0% to 38%. The transmissivity decreases by 3 orders of magnitude, as normal displacements increase, due to the change of void space geometry.

(2) When the fracture shears from 0.0 to 3.0 mm, the aperture shows a growing momentum, but with a much smaller gradient. And the performance of hydraulic conductivity declines by 1.5–2.0 orders of magnitude.

(3) By incorporating the contact effect, a new practical model (14) is proposed for evaluating the hydraulic aperture, using fractal dimension and contact ratio. The evaluation of effects of fracture geometry on the Darcy flow demonstrates that the fractal dimension and contact ratio play an important role in permeability of roughness rock fracture, while the influence of contact ratio is more dramatic.

(4) Verification of the new model is performed by comparing with five existing models and illustrates that the estimated values of hydraulic aperture fit well with the numerical results. The new model is found to have good accuracy for prediction of the hydraulic conductivity.

(5) Further studies will focus on modifying the proposed model by applying it to more natural fractures with different lithological rock and surface characteristics and high Reynolds number.

Conflicts of Interest

The authors declare that there are no conflicts of interest regarding the publication of this paper.

Acknowledgments

This research is financially supported by the National Science Foundation of China (nos. 51679173, U1765207, and 51709207) and the Natural Science Foundation of Hubei Province (no. 2016CFA083). This support is gratefully acknowledged.