Effect of Shearing on Non-Darcian Fluid Flow Characteristics through Rough-Walled Fracture

Abstract

The heterogeneous fracture geometry induced by the presence of roughness and shearing complicates the fracture flow. This paper presents a numerical investigation of the non-Darcian flow characteristics of rough-walled fractures during shear processes. A series of fracture flow simulations were performed on four types of fractures with different joint roughness coefficients (JRCs), and the different shear displacements were imitated by degrees of mismatch on two fracture surfaces. The results show that the disorder of fracture geometries and the increase in flow rate are the main causes for the emergence of an eddy flow region, which can significantly reduce the fracture conductivity and change the fracture flow from linear to nonlinear. The Forchheimer equation provides a good model for the nonlinear relationship between the hydraulic gradient and the flow rate in the fracture flow. When the shear displacement or JRC increased, the linear permeability coefficient $k_v$ decreased, while the nonlinear coefficient $\beta$ increased. A three-parameter equation of $\beta$ was used to examine the inertial effect induced by the fracture roughness JRC and the variation coefficient of aperture distribution ${\sigma }_s/e_m$. The critical Reynolds number was a combined effect of aperture, viscous permeability, and inertial resistance, assuming the flow becomes non-Darcian when the inertial part is greater than 10%.

1. Introduction

The fluid flow processes in fractured aquifers are of great importance in many engineering practices, such as water flow in dam foundation, underground mining extraction, enhanced geothermal systems, and fluid waste disposal. Natural rock masses commonly comprise matrix blocks with fractures that provide primary migration pathways for fluids. Understanding the fluid flow behaviors in rock fractures is essential for simulation of fluid flow in fractured rock masses.

Most fluid flow models were developed based on the assumption of the linearity of Darcy’s flow, which can only be anticipated for laminar flow by diminishing the inertia effects [1,2,3]. The inertial losses cannot be neglected with regard to viscous forces arising from changes in flow flux or direction and localized eddy formation [4,5]. In this case, the pressure drops more than the proportional increase in the flow velocity; this is known as the non-Darcian fluid flow. To describe such viscous and inertial losses, the Forchheimer equation is widely used [6,7,8].

$$-\nabla p=Av+B{v}^2$$

(1)

where $\mathit{\nabla }p$ is the pressure gradient, $v$ is the volumetric flow velocity, and $A$ and $B$ are the coefficients describing viscous and inertial losses, respectively. To incorporate the Forchheimer equation into the numerical modeling of non-Darcian flow, it is essential to determine the coefficients $A$ and $B$.

The nonlinear flow behavior in a rock fracture influenced by shearing has been investigated using numerical simulations [9,10,11,12] and laboratory experiments [13,14,15,16,17]. Geometry alterations related to shearing and surface roughness have resulted in an increase in the overall resistance [12]. Several geometric parameters, such as fractal dimensions [18], the high-order statistical parameters of the roughness height [19], the effective slope [20,21], the roughness decomposition [5], and the joint roughness coefficient (JRC) [22], have been used to represent the wall roughness and to better understand the effect of surface roughness on the fracture flow. However, the surface roughness alone is not sufficient to represent the influence of geometry alteration. Shearing combined with roughness can result in a highly complex geometry alteration and can cause significant change to the fracture channel. The heterogeneity of the fracture channel can block the streamlines and cause backflows or eddies that indirectly increase the resistance of water flow in fractures and reduce the fracture conductivity [23]. The heterogeneous fracture channel eventually changes the two coefficients in the Forchheimer equation that describe energy losses [24]. Thus, the mechanistic understanding and parametric study of non-Darcian flow behavior for rough-walled fractures with shearing remains an open issue.

In this paper, we evaluated the effect of shearing on non-Darcian flow characteristics in rough-walled fractures. A series of numerical simulations of fracture flow were performed on four standard fractures, with different JRCs subject to different shear displacements. The evolution of the geometric characteristics of the fracture during shearing was estimated and discussed quantitatively. Fluent software was adopted to solve the fluid flow problem in these rough-walled fractures, and the occurrence and development of eddies were evaluated to gain microscopic insight into nonlinear flow processes. The relationship between the pressure gradient and flow rate was examined using the Forchheimer equation, and the impacts of shear displacement and JRC on the Forchheimer coefficients and the critical Reynolds number were analyzed.

