Numerical simulation of two-dimensional cyclic water injection in rock fractures

doi:10.21203/rs.3.rs-3336766/v1

Download PDF

Research Article

Numerical simulation of two-dimensional cyclic water injection in rock fractures

https://doi.org/10.21203/rs.3.rs-3336766/v1

This work is licensed under a CC BY 4.0 License

Version 1

posted

You are reading this latest preprint version

Understanding the shear-seepage behaviors in rock fractures is crucial in excavation of tunnels through fractured surrounding rocks and mining explorations, however, the cyclic injection mechanism of fractures is not yet fully understood. This study utilized the COMSOL Multiphysics software to simulate the cyclic water injection in rock fractures with various injection amplitudes and frequencies. We documented that when subjected to cyclic water injection, the fluid flow in the rock fractures exhibited periodic fluctuations. As water injection amplitude increases, the flow pattern within the low-flow region become more turbulent, creating significant disturbances to the overall flow pattern. Additionally, the permeability of rock fractures varies periodically with the inlet flow rate, resulting in a permeability-flow rate relationship that forms an elliptical closed curve over time. Moreover, the amplitude of fluctuation in permeability increases as the water injection frequency and amplitude increase, illustrating that water injection within a specific range could enhance permeability in rock fractures. Furthermore, we proposed an approach to identify and quantify the primary channel and eddy flow areas.

Fracture flow

Numerical simulation

Cyclic injection

Permeability

Main watershed

The seepage law of two-dimensional flow channel was explored by combining real rough rock fractures modeling and numerical simulation of fluid dynamics.
The variations of permeability and flow rate of the fracture surface of constant water injection, variable frequency injection and variable amplitude water injection were systematically studied.

A simple and effective method that can be used to distinguish between the main watershed and the eddy region and quantitatively calculate the size of various regions is proposed.

With the vigorous advancement of geothermal energy exploitation, shale gas development, large-scale water conservancy hub projects, deep tunnels, and abandoned nuclear fuel underground storage projects, understanding the hydraulic characteristics of fractured rock masses is significant for the safe construction and operation of these projects (Falcone et al. 2018; Zhang and Zhao 2020; Alexander et al. 2011; González 2016; Rutqvist and Stephansson 2003; Leung and Zimmerman 2012). The connected fracture networks are widely distributed in rock having rough surface morphologies, which the connected fractures are served as the primary channels for fluid permeation (Liu et al. 2020; Huang et al. 2019). Although many scholars have proposed formulas for predicting the hydraulic characteristics of rock fractures considering multiple factors, a unified conclusion has not yet been reached due to the complex and stochastic seepage conditions (Zhao et al. 2021; Ma et al. 2013). For example, within eddy area, the correlation between flow rate and pressure drop is nonlinear and cannot be described using the commonly applied cubic law (Kohl et al. 1997). Additionally, the opening and closing of fractures resulting from normal or shear stress will significantly affect the nonlinearity (Zou et al. 2017; Li et al. 2019). Although some scholars have carried out laboratory tests on shear seepage of rock fracture (Xia et al. 2020; Esaki et al. 1999; Rong et al. 2016), there were few studies on corresponding numerical simulation, and little numerical simulations demonstrated the variations of water flow channels. Moreover, the impact of the amplitude and frequency of water injection on the hydraulic properties remains unclear, highlighting the need for further research to understand the fluid flow characteristics within rock fractures.

Currently, researches on rock fractures seepage are primarily being conducted through laboratory tests and numerical simulations both domestically and internationally. Scholars have conducted constant water injection experiments on rock fractures seepage using rock direct shear seepage apparatus (Ji et al. 2020; Hofmann et al. 2018; Xia et al. 2020; Xiao et al. 2013), mainly focusing on researching the impact of normal and shear loads during the shearing process on fracture seepage characteristics (Liu et al. 2020; Zou et al. 2017). Esaki et al. (1999) conducted coupled shear-seepage tests on granite under different constant normal loads and found that the variations of granite fractures permeability and shear displacement were similar under small shear displacement. Chen et al. (2021) utilized a newly developed experimental apparatus coupling seepage-stress and observed that shearing can alter the asperity of fractures, leading to variations in the permeability magnitude. Rong et al. (2016) conducted indoor shear-seepage experiments on rock joints with various roughness levels. The findings indicated that the deformation resulting from shearing had a significant impact on the linear coefficient “a” and nonlinear coefficient “b” in Forchheimer’s law. With the contraction and expansion of the fracture aperture due to shearing, the critical Reynolds number (Rec) exhibits a pattern of decreasing initially before increasing. Passelegue et al. (2018) conducted shear-seepage tests on granite under triaxial stress conditions. The results showed that the peak fluid pressure which caused fault activation was related to the fluid injection rate, and the higher the injection rate, the higher the peak fluid pressure that caused fault activation. Liu et al. (2019) conducted a digital speckle pattern observation test on the initiation and propagation of fractures in rock mass. By analyzing the evolution of the shale displacement field and strain field during hydraulic fracturing propagation, the dynamic expansion characteristics of rock fracture were revealed. Ji et al. (2021) investigated the slip characteristics of natural granite fractures by using constant water injection and triangular waveform cyclic water injection. These findings show that the peak slip rate of rock fractures during cyclic water injection was lower than that during constant water injection. Zhuang et al. (2020) adopted six different water injection fracturing schemes with a true tri-axial test device and evaluated their effects on the hydraulic characteristics of rock fractures. They found that cyclic pulse pressurization (CPP) is superior to the other five schemes and can improve hydraulic properties and reduce the seismic activities. In the aforementioned experiments, the effects of water injection amplitude and frequency on the seepage regime are not mentioned/analyzed.

