Mesoscopic Deformation of a Hydraulic and Mechanical Aperture of a Single Fracture under Normal Stress

The deformation of a fracture significantly affects the mechanical and hydromechanical behavior of a fractured rock mass. Numerous efforts have been made to describe the interactions between the fluid flow characteristics and structural deformation of fractures at the mesoscopic scale, and various geological engineering applications will benefit from the results of these studies. Some examples include oil/gas storage and exploitation [1, 2], geological sequestration of CO₂ [3, 4], and geological disposal of nuclear waste [5, 6]. Fluid flow in a single fracture is essential for research on fracture-matrix flow and large-scale field experiments [7-12], which is meaningful to quantitatively analyzing the change in the hydraulic and mechanical aperture in a fractured rock mass under normal stress.

Numerous experimental and numerical studies have shown that the deformation of a fracture aperture under normal stress is highly nonlinear because the contact area increases [13, 14]. Because the flow path changes with the deformation of the fracture aperture under normal stress, the hydraulic aperture of a fracture also has a nonlinear relationship with the normal stress, which indirectly affects the characteristics of the seepage mechanism [15, 16]. Based on the cubic law, the hydraulic aperture can be written as follows [17, 18]:

where Q is the discharge (L³T⁻¹), w is the fracture width (L), e_h is the hydraulic aperture (L), P is the fluid pressure gradient (ML⁻²T⁻²), and µ is the dynamic viscosity coefficient (ML⁻¹T⁻¹). The changes in the mechanical aperture and the fracture morphology must be considered in order to accurately estimate the hydraulic aperture under normal stress. Many researchers have proposed models for predicting hydraulic apertures [19-23]. Mechanical apertures have been measured using many methods and techniques in field tests and laboratory experiments, including digital image correlation [24], three-dimensional laser scanning [25], and electronic speckle pattern interferometry [26]. In principle, the relationship between the fluid flow mechanism’s characteristics and the normal stress is deduced indirectly based on the constitutive law of the fracture closure [27, 28], the hydraulic aperture prediction model, and the cubic law. Since the hydraulic aperture can reflect the geometric properties of a single fracture, the intrinsic permeability or transmissivity was calculated as follows [29]:

where k is the intrinsic permeability (L²), and T is the transmissivity (L³). The relationship between the intrinsic permeability and the normal stress is also available. Many previous qualitative and quantitative studies on the macroscopic scale could not clearly explain the reason for the occurrence of stress-induced changes in the fluid flow characteristics. The mesoscopic parameters reveal the complicated preferential and channel flow behaviors caused by normal stress loading and characterize the changes in these parameters at the macroscopic scale. Conversely, conducting research on the macroscopic-to-mesoscopic scale is much more difficult.

Aiming at the above issues, the closure deformation of a single fracture under normal stress was characterized based on the strain, which is defined as follows:

where ε_m is the linear strain of the fracture (-), e₀ is the initial aperture of the fracture (L), e is the closure of the fracture aperture due to the normal stress (L)], and e₀ – e is the variation in the aperture during the fracture closure progress. The advantage of using the fracture’s strain compared with the fracture’s aperture resides in the following aspects. (1) In many rock engineering projects, strain measurement techniques, for example, distributed fiber-optic strain sensing [30], linear variable differential transformers [31], and strain gauges [32], have been widely used in field tests and have been demonstrated to have a high strain accuracy on the order of the microstrain (με). (2) Based on Equations (2) and (4), the relationship between the intrinsic permeability and the strain in terms of the hydraulic aperture is as follows:

where e_h0 is the initial hydraulic aperture of the fracture (L), and ε_h is the strain in terms of the hydraulic aperture (-). (3) Because ε_h–e_h exhibits a linear trend (see Equation (5)), the results of the research conducted on the hydraulic aperture can for the most part be applied directly to the strain in terms of the hydraulic aperture. On this basis, the deformation of the mechanical and hydraulic aperture (e_m–e_h) can be directly associated with the strain in terms of the mechanical/hydraulic aperture (ε_m–ε_h). On the macroscopic scale, it is harder to measure the aperture than the strain.