2. Numerical Study

2.1. Fracture Model

Four fracture profiles from 10 standard fractures with different JRCs [25] were used for the roughness characterization and fluid flow modeling (Figure 1). The length was 100 mm and the JRC values were calculated using Equations (2a) and (2b) [26]. The upper surface was created by lifting up the lower profile for a vertical separation of 1.5 mm, leading to a tiny parallel channel for fluid flow.

$$Z_2={\left[\frac{1}{L}{{\displaystyle \int }}_0^L{\left(\frac{dy}{dx}\right)}^{2}dx\right]}^{1/2}={\left[\frac{1}{L}{\displaystyle \sum }_{i=1}^{N-1}{\left(\frac{y_{i+1}-y_i}{x_{i+1}-x_i}\right)}^2\right]}^{1/2}$$

(2a)

$$\mathrm{JRC}=32.2+32.47log{Z}_2$$

(2b)

To generate fracture models with different shear displacements, the upper surfaces were shifted in the x direction up to 2.0 mm, with the lower surface fixed. As shown in Figure 2, the red dotted line was the upper surface for the fracture model with shear displacements of 2.0 mm. It should be noted that shear mismatching simulations have been used in many studies [10,12,27,28]. Compared to the initial condition, the values of vertical aperture along the fracture length changed.

The identification of the geometric profile of the fracture channel was required for quantitative evaluation of fluid flow in rock fractures. Figure 3, Figure 4, Figure 5 and Figure 6 show that Gaussian distributions appeared to be appropriate for describing the frequency of vertical apertures with correlation coefficients of more than 0.8. The detailed statistical information for vertical apertures is listed in

Table 1. Statistical parameters for Gaussian fitting of vertical aperture with different fractures.

JRC	Shear Displacement ds (mm)	$\mathbf{Mean}\mathbf{Vertical}\mathbf{Aperture}{\mathit{e}}_{\mathit{m}}\left(\mathbf{mm}\right)$	$\mathbf{Standard}\mathbf{Deviation}{\mathit{\sigma }}_{\mathit{s}}\left(\mathbf{mm}\right)$	$\mathbf{Variation}\mathbf{Coefficient}{\mathit{\sigma }}_{\mathit{s}}/{\mathit{e}}_{\mathit{m}}$
5.34	0.5	1.497	0.0665	0.0444
	1.0	1.494	0.1193	0.0799
	1.5	1.491	0.1644	0.1103
	2.0	1.488	0.2040	0.1371
10.02	0.5	1.488	0.0947	0.0636
	1.0	1.477	0.1716	0.1162
	1.5	1.466	0.2341	0.1597
	2.0	1.454	0.2874	0.1977
13.51	0.5	1.494	0.1214	0.0813
	1.0	1.489	0.2211	0.1485
	1.5	1.483	0.3068	0.2069
	2.0	1.477	0.3861	0.2614
16.34	0.5	1.490	0.1510	0.1013
	1.0	1.481	0.2789	0.1883
	1.5	1.471	0.3876	0.2635
	2.0	1.461	0.4836	0.3310

, and the mean vertical aperture and standard deviation for a zero scenario are 1.5 mm and 0 mm, respectively. As shear displacement increased, the mean vertical apertures showed little change, while the standard deviation increased. The increase in JRC also resulted in an increase in standard deviation. The dimensionless variation coefficient ${\sigma }_s/e_m$ contained information on the channel size and the distribution level of oscillation.

2.2. Numerical Method

The flow of the isothermal and incompressible Newtonian fluid in a single fracture is governed by the Navier–Stokes equations, expressed as:

$$\frac{\partial \rho }{\partial t}+\frac{\partial \left(\rho v_i\right)}{\partial x_i}=0$$

(3a)