Thanks to advances in computational fluid dynamics (CFD), numerous scholars have conducted precise numerical simulations of seepage in rock fractures, both at the two-dimensional and three-dimensional scales (Zhong et al. 2021; Huang et al. 2021; Cardenas et al. 2009; Chaudhary et al. 2013; Xing et al. 2021). Notably, two-dimensional fracture models can serve as a reference for three-dimensional fracture models analysis (Lee et al. 1996). Dang et al. (2019) implemented a numerical simulation method that incorporates FLAC3D-Fluent, allowing for the analysis of the deformation effect of fracture surfaces and the assessment of how hydraulic apertures impacts the permeability of rock fractures. Zou et al. (2015) introduced a theory of wavelet decomposition for roughness. Through simulations of fluid flow in fractures using various models of different roughness, they discovered that joint roughness has a profound impact on the generation and shifting characteristics of eddy currents, which are more likely to emerge at the edge of the fracture (Quinn et al. 2010). Xie et al. (2015) used the COMSOL software to solve the N-S equations and performed computational analysis of the single fracture seepage characteristics at various shear displacement. Their study yielded insight into the variance of equivalent hydraulic aperture, mechanical aperture, and volume flow rate with shear displacement. Moreover, they examined the distribution and flow characteristics of fluid volume within the fracture under different shear directions. Dong et al. (2022) utilized the CDP model in ABAQUS to simulate the damage incurred to the surrounding rock mass of boreholes under pulse water pressure. Their research indicated a cumulative damage phenomenon of the perforation surface under pulse cyclic loads, and revealed that cracking due to pulse cyclic loading could reduce the initial fracturing pressure. Sun et al. (2014) performed a three-dimensional numerical simulation on the hydraulic fracturing of layered rock masses, identifying the distribution of stress fields, pore pressures, damage to fracture surfaces, aperture, and the temporal-spatial distribution characteristics of fracture propagation. Although they have conducted in-depth and detailed numerical simulation studies from different dimensions and different factors affecting shear seepage, few numerical simulation studies have been conducted on dynamic water injection.

We generated a detailed 3D model of an actual rock fracture, followed by selecting a two-dimensional flow channel for numerical study, besides, monotonic water injection and cyclic water injection, employing varied amplitudes and frequencies were conducted on a rock fracture model. The established numerical model, employed through the COMSOL software, includes the N-S equation for momentum conservation and the continuity equation for mass conservation as a finite element method. These simulations assessed the impact of cyclic water injection frequency and amplitude on the flow and permeability of rock fractures, providing insight into the seepage laws and evolutionary characteristics in engineering practice, such as tunnel excavation in water-rich areas and geothermal exploitations.

2.2 Fracture seepage equation

(1) Fluid control equation

The N-S equation reflects the general law of fluid flow in rough fractures. For isothermal, steady-state, and incompressible single-phase Newtonian fluid flows, the N-S equation can be expressed as:

$$\rho \left( {{\mathbf{u}} \cdot \nabla } \right){\mathbf{u}}=\mu {\nabla ^2}{\mathbf{u}} - \nabla P$$

According to the law of mass conservation, the continuity equation is:

$$\nabla \cdot {\mathbf{u}}=0$$

Where in the above formula, u = [u, v, w] represents the velocity vector of fluid in three directions; P is the total water pressure; ρ is the density of the fluid in the fractures; µ is the fluid dynamic viscosity coefficient.

When the fracture aperture is small and the Reynolds number is very small (Re≪1), the fluid inertial force will be small compared with the viscous force, so this inertial force can be ignored when solving the N-S equation, that is, the inertia term in the N-S equation can be ignored, and the N-S equation can be converted into the Stokes equation.

$$\mu {\nabla ^2}{\mathbf{u}}=\nabla P$$

Given the two primary issues in fracture seepage research - namely, shear seepage and variable frequency seepage - different control equations are utilized. Shear seepage pertains to the steady-state problem and is thus independent of time. Meanwhile, variable frequency seepage pertains to a transient problem, dependents on time, and features a control equation that incorporates a time term.

The steady-state control equation for the laminar flow module is:

$$\rho \left( {{\mathbf{u}} \cdot \nabla } \right){\mathbf{u}}=\nabla \cdot \left[ { - p{\mathbf{I}}+{\mathbf{K}}} \right]+{\mathbf{F}}$$