In addition to using the fracture strain to characterize the change in the hydraulic aperture under normal stress, the contact area, which is the area of total closure for a fracture, is also important [33-36]. Brown [33] digitally created two parallel rough fracture surfaces and moved them closer together, resulting in the contact area being distributed in the fracture’s void space. A fracture with a fractional contact area of 0.4%–7% still contains some fluid because the fracture cannot close completely under normal stress. Zimmerman [36] investigated the intrinsic permeability of a fracture containing a contact area with the randomly oriented inclusion of an elliptical shape derived based on Maxwell’s effective medium theory, and they found that

where k₀ is the nominal permeability (L²), c is the fractional contact area (-), $β=(1+α)2/4α$⁠, and α is the aspect ratio of the contact area (-). They concluded that the location of the contact area does not play an important role in the permeability until c >0.5 because the isotropic permeability tensor causes a randomly arranged contact area. In an early experimental attempt, Nemoto et al. [34] used pressure-sensitive film to directly measure the contact area. The ratio of the permeability in different directions increased with increasing contact area, indicating that the anisotropic flow became significantly emerging with the contact area percentages >40%. Many researchers now argue that the use of pressure-sensitive film overestimates the contact area due to the thickness of the material. Recently, Pirzada et al. [37] used X-ray transparent triaxial and X-ray microcomputed tomography to obtain the internal structure of a fracture. Xiong et al. [38] characterized a nonlinear flow regime that was affected by the contact area using the non-Darcy flow factor β, the critical Reynolds number Re_c, and the Forchheimer nonlinear coefficient B, and their results suggest the contact area is positively correlated with β and B and negatively correlated with Re_c. The engineering stability of deep rock mass engineering largely depends on the seepage characteristics of the rock joints. Therefore, changes in the stress in a deep rock mass lead to the opening or closing deformation of the joints, which affects the seepage characteristics of the joints. In the case of changes in the ground stress, structural surface deformation, which affects the distribution and movement of underground energy materials, will also lead to instability of underground engineering structures, such as dam seepage and coal mine seepage. Research on the flow channels controlled by the contact area has improved our understanding of the fluid flow behavior in fractures, helping predict the occurrence of nonlinear flow behavior in geological engineering.

The above contents show that fracture deformation under normal stress is characterized by fracture strain and the fractional contact area, both of which play definitive roles in controlling the flow, so it is speculated that ε_m–ε_h (the strain in terms of the mechanical/hydraulic aperture) and c (the fractional contact area) have the potential to qualitatively or quantitatively describe the mesoscopic deformation of a mechanical and hydraulic aperture. Therefore, in this study, we investigated: (1) the change in the fracture morphology distribution in a single mesoscopic fracture under normal stress, (2) the change in the fluid flow characteristic distribution in a mesoscopic single fracture under normal stress, and (3) the interdependence between the flow and stress in a single mesoscopic fracture under normal stress. To answer these questions, two fracture surfaces were created using the randomized form of the Weierstrass function in MATLAB, and COMSOL Multiphysics software was used to establish a mesoscopic rough-walled fracture model and to conduct the numerical simulations.

Because the roughness of a fracture surface is characterized by fractal geometry, the Weierstrass function was used to generate the two fracture surfaces [39, 40]. We used the randomized form of the Weierstrass function to define the corner coordinates of the mesoscopic plane.

where z_ij (x_i, y_j) are the corner point coordinates of the mesoscopic plane, i and j are the labels of the corner points of the mesoscopic plane, C_n is a mutually independent random number with a standard normal distribution, n is the number of random numbers, D and λ are fractal variables, and A_n and B_n are independent random numbers with a uniform distribution of [0, 2π]. The two digital surfaces were created using the randomized form of the Weierstrass function in MATLAB. The rough-walled fracture model was constructed using two digital surfaces in COMSOL Multiphysics. The two digital surfaces and the rough-walled fracture model are shown in Figure 1. For the numerical simulations, a mesh of about 2.4 million tetrahedral elements was used to discretize the mesoscopic rough-walled fracture model.

Many normal stress-induced fracture closure models have been used to describe the effect of the normal stress on the fluid flow in a single mesoscopic scale fracture. Wang et al. [16] proposed a coupling relationship between the fluid flow characteristics and confining pressures based on Goodman’s hyperbolic model to express fracture closure under normal stress [28]. Many analytical and empirical models have been used to represent the relationship between normal stress and fracture closure, and the Bandis model was used in this study [27]:

where ΔV_f is the fracture closure deformation; V_m is the maximum aperture closure deformation; σ_n is the normal stress; and K_ni is the initial normal stiffness of the fracture.

For the rough-walled fracture model, the initial mean aperture was set as 0.327 mm, and the elastic modulus E of shale was set as 1.25 × 10¹⁰ Pa. To describe the effect of the normal stress on the fracture deformation and the fluid flow characteristics, the normal stress was varied from 0.0570 to 1.4770 MPa. The model parameters and simulated cases are presented in Table 1.

The well-known Navier-Stokes (N–S) equations and the continuity equation have been widely used to describe incompressible steady-state fluid flow through a single mesoscopic scale fracture as follows:

where $u→$ is the flow velocity (LT⁻¹); ρ is the fluid density (ML⁻³); and $F→$ is the body force vector (MLT⁻²).

To clarify the change in the fracture morphology distribution of a single mesoscopic fracture under normal stress, the contact area and fracture aperture distribution during the fracture closure were obtained through model simulations. Figure 2 qualitatively shows the changes in the contact area distribution and size in six cases with normal stresses varying from 0.1331 to 1.4770 MPa. The contact area can be obtained by sweep and plane interception on COMSOL Multiphysics, and the size and proportion of the contact area can be obtained by integral. At a certain closing distance of the fracture surfaces, geometric contacts at the corresponding position will occur. The aperture after deformation is equal to the original aperture minus the closing distance. If the original aperture is less than or equal to the closing distance, this point presents a contact state. As the normal stress applied to the fracture gradually increased, the upper and lower surfaces moved closer together, and the total closure area increased.

Further quantitative analyses are presented in Figure 3; the contact area A_c (blue solid line), the projection area of the fracture A_f (green solid line), and the fractional contact area c (red dash line) all changed significantly when the normal stress σ_n was greater than 0.45 MPa. The total projection area is equal to the sum of A_c and A_f, which is subject to the size of the model and was equal to 100 mm² in this study. Thus, A_f decreased as the two fracture surfaces moved closer together because σ_n increased. c is defined as the ratio of A_c to A_f [36], which decreases with increasing σ_n. According to c, the fracture closure can be divided into three stages: (1) when c = 0, the fracture is completely open, that is, the total closure area is zero; (2) when 0 < c < 1, the fracture is partially closed; and (3) when c = 1, the fracture is completely closed.

To obtain the mechanical mean aperture under different normal stresses, different normal stress-induced distributions of a mechanical aperture in a rough-walled fracture were considered. For the flow behavior, the distribution of the mechanical aperture determines the void space in the fracture, and thus, it partially controls the state of the fluid. Figure 4 illustrates the distribution of the mechanical aperture under seven different normal stresses (the mechanical aperture ranges from 0 to 0.6 mm, where an aperture equal to zero indicates that the sides of the fracture are in contact, and the aperture is zero at this location). As the normal stress on the rough-walled fracture increased, the two walls of the fracture moved closer together, and the distribution of the aperture became more concentrated. After σ_n reached 0.2602 MPa, areas with apertures of zero (blue areas), that is, where the fracture walls were in contact, began to appear, and the zones with large apertures (magenta areas) began to decrease. When σ_n increased to 0.7176 MPa, the locations with apertures of zero coalesced into the area of total closure, and the aperture in these spots was almost less than 0.4 mm. As σ_n increased to 1.0607 MPa, the area of total closure accounted for a larger proportion of the fracture area, and these areas were surrounded by small aperture zones (light blue areas). The increase in σ_n resulted in a discrete distribution of the aperture. This directly led to heterogeneous transmissivity and indirectly led to complicated preferential and channel flow behaviors [41]. Consequently, the higher normal stress enhanced the complicated preferential and channel flow behaviors.