$$\frac{\partial \left(\rho v_i\right)}{\partial t}+\frac{\partial \left(\rho v_i{v}_j\right)}{\partial x_i}=-\frac{\partial p}{\partial x_i}+\frac{\partial }{\partial x_i}\left[\left(\mu +{\mu }_t\right)\left(\frac{\partial v_i}{\partial x_j}+\frac{\partial v_j}{\partial x_i}\right)\right]$$

(3b)

where $t$ is time, $v_i$ and $v_j$ are the velocity components, $x_i$ and $x_j$ are the coordinates, $p$ is the pressure, $\mu$ is the viscous coefficient, and ${\mu }_t$ is the turbulent viscous coefficient. This set of equations can be solved numerically using the commercial CFD software Fluent 6.3. A realizable k-ε model, which has shown good prediction results for phenomena such as rotation, boundary layer with significant backflow, and separation [29], was used to solve the flow problem. The turbulent kinetic energy and dissipation rate equations of the realizable k-ε model are as follows [29]:

$$\frac{\partial \left(\rho k\right)}{\partial t}+\frac{\partial \left(\rho k{v}_i\right)}{\partial x_i}=\frac{\partial }{\partial x_i}\left[(\mu +\frac{{\mu }_t}{{\sigma }_k})\frac{\partial k}{\partial x_j}\right]+G_k-\rho \epsilon$$

(4a)

$$\frac{\partial \left(\rho \epsilon \right)}{\partial t}+\frac{\partial \left(\rho \epsilon v_i\right)}{\partial x_i}=\frac{\partial }{\partial x_i}\left[\left(u+\frac{u_t}{{\sigma }_{\epsilon }}\right)\frac{\partial \epsilon }{\partial x_j}\right]+\rho C_{1}E\epsilon -\rho C_2\frac{{\epsilon }^2}{k+\sqrt{\nu \epsilon }}$$

(4b)

$$C_1=\mathit{max}\left(0.43,\frac{\eta }{\eta +5}\right)$$

(4c)

$$\eta ={\left(2E_{ij}{E}_{ij}\right)}^{1/2}\frac{k}{\epsilon }$$

(4d)

$$E_{ij}=\frac{1}{2}\left(\frac{\partial v_i}{\partial x_j}+\frac{\partial v_j}{\partial x_i}\right)$$

(4e)

where $k$ is the turbulent kinetic energy, $\epsilon$ is the dissipation rate, $G_k$ is the turbulent energy generated by the average velocity gradient, $\nu$ is the kinematic viscosity, $E$ is the mean strain rate, ${\sigma }_k$ is the turbulent Prandtl number of turbulent kinetic energy, and ${\sigma }_{\epsilon }$ is the turbulent Prandtl number of the dissipation rate. In Fluent, ${\sigma }_k$, ${\sigma }_{\epsilon }$, and $C_2$ were set as default constants, where ${\sigma }_k=1.0$, ${\sigma }_{\epsilon }=1.2$ and $C_2=1.9$ [30].

The momentum turbulent kinetic energy and turbulent dissipation rate were solved using second-order upwind schemes, and pressure–velocity coupling was performed using an embedded simple algorithm. The upper and lower fracture surfaces were set as walls with no slip and no flow. A constant velocity inlet with a uniform velocity profile was applied on the left boundary. The opposite end of the channel was modeled as a constant pressure outlet under atmospheric conditions. The density of the water was 998.2 kg/m³, and the dynamic viscosity coefficient was 1.003 × 10³ Pa·s (at 20 °C). The water flow in each fracture was modeled at five different shear displacements (0, 0.5, 1.0, 1.5, and 2.0 mm) and nine different injection velocities (0.001, 0.005, 0.01, 0.05, 0.1, 0.25, 0.5, 0.75, and 1.0 m/s). A total of 180 simulations were performed.

3. Results and Discussion

3.1. Emergence of Eddies

The Reynolds number ($Re$) provides a ratio of inertial forces to viscous forces. It is defined as follows for flow in fractures [7]:

$$Re=\frac{\rho v{e}_m}{\mu }$$

(5)