The transient control equation for the laminar flow module is:

$$\rho \frac{{\partial {\mathbf{u}}}}{{\partial t}}+\rho \left( {{\mathbf{u}} \cdot \nabla } \right){\mathbf{u}}=\nabla \cdot \left[ { - p{\mathbf{I}}+{\mathbf{K}}} \right]+{\mathbf{F}}$$

$${\mathbf{K}}=\mu \left[ {\nabla {\mathbf{u}}+{{\left( {\nabla {\mathbf{u}}} \right)}^T}} \right]$$

$$\rho \nabla \cdot {\mathbf{u}}=0$$

where Eq. (7) is the fluid continuity equation. where ρ is the fluid density; p is the fracture fluid pressure; µ is the fluid dynamic viscosity coefficient; t is time; u is the flow rate vector; F is the fluid volumetric force; I is a unity matrix; p and u are the dependent variables in the control equation; ▽ is the Hamiltonian operator (differential operator). From the above control equations, the N-S equation is the core of the numerical simulation of variable frequency seepage.

(2) Reynolds number

The Reynolds number Re, defined as the ratio of the fluid inertial force to the viscous force, is calculated in the fracture as follows:

$${R_e}=\frac{{pQ}}{{\mu w}}$$

$$Q={A_h}u=w{e_h}u$$

Where Q is the volume flow; is the fracture width, take 1 in the two-dimensional model.

(3) Permeability

At lower flow velocities, rock fracture seepage obeys Darcy's law, which is:

$$Q=\frac{{k{A_h}}}{\mu }\nabla P$$

where Q is the volume flow, k is the fracture permeability, ▽P is the pressure gradient at both ends of the joint, A_h is the seepage cross-sectional area, and µ is the dynamic viscosity coefficient of the fluid. Darcy's law states that the permeable flow rate in a medium is proportional to the head difference and the cross-sectional area of the overflow, and inversely proportional to the length L of the seepage path, so the Darcy permeability k can be expressed as:

$$k=\frac{{\mu QL}}{{{A_h}\Delta P}}=\frac{{\mu QL}}{{w{e_h}\Delta P}}$$

where △P is the pressure difference between the two ends of the joint, L is the length of the seepage path, and e_h is the equivalent hydraulic opening. The relationship between permeability K and permeability k is:

$$K=\frac{{k{\gamma _w}}}{\mu }$$

To calculate the permeability of a rock joint, the equivalent hydraulic opening must be determined using the cube law. Typically, fluid flow through a fracture can be simplified as occurring between two smooth parallel plates with zero fracture roughness coefficient (JRC = 0). With Darcy’s law assumed, the cube law is derived using this parallel plate model. Specifically, the volume flow rate Q is proportional to the cube of the ideal aperture size e. This leads to the cubic law formula for seepage through rough joints, as follows:

$$Q= - \frac{{w{e_h}^{3}}}{{12\mu }}\nabla P$$

The relationship between permeability and equivalent hydraulic aperture obtained by combining Darcy's law is:

$$k=\frac{{e{h^2}}}{{12}}$$

According to equations (11) and (13), the fracture permeability k can be solved.

3.1 Modeling process

(1) Preparation of rock specimens: The test utilized a complete basalt specimen with dimensions of 130 mm (L) × 130 mm (W) × 130 mm (H). The specimen was split into two halves in the middle to create two rough fracture surfaces, as illustrated in Fig. 1a. The space between the upper and lower rough surfaces forms channel for fluid flow.

(2) 3D optical scanning: 3D scanning of rock fracture surface to obtain spatial point data of upper and lower rough fracture surface, which is used for the establishment of channel geometric model. As shown in Fig. 1b, the two cameras of the 3D scanner collect asperities spatial information, and then the computer calculates the 3D coordinates of the target point. Then, a complete 3D rough fractures numerical model is established.

(3) Two-dimensional modeling: The fracture plane coordinates obtained from 3D scanning were converted into ASCII files and imported into FLAC3D to reconstruct the three-dimensional rough fractures, visualized in Fig. 1c. The upper rock was horizontally sheared 5.0 mm while the lower rock is fixed, since the fracture surface is rough, the upper surface will appear dilation and compaction during shearing. The three-dimensional numerical model was processed in FLAC3D, and a two-dimensional section perpendicular to the rock fracture surface was obtained, as shown in Fig. 1d. The model’s left and right boundaries were the fluid inlet and outlet, respectively, while the upper and lower boundaries consisted of the rock side walls. The hydraulic head was applied at the inlet, while a free surface boundary condition was applied at the outlet. The model included an inlet opening of 1.42 mm, an outlet opening of 0.99 mm, and a fluid channel length of 21.01 mm.

(4) Mesh generation: This paper selects a free triangle grid to divide the fluid. In addition, according to the flow characteristics of the boundary layer of the rock wall, a quadrilateral grid is used for the upper and lower walls, and the fluid flow along the x direction. as shown in Fig. 2a and Fig. 2c.

3.2 Initial parameters and boundary conditions