Figure 5 presents the probability density function (PDF) of the fracture aperture under seven different normal stress values. The green curves in Figure 5 are the normal PDFs, and Std denotes the standard deviation of the fracture aperture. The locations with apertures of zero were not integrated into the mechanical aperture statistics. The distributions of the mechanical aperture under different normal stresses are well approximated by the normal distributions with different mean and standard values. As σ_n increases, the mean decreases, and the Std generally increases. Under normal stress loading, the distribution of the mechanical aperture maintained a normal distribution. Similarly, the rule of the mechanical aperture maintaining a normal distribution did not change with the contact area based on Figures 2, 3, 5. As the normal stress increased, the mesoscopic rough-walled fracture model underwent a nonlinear contact progress, which depended on the roughness of the two digital surfaces and the initial aperture distribution of the fracture. The closure of the fracture increased with increasing normal stress, and the small aperture zones (light blue areas) in the model closed first (Figure 4). Thus, the history of the normal stress applied to the fracture is quantitatively expressed by the initial distribution of the aperture. The continuous progress of the fracture closure is determined by the fracture morphology distribution, which has important theoretical value in research on fracture deformation.

The simulated velocity field and streamlines in the X-direction are shown in Figure 6 for the mesoscopic rough-walled fracture models, with σ_n of 0.0570–1.4770 MPa. Based on the results of the mesoscopic rough-walled fracture models, the velocity distribution is significantly affected by the normal stress σ_n. Based on the influence of σ_n on the fracture morphology distribution (Figures 2, 3,4, and 5), the fluid flow under different σ_n was analyzed. For σ_n <0.2602 MPa, the fractional contact area c was equal to zero, and the fracture was completely open (blue spots), that is, the two walls were not in contact. The low-velocity zones (blue areas in Figures 6(a) and 6(c)) were surrounded by areas with a relatively smaller aperture (light blue areas in Figure 4), and the zone with reverse/eddy flow was not clearly visible on the streamlines (Figures 6(b) and 6(d)). When σ_n increased to 0.7176 MPa (Figures 6(e), 6(g), and 6(i)), c was less than 0.5, and areas with apertures of zero (blue spots in Figure 4), that is, where the walls were in contact, appeared. The relatively smaller aperture areas (light blue spots in Figure 4) were surrounded by the low-velocity zones (blue areas in Figures 6(a) and 6(c)), and they gradually closed and coalesced with the areas of total closure (blue areas in Figure 4). The zones with reverse/eddy flow were distributed behind the areas of total closure (Figures 6(f), 6(h), and 6(j)). When σ_n increased from 1.0607 to 1.4770 MPa, the streamlines revealed that the fluid flow became more concentrated due to the limitations imposed by the areas of total closure (Figures 6(l) and 6(n)). In the channels between two areas of total closure, as a result of the decrease in the cross-sectional area of the flow, higher velocity zones were formed (red areas in Figure 6), which further promoted the formation of local eddies (Figures 6(k) and 6(m)). It should be noted that the fluid flow in the mesoscopic rough-walled fracture models exhibited complicated preferential and channel flow behaviors due to the increase in the contact area under normal stress loading.

Figure 7 shows the fluid flow in the mesoscopic rough-walled fracture from the perspectives of the mean velocity u_m (green line) and velocity difference $|Δu|$ (red line), where $|Δu|$ is defined as follows:

where $|Δu|$ is the velocity difference (LT⁻¹), u₁ is the outlet velocity (LT⁻¹), and u₂ is the inlet velocity (LT⁻¹). For fluid flow in the X-direction, $|Δu|$ and σ_n are positively correlated, while u_m and σ_n are negatively correlated, implying that the decrease in the velocity in the mesoscopic rough-walled fracture models increased with increasing σ_n. In contrast, the decrease in u_m was due to the appearance and increase in the low-velocity zone/area of total closure under the normal stress applied to the two surfaces of the mesoscopic fracture.

Furthermore, Figure 8 illustrates the relationship between $|Δu|$ and σ_n for fluid flow along the Y-direction in the mesoscopic rough-walled fracture models, and the difference in $|Δu|$ between the X- and Y-direction is obvious, suggesting that normal stress-induced flow anisotropy exists. When σ_n is 1.0607 MPa, $|Δux|$ = 1.388 m/s and $|Δuy|$ = 1.3572 m/s, which agree with the fact that the velocity differences in the different flow directions are close. In contrast, the tortuosity of the streamlines is different between the partial fluid fields of the fluid flow along the X- and Y-direction.