where $e_m$ is the characteristic dimension of the flow system that is equal to the mean aperture of the fracture, and $v$ is the characteristic flow velocity that is approximately equal to the inlet velocity.

Figure 7 shows the x-velocity distribution (x-velocity is the velocity component along the flow direction) and local eddy development at x = [75 mm, 85 mm] with inlet velocities of 0.001, 0.01, 0.1, and 1.0 m/s, corresponding to the Reynolds numbers of 1.5, 15, 150, and 1500, respectively. For a given fracture, eddy flow occurred with an increase in the Reynolds number due to the local effects of sharp-cornered asperities, which enhanced the tortuosity by changing local flow direction. The locations and sizes of the eddies were identified by checking the direction of streamlines. The fluid flowed slowly and backward near the margin of the local wider zones due to the inertial effect of fluid flow. Eddies were developed in local narrow-to-wide regions in the fracture channel, which is consistent with previous experimental observations [12]. When the Reynolds number was small, the pattern of streamlines indicated that the overall fluid flow was close to laminar flow in the fracture channel. With increasing Reynolds numbers, the size, number, and shape of eddies increased dramatically.

The temporal and spatial distribution of these eddy flow regions also changed with different JRCs and shear displacements. Figure 8 and Figure 9 describe the local eddy development with increasing JRCs and shear displacements, respectively. Because increasing the JRC and shear displacement drove the fracture to a high level of mismatch, the local regions changed from a relatively narrow channel to a wide region and thus produced more and larger eddies. In fractures with larger JRCs and shear displacements, eddy flow phenomena increased and the overall fluid flow became nonlinear, with more complex velocity distribution. The further occurrence and expansion of eddies significantly narrowed the effective flow channel, reduced the fracture permeability, and changed the fracture flow from linear to nonlinear.

3.2. Correlation between Inlet Flow Velocity and Hydraulic Gradient

The correlation of the average flow velocity to the hydraulic gradient between the inlet and outlet boundaries is important to describe the macroscopic fluid flow behavior. Figure 10 shows the fitting curve of the hydraulic gradient based on Equation (1) as a function of average flow velocity. The hydraulic gradients varied nonlinearly with average flow velocity, and the relationship was represented using the Forchheimer equation with a correlation coefficient greater than 0.999. A higher JRC or shear displacement resulted in a greater degree of mismatching and stronger nonlinearity. When the JRC value was small, there was less effect of shearing on the non-linearity of the pressure drop. It should be noted that the term “nonlinearity” of flow does not mean that the flow regime was turbulent. The small difference between the curves of JRC = 13.51 and JRC = 16.34 at large shear displacement may have been due to the non-stationary nature of fracture profiles with large JRCs, where the vertical aperture did not strictly obey the Gaussian distribution [5].

3.3. Forchheimer Coefficients

The Forchheimer equation (Equation (1)) proved to be adequate for the macroscopic description of water flow in rock fractures. The coefficients A and B are commonly written as:

$$A=\frac{\mu }{k_v}$$

(6a)

$$B=\beta \rho$$

(6b)

where $\beta$ is the non-Darcian coefficient or inertial resistance coefficient, $k_v$ is the intrinsic permeability or viscous permeability, $\rho$ is the fluid density, and $\mu$ is the fluid dynamic viscosity.

Based on the best correlation between the flow flux and the hydraulic gradient using the Forchheimer equation, the coefficients $\beta$ and $k_v$ were calculated. Figure 11 and Figure 12 show the correlation of $\beta$ and $k_v$, respectively, to shear displacement with different JRCs. With an increase of shear displacement, the linear coefficient $k_v$ decreased but the nonlinear coefficient $\beta$ increased, further confirming the presence of nonlinear pressure drops. The greater the value of the JRC, the stronger the effect of change of shear displacement on $\beta$.

The emergence of nonlinear flow in rough-walled fractures during the shear process was strongly related to the degree of mismatch in the two fracture surfaces [8] and can be characterized as a dimensionless variation coefficient of mean aperture ${\sigma }_s/e_m$. Because the non-Darcian coefficient $\beta$ was measured using L⁻¹ and mean aperture measured $e_m$ using L, the following equation was used to measure the variation of $\beta$:

$$\beta e_m=JR{C}^{{\varsigma }_1}\left[{\varsigma }_2+{\varsigma }_3{\left(\frac{{\sigma }_s}{e_m}\right)}^{{\varsigma }_4}\right]$$

(7)

where ${\varsigma }_1$, ${\varsigma }_2$, ${\varsigma }_3$, and ${\varsigma }_4$ are dimensionless parameters that can be determined by the numerical data for $\beta$, as shown in Figure 12. The most accurate values of the four coefficients were obtained using the Levenberg-Marquardt (LM) optimization algorithm as follows: ${\varsigma }_1=0.7123$, ${\varsigma }_2=4.0882$, ${\varsigma }_3=52.2052$, and ${\varsigma }_4=1.8087$. Figure 13 shows a comparison between the prediction using Equation (7) and the numerical data, indicating a reasonable agreement. The correlation coefficient of 0.944 is high as the four coefficients are globally optimized for the data of four JRCs. The equation for $\beta$ can quantify the inertial effect induced by the fracture roughness and spatial heterogeneities of aperture, and provide some insight into the non-Darcian flow behaviors in rough-walled fractures.

3.4. Critical Reynolds Number

The critical Reynolds number is defined as the critical state in which the flow changes to a non-Darcian zone, where the inertial pressure drop reaches $\alpha$ percentage of the overall pressure drop. $\alpha$ is given as:

$$\alpha =\frac{B{v}^2}{Av+B{v}^2}=\frac{1}{1+\frac{\mu }{\beta \rho kv}}$$

(8)

The critical Reynolds number can be determined by choosing a proper value of $\alpha$. Considering that $\alpha$ is usually taken as 10% [7,31], the critical Reynolds number can be written as:

$${\mathit{Re}}_c=\frac{e_m}{9\beta k_v}$$

(9)

Equation (9) indicates that the transit of fracture flow from linear to nonlinear is a combined effect of aperture, viscous permeability, and inertial resistance. The critical Reynolds number as a function of the variation coefficient ${\sigma }_s/e_m$ and JRCs is plotted in Figure 14. The fluid flow is easier to change into the nonlinear flow regime with a smaller critical Reynolds number. With an increase in the variation coefficient ${\sigma }_s/e_m$, the critical Reynolds number decreased at first and then increased. The fracture with higher JRCs led to the enhancement of flow nonlinearity and earlier onset of nonlinear laminar flow during shearing. Due to changes in the variation coefficients governed by the shear process and JRCs, significant variations in the coefficients $k_v$ and $\beta$ and the critical Reynolds number were observed. It is clear that surface roughness and mismatching played important roles in the variation of the critical Reynolds numbers and the nonlinearity of the fluid flow through rough-walled rock fracture. It should be noted that the increase in the critical Reynolds number in accordance with increasing heterogeneities of aperture may be due to a significant reduction in fracture conductivity [15].

4. Conclusions

We performed a series of numerical simulations on four standard fractures with different shear displacements to examine the influence of shearing and fracture roughness on non-Darcian flow characteristics. The results showed that increasing the JRC and shear displacement resulted in a high degree of mismatch in the fracture, and eddy flow occurred and expanded, causing a reduction in fracture conductivity and an enhancement of the local complexity of nonlinear fracture water flow. The relationship between the pressure gradient and average flow velocity revealed that the Forchheimer equation offers a good description for the non-Darcian flow behavior in rough-walled fractures with shearing. The effect of shearing and the JRC on the Forchheimer coefficients was examined. A three-parameter equation as a function of the JRC and the variation coefficient of aperture distribution ${\sigma }_s/e_m$ was developed for empirical parameterization of the non-Darcian coefficient $\beta$. The equation has the advantage of dimensional consistency and fits well with the numerical results. Based on the Forchheimer equation and taking 10% of the nonlinear effect as a critical state to distinguish between linear and nonlinear flow, the critical Reynolds number for the specific flow cases was estimated by the combined effect of aperture, viscous permeability, and inertial resistance at the macroscopic level.