The following simulation was conducted in the COMSOL software. The density of the fluid is 1.0 g/cm³, the dynamic viscosity coefficient is 0.001 Pa·s, and the fluid is set to incompressible at room temperature (25º), ignoring thermal and gravitational effects. The initial flow rate of the fluid u₀ = 0 m/s, the initial pressure is standard atmospheric pressure, and there is no slippage on the side wall. The fluid inlet is located on the left side of the fracture flow channel and is controlled by the average flow Q, which fluctuates with time, and the control equation is shown in the formula (15):

$$Q={Q_0}+A\sin \left( {\omega t+\theta } \right)$$

where Q is the inlet boundary flow, Q₀ is the initial flow, A is the flow amplitude, ω is the angular frequency, t is the time, and θ is the phase. Set the initial flow rate Q₀ to 100 ml/s, the initial phase θ₀ = 0, and the conversion relationship between angular frequency ω and frequency f is shown in Eq. (16):

$$\omega =2\pi f$$

The outlet is located on the right side of the flow channel, and the outlet is free to flow set as the standard atmospheric pressure.

3.3 Frequency conversion seepage calculation process

The flow rate in the fracture channel model is regulated by a sinusoidal function (Eq. 15). The evolution of flow rate and pressure under three different boundary conditions, consisting of constant water injection, variable amplitude injection, and variable frequency water injection, are obtained and presented in Table 1. Group 1 is constant water injection, keeping the inlet flow rate of 100 ml/s; The frequency of cyclic water injection in Group 2 remained constant at 1 Hz, and the amplitudes were 10 ml/s, 20 ml/s, and 50 ml/s, respectively. The water injection amplitude in Group 3 was set to 20 ml/s, and the frequencies were 0.05 Hz, 0.1 Hz, 0.25 Hz, 0.5 Hz, 0.75 Hz, and 1 Hz, respectively. In the transient laminar flow calculation of COMSOL, the variation in fluid velocity and fluid pressure of a single flow channel over a period of 5 cycles is solved.

In the model, the inlet boundary can be specified by flow rate, pressure, or mass flow. Through multiple tests, the method of controlling the inlet flow and outlet pressure was found to offer the best convergence effect. The calculation was carried out for a solution range of 5 cycles.

Table 1

Numerical simulation scheme of periodic water injection.
Group	Inlet boundary			Outlet boundary	Solving time
	Flow rate (ml/s)	Frequency (Hz)	Period(s)	Relative pressure (Pa)	Total time (s)
Constant injection	100	—	—	0	22 min 38s
Variable amplitude injection	100 ± 10	1	1	0	25 min 24s
	100 ± 20	1	1	0	25 min 55s
	100 ± 30	1	1	0	26 min 12s
Variable frequency injection	100 ± 20	0.05	20	0	24 min 44s
	100 ± 20	0.1	10	0	24 min 47s
	100 ± 20	0.25	4	0	25 min 16s
	100 ± 20	0.5	2	0	25 min 21s
	100 ± 20	0.75	1.33	0	25 min 49s
	100 ± 20	1	1	0	25 min 55s

The variable amplitude injection scheme adopts sinusoidal waveform, the initial flow rate of water injection is 100 ml/s, and finally the periodic water injection flow curves under different amplitudes are shown in Fig. 3. The variable frequency injection scheme also adopts a sinusoidal waveform with a constant amplitude of 20 ml/s. Set six different water injection frequencies of 0.05 Hz, 0.1 Hz, 0.25 Hz, 0.5 Hz, 0.75 Hz, and 1 Hz, as shown in Fig. 4.

The analysis of the numerical simulation results for frequency conversion seepage flow in fracture focuses primarily on the qualitative analysis of flow rate distribution characteristics. In this study, the entry flow line is set to 30. Based on the observed streamline distribution characteristics, the seepage characteristics of the fracture can be understood to a certain extent. The expression of the streamline is presented as follows:

$$u=\sqrt {u_{x}^{2}+u_{y}^{2}}$$

In Eq. (17), u_x and u_y represents the flow rate in the x direction and the y direction, respectively.

4.1 Constant water injection scheme

Constant water injection refers to the relative seepage pressure (0 Pa) and fluid rate (100 ml/s) remain all constant. As shown in Fig. 5, upon observing the distribution of streamlines, we observed that the narrow watershed appears darker, indicating a higher flow rate, the fluid does not seem to follow the shortest path through the fracture due to the joint effect of side wall obstruction and a relative larger fluid rate, resulting in a unique flow trajectory. The upper fracture was observed to have a high flow rate resulting in a relatively smooth streamline, whereas the flow rate at the bottom was slower, measuring less than 0.06 m/s. Near the low flow rate area’s edge wall, a clear eddy was observed, the emergence of eddy current will lead to the loss of kinetic energy and reduce the permeability of fractures. Further examination of the flow line showed that the fluid flowed clockwise in this eddy area, with a notable dividing line between it and the main flow channel.

4.2 Variable amplitude water injection scheme