Under normal stress loading, the contact points formed by the contact of the upper surface and lower surface coalesced to form the area of total closure, which promoted preferential and channel flow behaviors in the fluid field. Thus, the σ_n-induced deformation of the fracture led to changes in the fluid flow characteristics, that is, Q, P, and k. Fidelibus [42] proposed that the relationship between the hydraulic aperture e_h and closure U_n is as follows:

where f_c is the coupling factor (L). Using Equations (4), (5), and (13), we can obtain the relationship between the strains in terms of the mechanical/hydraulic aperture ε_m–ε_h as follows:

where $c=em0eh0$⁠. Furthermore, the intrinsic permeability k can be obtained as follows:

The relationship ε_m–ε_h (Equation (14)) was used to fit the data obtained from the mesoscopic rough-walled fracture models (Figure 9). In Figure 9, the blue line is the ε_h of the flow in the X-direction, and the gray line is the ε_h of the flow in the Y-direction. Based on the regression results, the ε_m–ε_h relationship is linear, and f_c is the same in the X- and Y-direction. Thus, the flow direction has no impact on f_c, and its effect on c is confirmed. The correlation coefficients (R²) of ε_m–ε_h for the X- and Y-direction flow are 0.9788 and 0.9979, respectively, indicating that the relationship ε_m–ε_h can be expressed as in Equation (14). For k, Figure 10 shows that the R² values of the flow in the X- and Y-direction flow are 0.9849 and 0.9989, respectively. The k-values of mesoscopic rough-walled fracture models under σ_n can be predicted well using Equation (15). The quantitative relationship in Equation (14) describes the deformation of the hydraulic and mechanical aperture of a mesoscopic rough-walled fracture, which reveals the changes in the flow and stress on the mesoscopic scale. Based on the mesoscopic scale ε_m–ε_h, the relationship between the intrinsic permeability k and the strain ε on the macroscopic scale is shown in Equation (15), which can be applied in performance assessment and optimization design in geology engineering. The changes in the fluid flow mechanism characteristics and fracture morphology under different normal stresses are described quantitatively in theory. In practice, this quantitative relationship predicts the changes in the flow and permeability in fractures under changing in situ stress.

To consider the interdependence between flow and stress in a single fracture, the mechanical aperture and hydraulic aperture data of experiments with different lithologies are used for verification (Figures 11 and 12 [21, 43]).

The hydraulic aperture and mechanical aperture data obtained from the experiments were fitted using the Equations (14) and (15) proposed in this paper. The fitting results are shown in Figures 11 and 12. Obviously, the equation model in this paper can effectively describe the interdependence between flow and stress under normal stress loading. Therefore, this comparison demonstrated that the proposed theory can be used to estimate the interdependence between the seepage behavior of rocks and stress with reasonable accuracy. The suggested theory between stress and flow in a fracture under normal stress loading is consistent with results from other experiments. Therefore, it has important theoretical guiding significance for the stress and flow behavior in a fracture of related practical projects.

In this study, we investigated the impact of the normal stress on the fracture morphology distribution and fluid flow and determined the relationship between the strains in terms of the mechanical/hydraulic aperture ε_m–ε_h. The main conclusions of this study are as follows.

1. During normal stress loading, the area of total closure increases as σ_n increases. The distribution of the mechanical aperture is approximated well by the normal distribution. The fluid flow exhibits complicated preferential and channel flow behaviors.

2. As σ_n increases, the mean value decreases, and the Std generally increases. Whether or not the distribution of the mechanical aperture conforms to a normal distribution is not determined by the applied normal stresses. The mean velocity increases, and the velocity difference decreases. In addition, normal stress-induced flow anisotropy exists.

3. Strains in terms of the mechanical/hydraulic aperture can be predicted using the proposed linear relationship, in which the parameter c is affected by the flow direction, but the parameter f_c is not. The relationship between the intrinsic permeability and strain is clear, and it has practical application value in various geological engineering applications.

This work was supported by the National Natural Science Foundation of China (Grant Nos. 41831289, 41772250, and 42072276). The authors also thank the anonymous reviewers for their helpful comments and suggestions.

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.

The data that support the findings of this study are available from the corresponding author upon reasonable request.