Abstract
The effects of roughness and normal stress on hydraulic properties of fractures are significant during the coupled shear flow test. Knowing the laws of fluid flow and solute transport in fractures is essential to ensure the nature and safety of geological projects. Although many experiments and numerical simulations of coupled shear flow test have been conducted, there is still a lack of research on using the full Navier-Stokes (N-S) equation to solve the real flow characteristics of fluid in three-dimensional rough fractures. The main purpose of this paper is to study the influence of roughness and normal stress on the fluid flow and solute transport through fractures under the constant normal stiffness boundary condition. Based on the corrected successive random addition (SRA) algorithm, fracture surfaces with different roughness expressed by the Hurst coefficient () were generated. By applying a shear displacement of 5 mm, the sheared fracture models with normal stresses of 1 MPa, 3 MPa, and 5 MPa were obtained, respectively. The hydraulic characteristics of three-dimensional fractures were analyzed by solving the full N-S equation. The particle tracking method was employed to obtain the breakthrough curves based on the calculated flow field. The numerical method was verified with experimental results. It has been found that, for the same normal stress, the smaller the fracture value is (i.e., more tough the fracture is), the larger the mechanical aperture is. The ratio of hydraulic aperture to mechanical aperture () decreases with the increasing of normal stress. The smaller the value, the effect of the normal stress on the ratio is more significant. The variation of transmissivity of fractures with the flow rate exhibits similar manner with that of . With the normal stress and value increasing, the mean velocity of particles becomes higher and more particles move to the outlet boundary. The dispersive transport behavior becomes obvious when normal stress is larger.
1. Introduction
In some rock masses with low permeability, fractures are the main channels for fluid flow. The hydraulic characteristics of fractures are of great concern in some rock mechanics and geotechnical applications, such as nuclear waste disposal, geothermal energy mining, and deep mineral mining. The fluid flows through the fractures while the solutes or particles also move with the fluid by the convection and diffusion mechanisms. Understanding the laws of fluid flow and solute migration behavior in fractures is essential to ensure the safety of these geological engineering.
The fracture surface of the natural rock mass is generally rough [1], causing the rock fractures to be composed of void spaces and contact areas, rather than many studies assumed that fractures are composed of two relatively smooth parallel plates. Fluid bypasses the contact area and flows through the void space with tortuosity. The characteristics of void space geometry and the contact area distribution have significant effects on the hydraulic properties of a fracture [1]. It has been demonstrated that the fluid flow and transportation process in a single fracture are heavily influenced by the roughness of the fracture surface, and the mechanical aperture of a rock fracture is usually larger than its hydraulic aperture [2]. Li et al. [3] implemented a series of coupled shear flow tests to analyze the influence of geometric features of fracture on rock mechanical behavior and proposed the empirical correlations to evaluate the effects of surface roughness and contact area on the behavior of fluid flow through rough fractures. Zou et al. [4] investigated the effects of wall surface roughness on fluid flow through fractures; the result indicated that the flow rate and the roughness of the fracture surface are the main reasons for the dynamic evolution of the vortex area in the flow field; when the flow rate is high, the fluid flow field is usually not only nonlinear but also constantly produces eddies in the boundary layer region of the rough fracture which will affect the solute transportation in fractured rock masses. Wang et al. [5] adopted 3D lattice Boltzmann method and combined wavelet analysis technique to investigate the impact of surface roughness on the nonlinear fluid flow 3D rock fractures. The result shows that the primary roughness mostly controls the pressure distribution and fracture flow paths at a large scale, whereas the secondary roughness determines the nonlinear properties of the fluid flow at a local scale.
For a single rock fracture, the roughness of the upper and the lower surfaces is the same under the initial conditions, and then, the deformation occurs with the normal stress and shear stress applied to the natural rock [6–8]. In order to study the effect of shear on fluid flow in fractures, some scholars have carried out a series of experimental and numerical simulation research on coupled shear flow. The laboratory research on the shear behavior of fractures is usually implemented under the condition of constant normal stress (CNL). However, for many field situations, the normal stress imposing on the fracture surface will change during the shearing process. The expansion of fractures is usually constrained by the closed environment which is represented by the constant normal stiffness (CNS) [9]. Therefore, it is necessary to study the shear characteristics of the fractures under the boundary condition of constant normal stiffness. Indraratna et al. [9] studied the effect of initial normal stress levels on shear behavior of joints under CNS conditions. The result revealed that the initial normal stress has significant effect on the shear dilation rate. Different initial normal stresses will lead to the variety of the distribution of voids and contact areas in the fractures under shear.
In rock mechanics and rock engineering practices, the fractures are usually simplified to two parallel plates in which the fluid flow follows the cubic law. However, the parallel models are inadequate to describe the hydraulic and transport properties of natural fractures with rough surfaces [10]. To further take into account the geometrical characteristics of fractures, the simplified forms of the full Navier-Stokes equations, such as the Reynolds equation [11–17] and the Stokes equations [18, 19], were used in the estimation of hydromechanical properties of rock fracture. The flow rate and Reynolds number in natural rock fracture are not always small, and the inertial effects usually increase with the complexity of void space geometry; therefore, the nonlinear term is not always negligible [4, 20–25]. Without considering the influence of inertia effects, the simplified form usually shows deviation from the actual situation. In order to accurately describe the fluid flow in the fractures, numerical simulation by solving the full N-S equations should be adopted, especially for rough natural fracture.
Many coupled shear flow experiments have been carried out. These studies mainly focused on using the simplified N-S equation to describe the fluid flow in a two-dimensional fracture. And there is still a lack of research on using the full N-S equation to solve the real flow characteristics of fluid in three-dimensional fractures. It should be noted that most of published research considered only simplified fracture surface topography or specific boundary conditions (e.g., constant normal stress boundary). A comprehensive study of the impact of normal stress and roughness on fluid flow and particle transport under the constant normal stiffness boundary condition was rarely presented in existing publication. In this study, several 3D fracture surfaces with different roughness coefficients were obtained at first. The fracture used for numerical simulation was composed of two fracture surfaces with the same roughness and the initial aperture was zero. A series of shear tests under the constant normal stiffness with different normal stress were then employed. Finally, a number of numerical simulations by solving the full 3D Navier-Stokes equations were adopted to investigate the fluid flow behavior and solute transportation through fractures.
2. Geometry Model of 3D Rough Fractures
2.1. 3D Self-Affine Fracture Surface
The natural rough rock fracture surface typically follows a self-affine fractal distribution [26–28], and the topography of a fracture surface can be characterized by fractional Brownian motion (fBm) [29]. Several methods such as the Fourier transformation, the randomized Weierstrass-Mandelbrot function, and the successive-random addition (SRA) were widely used to model the fBm [30–33]. In this study, the efficient SRA algorithm, which is easy to understand and use, is adopted to generate rough fracture surface.
In fBm, the surface asperity height is defined as a random and single-valued function of two independent spatial variables, and . The stationary increments, , over the distance displays a Gaussian distribution with zero mean and variance [5, 34]. The statistical self-affinity of fBm can be expressed as follows: where represents the mathematical expectation, and is the roughness exponent or Hurst exponent varying from 0 to 1 and related to the 3D fractal dimension () by .
Liu et al. developed a corrected SRA algorithm to generate the 3D self-affine fracture surfaces which overcomes the problems associated in the traditional SRA algorithms that provide stochastic fractal distributions of questionable scaling and correlation properties [30, 35]. In this study, a series of rock fracture surfaces, with , 0.55, and 0.6 (the corresponding joint roughness coefficient () was calculated using Equation (2) [36] and the values are 18.2940, 12.3983, and 6.5396), respectively, were generated by employing Liu’s algorithm. These surfaces are both 204.8 mm in length and width. As shown in Figure 1, the fracture surface asperity becomes flatter with the increase of . where is the number of intervals, is the interval, and is the height of surface.
2.2. Aperture Distribution of Fractures under Different Normal Stress
A lot of research revealed that the aperture which is defined as the distance between the two fracture surfaces has a significant influence on the hydromechanical characteristics of the rock fractures [37–41]. In this study, it is assumed that the two fracture surfaces are in contact with each other at the initial condition. The aperture between the surfaces was assigned zero. When shearing occurs, the fracture surfaces separated in the horizontal and vertical directions. Under a specific shear displacement, the aperture distribution can be calculated by the following formula [37]:
In the above, is the asperity height of the generated 3D self-affine fracture surfaces, is the shear displacement, and is the normal displacement or dilation caused by the shear displacement.
Research on the fracture shear behavior in the laboratory was usually carried out under constant normal stress boundary conditions, where the normal stress always remains constant and the rock joint dilates freely during shearing [9]. However, some researchers have pointed out that the joint dilation may be constrained by a restricted environment, which usually represents a constant normal stiffness condition (CNS) in engineering practice and the CNS boundary conditions are more applicable to many areas [9, 42–45]. In this study, the following analytical model of dilation and shear displacement under CNS conditions proposed by Indraratna et al. [9] was used to calculate the dilation caused by shear displacement. where is the shear displacement and is the dilation rate that change with the ratio of shear displacement to peak shear displacement () for a fracture subjected to shear under CNS boundary condition.
The value of dilation rate can be calculated by the following equation.
Here, and are decay constants which can be calculated from experimental data. is peak shear displacement. The value of at which the dilation rate is assumed to begin is about 0.3 for rough fracture [14, 46]. is the peak dilation rate and can be obtained from the following equation.
In the above, is the joint roughness coefficient, is the compressive strength of the joint surface, is the damage coefficient that is given a value of 1 or 2 under low normal stress or high normal stress, is the CNS at an external boundary, is the initial normal stress, is the initial joint normal stiffness at zero normal stress level, and is the maximum closure of the joint. Related parameters which are used in numerical simulation are listed in Table 1.
Either the joint roughness coefficient or the normal stress can significantly affect the aperture distribution under the shearing. Limited by the challenging computational capacity of solving the N-S equations for the 3D roughness fractures with high precision, only a square area with a side length of 27.6 mm on the surface which is in a certain range of mm and mm was cut out to form the 3D models. Fractures exhibit poor connectivity, and no obvious fluid flow can be observed with small shear displacement. Finally, 5 mm was selected as the shear displacement, with which some of the main flow channel formed and the contact area was more concentrated. Figures 2–4 show the aperture distribution of different roughness fractures with the shear displacement of 5 mm under different normal stress.
3. Fluid Flow and Solute Transport Simulation
3.1. Fluid Flow Simulation
3.1.1. Governing Equation
The governing equation for a single fluid flow in a fracture is Navier-Stokes equations, which are derived from Newton’s second law and are strict statements of the momentum conservation. For the steady-state and incompressible single Newtonian fluid, the Navier-Stokes equation can be expressed in a vector form as follows: where is the flow velocity which contains three velocity components in 3D rock fracture, is the kinematic viscosity of fluid defined as , is the fluid density, is the pressure, and represents the body forces acting on the fluid.
The N-S equations are a set of nonlinear partial differential equations coupled with velocity and pressure fields [20]. The left of the N-S equations contains a nonlinear convection acceleration term, and the right has a viscous diffusion term, so the governing equation of viscous fluid is a nonlinear convection diffusion equation. It is difficult to use theoretical methods to obtain accurate solutions for complex three-dimensional flows, except for a few simple fluid flows which can obtain analytical solutions. In this study, the commercial FEM software of COMSOL Multiphysics was employed to simulate the fluid flow through the roughness fractures under the shearing. The density and viscosity of water at 10°C were taken as and .
3.1.2. Boundary Conditions
In this study, the axis direction was selected as the main fluid flow direction. The two boundaries at and were set as the inlet boundary and outlet boundary. The combination of flow rate and pressure boundary conditions was used. Four laminar entrance inflow conditions with different constant flow rates were adopted at the inlet boundary to investigate the effect of inlet flow condition on the hydromechanical properties of rough fractures and the solute transport in fractures. At the same time, set the pressure at the outlet boundary to zero. The remaining upper, lower, and the two side boundaries were set to no flow and no slip, where the fluid velocity relative to the walls’ velocity is zero. Such boundary conditions are consistent with many laboratories’ experiment conditions [20].
3.2. Solute Transport Simulation
The particle tracking approach was adopted to simulate solute transport based on the flow result obtained from steady-state flow field. In this study, only advection process was considered, and the particles were driven by the fluid and transported along the fluid flow path. The random dispersion due to diffusion of the solute particles within the fluid in fractures and other retardation mechanisms such as sorption or decay were not taken into account [17]. The particle tracking method in the fluid was used to compute the motion of particles in a background fluid. The particle momentum comes from Newton’s second law, which states that the net force on a particle is equal to its time rate of change of its linear momentum in an inertial reference frame. The particles in the fluid are driven by drag force and its momentum can be described as follows:
Here, is the particle mass, is the particle velocity response time, is the velocity of the particle, is the fluid velocity, and is the drag force.
Once the flow velocity was calculated element by element by solving the N-S equation, all the particles travel following the streamlines. Particles were initially placed at one edge of a FEM element. They follow the velocity of each element. The travel time in each element was given by magnitude of the velocity of the element as follows [47]:
In the above, is the travel time of particle in tracking step , is the travel distance of particle inside the element corresponding to the tracking step , and is the magnitude of the velocity vector of the element corresponding to the tracking step .
The total residence time for a particle is the sum of the time of all the elements that the particle is passing through. where is the total travel time of particle , and is the number of tracking step for particle .
The number of particles injected at each position along the entrance of the fracture is proportional to the velocity. This means that more particles will be attracted to places with higher flow velocity. From the particle tracking simulation, the breakthrough curve that represents the relationship between the percentage of particles collected at the outlet boundary and time can be obtained using the tracking time of each particle. The Peclet number can be defined in terms of the variance and mean travel time using the following equation [48]. where and are the variance and mean travel time, respectively.
4. Result and Analysis
Three self-affine fracture surfaces were generated with , 0.55, and 0.6, respectively. The fracture models were composed of two fracture surfaces with the same roughness coefficient at the initial state. For fracture models with different , a shear displacement equaling to 5 mm was used and three normal stress that equal to 1 MPa, 3 MPa, and 5 MPa were loaded on the models with CNS boundary in the shear simulation. For each model, the normal displacement was calculated using the Equation (4). The final void distribution of each model was different due to the application of different normal stresses. Four flow rates that equal to , , , and were adopted at the inlet boundary. The corresponding Reynolds (Re) numbers are 0.4, 4, 40, and 400 which exist in ground flow in rock engineering, hydraulic engineering, and laboratory experiments [20].
4.1. Characteristics of Aperture Distribution
Distribution characteristics of aperture have an important influence on the simulation of fluid flow in fractures. The evolution of fracture aperture can be calculated with the values of asperity height and the normal displacement during shear. Figure 5 shows the relationship of the variations of normal displacement and the normal stress for fractures with different roughness coefficients. Both roughness and normal stress can significantly affect the normal displacement. For all fracture models, the normal displacement () generally decreased as the normal stress () increased. At the same time, as the magnitude of increased, the normal displacement () became smaller. In other words, the rougher the fracture surface was, the greater the normal displacement () could be obtained under the same normal stress during the shear. This can be attributed to the fact that the peak dilation rate decreases with an increasing normal stress but increases with the increasing which represent the roughness coefficient of surface. The larger the value is, the rougher the surfaces are [9]. Many studies have shown that the aperture of fractures usually obeys a normal distribution [20, 49]. According to the obtained apertures, Gaussian fitting was performed on the aperture of each fracture model. Some statistical parameters (the mean aperture and the standard deviation) of local mechanical aperture under different normal stress with different are listed in Table 2. Figure 6 is the aperture density distribution of fracture under normal stress 1 MPa with value equaling to 0.5. Both mean value and its standard deviation are decreased with the increment of normal stress and increased with the roughness.
Figures 2–4 show the spatial distribution of fracture aperture. The local apertures were discretely distributed in space, but showed good connectivity. For a fracture with higher roughness at a lower normal stress, there was no contact area between opposite fracture surfaces. In contrast, some contact areas (small white squares representing zero aperture) were appeared gradually with lower roughness at a higher normal stress. Small aperture zones were shown around the contacting area and fluid flow through the main channel with higher aperture.
4.2. Hydraulic Characteristics of the Fracture
Generally, the calculated hydraulic aperture was lower than the mean mechanical aperture due to the existence of contact areas within the fractures and the tortuosity of streamlines. The complexity of the void geometry decreased as the shear displacement increased. At the same time, more flow channels emerged accompanied with the decrease of contact ratio. When the shear displacement became larger, the void geometry became more like parallel plates model with large mechanical aperture. While choosing 5 mm as the shear displacement, the calculated ratio of hydraulic aperture to mechanical aperture is basically from 0.8 to 0.975 during the shear of both three fractures (Figures 7–9). As shown in Figure 7, the rougher fracture with equaling to 0.5 exhibits larger decrease comparing to other fractures with lower value since the rougher fracture could produce more complex void geometries during the shearing. Lower normal stress produces larger shear dilation and results in larger local aperture. It can be seen from Figure 7 that the value of decreases with the increment of .
The variation of transmissivity of fractures with the flow rate exhibited similar manner with that of , since the transmissivity is proportional to . Figure 10 elucidates the relationship of transmissivity and flow rates with under different normal stress. It can be observed that the measured transmissivities under larger normal stress have smaller values and decrease with the flow rates, which indicates that transmissivity not only changes significantly with the increasing of flow rates but also affected by the normal stress. Figure 11 shows the relationship of transmissivity and flow rates with .
Experiment data from Xiong et al. [1] were used to verify the proposed numerical method. Various flow rates from Xiong’s experiment were imposed on the inlet boundary of fracture with value 0.6. The normal stress was set as 1.5 MPa and the shear displacement adopted was 8 mm according to Xiong’s experiments. Due to lack of original experimental data on the fracture geometry, the fracture was randomly generated with . Figures 12 and 13 show the results from experiments and simulations, respectively. It could be seen that the simulation results show the same trend with the experimental data and the values of transmissivity are close. The difference between experimental and simulation values is very likely caused by the difference of fracture geometry in the simulations and the experiments.
4.3. Solute Transport Simulation
In the previous section, flow velocity fields were calculated by using the N-S equation. Particles transported following different velocity trajectories with different numbers. The effects of flow velocity through the fractures on the particle movement characteristics were investigated in this section.
According to the simulation results of particle tracking, the breakthrough curves were obtained. As shown in Figure 14, the breakthrough curves shift to the left without changing their shape when normal stress becomes larger. The mean velocity of particles becomes higher, and more particles move faster to the outlet boundary with larger normal stress. When , a long tail appears on the breakthrough curve and the number of particles finally reaching the outlet decreases with the increment of normal stress, which can be explained by the occurrence of contact areas. For the fracture with , the contact areas gradually appear with the increment of normal stress, which leads to the preferential flow field. Particles travel through the preferential flow. Such preferential transport behavior will result in earlier arrival of particles through the high-velocity zone in fracture, but heavy tailings of particles spread in the low-velocity zones around the contact areas [50]. The very low transmissivity region increases with the normal stress and more particles get trapped in this zone. With increasing roughness, the inclination of the breakthrough curves becomes larger. The preferential flow phenomenon becomes more obvious with larger value and particles travel through the channels with higher velocity.
(a)
(b)
(c)
The mean and standard deviations of particle travel time were calculated using results obtained by solving the N-S equation at the identical percentage of particle collection at the outlet (Table 3). Using the statistical parameters of travel time, the Pelect numbers (Pe) for characterizing the dispersion of the transport can be obtained by Equation (12). The relationship between Pe and normal tress presents the similar evolution trend under different flow rates. Figure 15 shows that the Pe numbers decrease with increasing normal stress and value. In other words, the dispersion () increases with increasing normal stress and value, since Pelect number (Pe) has an inverse relationship with dispersion (, where is the fracture length). One interpretation for this behavior may be that the preferential flow channels gradually appeared with the increasing normal stress and value. The formation of preferential flow results in more discretization velocity field and finally causes the dispersive transport behavior of particles.
5. Conclusion
In this paper, several 3D self-affine fracture surfaces were generated at first by using the SRA method and then the effects of surface roughness and normal stress on fluid flow and solute transport in fractures were investigated through a series of coupled shear flow simulation. The fluid flow was simulated by solving the N-S equations, and the solute transport was simulated by the particle tracking method with the fluid velocity fields predicted by solving the N-S equations. The results revealed that the roughness and the normal stress may have significant influence on the fluid flow field and the residence time of particle, which is an important issue that needs to be considered in the safety assessment of radioactive waste repositories in fractured crystalline rocks.
Roughness may affect the tortuous degree of flow channel and make a sudden change of local mechanical aperture which will make a profound impact on fluid flow and particle transport. The normal stress mainly affects the fracture aperture which can change the flow field in turn and finally affects the travel path and residence time of particles.
Data Availability
All the data used in the present study are available by contacting the corresponding author and/or the first author. Their email addresses are [email protected] (Fenhua Ren) and [email protected] (Min Wang), respectively.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
This study has been supported by the Key Laboratory of Western Mine Exploitation and Hazard Prevention, Ministry of Education (SKLCRKF1901) and the National Natural Science Foundation of China (51804163). These supports are gratefully acknowledged.