We select the last cycle, that is, the fifth cycle, for analysis to avoid interference caused by aperiodic factors. The initial flow value, flow peak, and flow valley values are the three characteristic time points of t = 4 s, t = 4.25 s, and t = 4.75 s of the fifth cycle, as shown in Fig. 6.

Figure 6 shows the variations in fluid flow rates within the fractures. The fluid appears to flow more rapidly in the narrow area of the fracture and the outlet, while slower flow rates are observed at the edge of the flow channel due to the boundary effects. While observing the flow rate across various time points, the flow rate inside the fracture increases from its initial value (t = 4 s) to its peak value (t = 4.25 s), allowing for the expansion of the dark high flow rate area. However, when the boundary flow decreases from its peak value (t = 4.25 s) to its valley (t = 4.75 s), the flow rate in the fracture decreases as a whole. Furthermore, the eddy region expands, and small eddies are separated from the larger eddies.

The comparisons between different water injection amplitudes indicates that, as the amplitude increases from 10 ml/s to 30 ml/s, there is an overall increase in flow rate within the fracture corresponding to the flow peak, and a decrease in the flow rate corresponding to the flow valley. Furthermore, an increase in amplitude also causes a change in flow pattern within the eddy region. The evolution of flow rate over time is found to be closely related to changes in boundary flow rates, while the influence of amplitude changes on flow rate distribution is mostly reflected in the eddy region.

Using Darcy’s law and Cube’s law, the law of fracture permeability over time is calculated based on Eq. (11), and the results are shown in Fig. 7. We observe that, starting from the second cycle, when the inlet boundary is constant water injection, permeability remains relatively constant. However, when the inlet boundary adopts variable frequency water injection, the permeability begins to fluctuate over time, showing stable and regular fluctuation from the second cycle. Observing the permeability curve at various water injection amplitudes, it manifests that the peak permeability tends to increase as the water injection amplitude increases while the valley permeability decreases as the water injection amplitude increases. Notably, the increase in water injection amplitude has a greater impact on peak permeability than on the valley permeability. As a result, the permeability amplitude increases with increasing water injection amplitude. The phenomenon is reasonable because larger amplitude, generating larger dilation and aperture volume, resulting in the increment of permeability.

Figure 8 displays a clear nonlinear relationship between inlet boundary flow and fracture permeability. The abscissa of the figure represents the inlet boundary flow, while the ordinate represents the fracture permeability. The circular pattern in the figure reflects the process of boundary flow rising and falling, and each anticlockwise cycle, in the direction of the arrow, forms a complete water injection cycle. Using the curve with an amplitude of 100 ± 30 ml/s in Fig. 9 as an example, the relationship between inlet boundary flow and fracture permeability is analyzed. Starting from a flow rate of 70 ml/s, the inlet boundary flow increases in the direction indicated by the red arrow, resulting in a decrease in permeability within the red area. As the flow rate increases to 126 ml/s, the permeability starts to increase again. When the flow rate reaches its peak at 130 ml/s and enters the blue area, the flow rate begins to decrease, while the permeability continues to increase. Ultimately, when the flow rate drops to 74 ml/s, the permeability reaches a peak of 12×10^− 8 m² before moving back into the red zone to complete a water injection cycle.

4.3 Variable frequency injection scheme

Under the conditions of six sets of periodic cyclic water injection at different frequencies, the evolution of flow rate with time is illustrated in Fig. 10.

The flow pattern of the fluid changes periodically over time, with the analysis progressing from left to right in time. As the flow rate increases from its initial value to its peak value, the flow rate inside the fracture increases concurrently with larger boundary flow. This is manifested through the expansion of the dark high flow rate area. Conversely, when the boundary flow decreases from the peak to the valley, the flow rate in the fracture decreases as a whole, resulting in the expansion of the eddy area. This phenomenon becomes more pronounced at higher frequencies and follows a similar pattern to that of variable amplitude water injection.

From the comparison between different frequencies, the water injection frequency increases from 0.05 Hz to 1 Hz, the flow rate of water injection changes more severely, the perturbation of the fracture flow is greater, the overall flow rate inside the fracture corresponding to the flow peak is larger, the flow rate corresponding to the flow valley value is smaller, and the increase of frequency causes the change of the flow pattern in the eddy area, the higher the frequency, the greater the flow rate, but the influence of frequency on the flow pattern in the fracture is smaller than the amplitude.

Four typical variable frequency water injection schemes (f = 0.05 Hz, f = 0.5 Hz, f = 0.75 Hz and f = 1 Hz) are selected to plot the permeability curve with time, the fracture permeability changes periodically. Comparing Fig. 11 with Fig. 4, the fluctuation of permeability and boundary flow does not rise or fall synchronously. Taking f = 1 Hz as an example, as shown in Fig. 4f, the flow curve rises from 0 s to 0.25 s, but the permeability curve in Fig. 10d gradually decreases. Nevertheless, when the flow rate decreased from a peak of 0.25 s to a valley of 0.75 s, the permeability increased from a valley of 8.9×10^− 8 m² to 11.1×10^− 8 m².

Figure 11 displays the circulation curve of permeability versus boundary flow, which expands longitudinally with the increase of frequency. For instance, the curve of 1 Hz, the inlet boundary flow increases along the direction indicated by the red arrow, starting from a flow rate of 80 ml/s. As the flow rate increases to 118 ml/s, permeability decreases within the red area. Nevertheless, when the flow rate reaches its peak at 120 ml/s and enters the blue area, the permeability continues to increase even as the flow rate decreases. Ultimately, when the flow rate drops to 84 ml/s, the permeability reaches a peak of 11.2×10^− 8 m², before moving back into the red zone to complete a water injection cycle.

Combining the analysis of Fig. 10 and Fig. 11, we observeed that the flow rate near the boundary does not change under constant water injection, resulting in a corresponding constant permeability shown as a point on the figure. However, under variable frequency water injection, the permeability of rock fractures changes periodically with the inlet flow, resulting in a permeability-flow curve with time that resembles a closed elliptical shape. The lateral width of the elliptical curve represents the variation range of flow, while the longitudinal height represents the varied range of permeability. Accordingly, an increase in the inlet flow amplitude results in the elliptical curve expanding outward from the same center. This expansion reflects that the amplitude of permeability increases with the increase in flow amplitude, with the impact on peak permeability being more pronounced than on the valley. With the increase of inlet flow frequency, the elliptic curve expands along the longitudinal height, reflecting that the permeability amplitude increases with larger flow frequency, and the influence on the peak permeability is more obvious than the valley.

4.4 Mainstream domain of rock seepage

We propose a method of identifying the mainstream domain and the eddy area which quantitatively calculate the proportion of the mainstream in the flow region, and explore the relationship between permeability and the domain area.

(1) Mainstream domain identification method

Affected by the geometry boundary, the fluid in the fracture does not flow uniformly, but is concentrated in a certain area. The mainstream domain can be separated from the eddy area with low velocities. Zhou et al. (2019) designed an automatic identification method for the eddy zone in fracture flow, which is used to detect the boundary between the recirculation and the mainstream region.

This paper proposes a simpler and effective method for identifying the mainstream region. Specifically, the area within the flow channel where flow rate is greater than the average flow rate is defined as the mainstream region, while the area where flow rate is less than the average flow rate is defined as the eddy region. This method can divide the mainstream in real-time within COMSOL software and extract it directly during post-processing

Figure 12 illustrates the flow area before and after the identification. The boundary between the two kinds of regions is clear, with the blue area representing the mainstream and the white area representing the eddy. Within the mainstream region, the streamlines parallel each other, extending from the inlet to the outlet, with minimal resistances resulting in faster flow rate. In contrast, the eddy region’s streamline forms closed loops, circulating within the white area and not participating in the mainstream area’s flow. Consequently, the eddy has negative effect on the fracture permeability.

(2) Calculation and analysis of the proportion of mainstream region

Figure 13 illustrates the evolution of permeability and the mainstream proportion over time at various water injection frequencies. As demonstrated in Fig. 13a, the permeability is closely linked to the area of the mainstream domain, with an increase in permeability resulting in a corresponding increase in the mainstream domain area. Similarly, a decrease in permeability leads to a decrease in the mainstream domain area. As a result, both variables fluctuate periodically within a stable cycle.

Zhou et al. (2019) proposed that an increase in the eddy region during steady-state flow results in additional kinetic energy dissipation, leading to a reduction in permeability. this study adopted the transient flow analysis, moreover, we conclude that an increase in the eddy area leads to a decrease in the mainstream domain area, resulting in additional kinetic energy dissipation and a reduction in fluid passing capacity. Consequently, a decrease in permeability occurs, perhaps explaining the periodic fluctuations of permeability over time. Essentially, this process represents transient, nonlinear flow.

Figure 14a illustrates the periodic fluctuations of permeability over time, with the fifth period emphasized in Fig. 14b. A phase difference exists between the two curves within this stable period. This phase difference between permeability and the main watershed gradually increases as the water injection frequency increase, which appears the similar phenomenon between dynamic normal stress and shear stress in previous rock fracture shearing experiments (Dang et al. 2016; Dang et al. 2022). This phenomenon arises due to the intensification of nonlinear flow within the fracture as the water injection frequency rises. Specifically, changes in the eddy zone will lag behind changes in the permeability due to this intensification, resulting in a time difference between the two variables. Such deviation increases accordingly with the injection frequency.

In this paper, we made fluid simulation in COMSOL to investigate flow characteristics under constant water injection, cyclic flow rate injection with variable amplitude and frequency in a sheared rock fracture. The main findings of this paper are:

(1) Under variable water injection frequency, flow trajectories change periodically with the inlet flow. Greater inlet flow amplitudes result in more significant disturbances to the flow pattern and more complex eddy phenomena in regions of low flow rate. In contrast, the frequency of the inlet flow has a minor impact on the flow morphology when compared to the amplitude effects.

(2) The permeability of rock fractures varies periodically with inlet flow subjected to cyclic fluid injection, moreover, permeability-boundary flow curve is a closed ellipse over time. Larger fluid injection amplitude and frequency results in the permeability-boundary flow rate elliptical curve expands outward from the same center, reflecting that the magnitude of permeability change increases with larger fluid injection amplitude and frequency. The impact of fluid injection amplitude on permeability is greater, and the injection amplitude have a more pronounced effect on the peak permeability.

(3) An effective method is proposed for identifying the flow region that can distinguish the mainstream and eddy regions in rock joints and quantitatively calculate the size of each region. Our findings reveal a correlation between the area of the mainstream and permeability, with changes to mainstream area lagging behind changes to the permeability. The phase difference between these two variables increases with larger injection.

ρ	[kg/m³]	Fluid density
P	[MPa]	Total water pressure
μ	[-]	Dynamic viscosity coefficient of fluid
Re	[-]	Reynolds number
Q	[ml]	volume flow
k	[md]	Permeability
A_h	[m²]	Cross-sectional area of channel
L	[m]	The length of the seepage path
e_h	[-]	Equivalent hydraulic opening
K	[-]	Permeability coefficient
Q₀	[ml]	Initial volume flow
A	[-]	Flow amplitude
w	[mm]	Fracture root mean square heights
ω	[rad/s]	Angular frequency
t	[s]	Time
θ	[°]	Phase
θ₀	[°]	Initial phase

Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Acknowledgments

This research is financially supported by

Alexander T, Baihly J, Boyer C, et al (2011) Shale gas revolution[J]. Oilfield Rev 23(3): 40-57
Cardenas MB, Slottke DT, Ketcham RA, et al (2009) Effects of inertia and directionality on flow and transport in a rough asymmetric joint[J]. J Geophys Res Solid Earth 114(B6)
Chaudhary K, Cardenas MB, Deng W, et al (2013) Pore geometry effects on intrapore viscous to inertial flows and on effective hydraulic parameters[J]. Water Resour Res 49(2): 1149-1162
Chen W Z, Sui Q, Wang L Y, et al (2021) Development and application of coupled stress-seepage shear flow test system for rock joint surface [J]. Ch J Rock Mech Eng 2021, 40(9): 1729-1738. (in Chinese)
Dang W, Konietzky H, Frühwirt T (2016) Direct shear behavior of a plane joint under dynamic normal load (DNL) conditions[J]. Eng Geol 213: 133-141
Dang W, Tao K, Chen X (2022) Frictional behavior of planar and rough granite fractures subjected to normal load oscillations of different amplitudes[J]. J Rock Mech Geotech Eng 14(3): 746-756
Dang W, Wu W, Konietzky H, et al (2019) Effect of shear-induced aperture evolution on fluid flow in rock joints[J]. Comput Geotech 114: 103152
Dong KX, Zhao AN, Wang SL, et al (2022) Influence of pulse parameters on hydraulic fracturing initiation pressure of tight reservoir [J]. J Liaoning Tech Univ (in Chinese)
Esaki T, Du S, Mitani Y, et al (1999) Development of a shear-flow test apparatus and determination of coupled properties for a single rock joint[J]. Int. J. Rock Mech. Min. Sci 36(5): 641-650
Falcone G, Liu X, Okech RR, et al (2018) Assessment of deep geothermal energy exploitation methods: The need for novel single-well solutions[J]. Energy 2018, 160: 54-63
González MYV (2016) Shale gas: Opportunities and challenges between Mexico and the United States[J]. Nat Resour J 56(2): 313-328
Hofmann H, Zimmermann G, Zang A, et al (2018) Cyclic soft stimulation (CSS): a new fluid injection protocol and traffic light system to mitigate seismic risks of hydraulic stimulation treatments[J]. Geotherm Energy 6(1): 1-33.
Huang N, Jiang YJ, Cheng YF, et al (2021) Experimental and Numerical Simulation Study on Seepage Characteristics of complex three-dimensional rough joint Network based on 3D printing Technology [J]. Rock Soil Mech, 42(6): 1659-1668. (in Chinese)
Huang N, Liu R, Jiang Y, et al (2019) Shear-flow coupling characteristics of a three-dimensional discrete fracture network-fault model considering stress-induced aperture variations[J]. J Hydrol 571: 416-424
Ji Y, Wanniarachchi WAM, Wu W (2020) Effect of fluid pressure heterogeneity on injection-induced joint activation[J]. Comput Geotech 123: 103589
Ji Y, Zhuang L, Wu W, et al (2021) Cyclic water injection potentially mitigates seismic risks by promoting slow and stable slip of a natural joint in granite[J]. Rock Mech Rock Eng 1-17
Kohl T, Evans KF, Hopkirk RJ, et al (1997) Observation and simulation of non‐Darcian flow transients in fractured rock[J]. Water Resour Res 33(3): 407-418
Lee CH, Chang JL, Hsu KT (1996) Investigation of hydraulic aperture at surface-exposed rock joints in situ[J]. Geotechnique 46(2): 343-349
Leung CTO, Zimmerman RW (2012) Estimating the hydraulic conductivity of two-dimensional fracture networks using network geometric properties[J]. Transp Porous Media 93: 777-797
Li Y, Sun S, Tang C (2019) Analytical prediction of the shear behaviour of rock joints with quantified waviness and unevenness through wavelet analysis[J]. Rock Mech Rock Eng 52: 3645-3657
Liu Q, Liang B, Li HY, et al (2019) Experimental study on the propagation characteristics of hydraulic fracturing on shale surface [J]. J Liaoning Tech Univ 6
Liu R, Huang N, Jiang Y, et al (2020) A numerical study of shear-induced evolutions of geometric and hydraulic properties of self-affine rough-walled rock fractures[J]. Int J Rock Mech Min Sci 127: 104211
Ma D, Miao XX, Chen ZQ, et al (2013) Experimental investigation of seepage properties of fractured rocks under different confining pressures[J]. Rock Mech Rock Eng 46(5): 1135-1144
Passelègue FX, Brantut N, Mitchell TM (2018) Fault reactivation by fluid injection: controls from stress state and injection rate[J]. Geophys Res Lett 45(23): 12,837-12,846
Quinn PM, Cherry JA, Parker BL (2020) Relationship between the critical Reynolds number and aperture for flow through single joints: Evidence from published laboratory studies[J]. J Hydrol 581: 124384
Rong G, Yang J, Cheng L, et al (2016) Laboratory investigation of nonlinear flow characteristics in rough joints during shear process[J]. J Hydrol 541: 1385-1394
Rutqvist J, Stephansson O (2003) The role of hydromechanical coupling in fractured rock engineering[J]. Hydrogeol J 11: 7-40
Sun KM, Wang S, Zhang SC (2014) Numerical simulation of hydraulic fracturing crack propagation in shale gas reservoir [J]. J Liaoning Tech Univ 33(1): 5-10. (in Chinese)
Xia CC, Yu QF, Qian X, et al (2020) Experimental study on shear-seepage characteristics of rock joints with constant normal stiffness [J]. Rock Mech Rock Eng 41(1): 57-66 (in Chinese)
Xiao WM, Xia CC, Wang W, et al (2013) Seepage analysis of rough joints based on Barton dilatancy model [J]. Ch J Rock Mech Eng 32(4): 814-820. (in Chinese)
Xie L Z, Gao C, Ren L, et al (2015) Numerical investigation of geometrical and hydraulic properties in a single rock joint during shear displacement with the Navier–Stokes equations[J]. Environ Earth Sci 73: 7061-7074
Xing K, Qian J, Zhao W, et al (2021) Experimental and numerical study for the inertial dependence of non-Darcy coefficient in rough single joints[J]. J Hydrol 603: 127148
Zhang Y, Zhao GF (2020) A global review of deep geothermal energy exploration: from a view of rock mechanics and engineering[J]. Geomech Geophys Geo-Energy Geo-Resour 6: 1-26
Zhao Y, Liu Q, Zhang C, et al (2021) Coupled seepage-damage effect in fractured rock masses: Model development and a case study[J]. Int J Rock Mech Min Sci 144: 104822
Zhong Z, Yan JX, Xu CQ, et al (2021) Study on size effect of permeability in rough single crack and its influencing factors [J]. Water Resour Hydropower Eng 52(4): 192-201
Zhou JQ, Wang L, Chen YF, et al (2019) Mass transfer between recirculation and main flow zones: Is physically based parameterization possible? [J]. Water Resour Res 55(1): 345-362
Zhuang L, Jung SG, Diaz M, et al (2020) Laboratory true triaxial hydraulic fracturing of granite under six fluid injection schemes and grain-scale fracture observations[J]. Rock Mech Rock Eng 53: 4329-4344
Zou L, Jing L, Cvetkovic V (2015) Roughness decomposition and nonlinear fluid flow in a single rock joint[J]. International Journal of Rock Mechanics and Mining Sciences 75: 102-118
Zou L, Jing L, Cvetkovic V (2017) Shear-enhanced nonlinear flow in rough-walled rock fractures[J]. Int J Rock Mech Min Sci 97: 33-45

No competing interests reported.

supplement.1.png

Download PDF

Version 1

posted

You are reading this latest preprint version

Numerical simulation of two-dimensional cyclic water injection in rock fractures

Status:

Version 1

Abstract

Figures

Highlights

1 Introduction

2 Governing equations

2.2 Fracture seepage equation

3 Numerical simulation

3.1 Modeling process

3.2 Initial parameters and boundary conditions

3.3 Frequency conversion seepage calculation process

4 Numerical simulation results

4.1 Constant water injection scheme

4.2 Variable amplitude water injection scheme

4.3 Variable frequency injection scheme

4.4 Mainstream domain of rock seepage

5 Conclusions

Abbreviations

Declarations

References

Additional Declarations

Supplementary Files

Status:

Version 1