Studies have focused on describing the interactions between the fluid flow characteristics and structural deformation of fractures at the mesoscopic scale, which is a scale between the macroscale and the microscale. In this work, a three-dimensional numerical simulation based on the Navier-Stokes equation was carried out to investigate the effect of normal stress on the fracture morphology distribution, the fluid flow characteristics distribution, and the interdependence between the flow and stress in a single mesoscopic fracture. Two fracture surfaces of a mesoscopic rough-walled fracture model were created. Results suggest that the nonlinear relationship between the normal stress and deformation due to the area of the total closure increases unevenly. Distributions of the mechanical aperture are approximated well by a normal distribution. Change in the fluid flow is due to the increase in the fractional contact area. The low-velocity zones are surrounded by relatively smaller apertures, which gradually close and join the areas of total closure. Under the limitation of the total closure areas of the two fracture surfaces, the appearance of channel flow behavior. Compared with the flow in the X- and Y-directions, normal stress-induced flow anisotropy occurred. The mesoscopic quantitative relationship between the strains in terms of the mechanical/hydraulic aperture was determined and proven. A macroscopic relationship between the intrinsic permeability and the strain was deduced, which enhances the evaluation and design of various geological engineering applications in which fracture deformation is considered.
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 CO2 [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 (L3T−1), w is the fracture width (L), eh is the hydraulic aperture (L), P is the fluid pressure gradient (ML−2T−2), and µ is the dynamic viscosity coefficient (ML−1T−1). 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 , three-dimensional laser scanning , and electronic speckle pattern interferometry . 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 :
where k is the intrinsic permeability (L2), and T is the transmissivity (L3). 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 (-), e0 is the initial aperture of the fracture (L), e is the closure of the fracture aperture due to the normal stress (L)], and e0 – 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 , linear variable differential transformers , and strain gauges , 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 eh0 is the initial hydraulic aperture of the fracture (L), and εh is the strain in terms of the hydraulic aperture (-). (3) Because εh–eh 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 (em–eh) 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  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  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 k0 is the nominal permeability (L2), c is the fractional contact area (-), , 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.  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.  used X-ray transparent triaxial and X-ray microcomputed tomography to obtain the internal structure of a fracture. Xiong et al.  characterized a nonlinear flow regime that was affected by the contact area using the non-Darcy flow factor β, the critical Reynolds number Rec, and the Forchheimer nonlinear coefficient B, and their results suggest the contact area is positively correlated with β and B and negatively correlated with Rec. 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.
2. Theoretical Background
2.1. Generation of Rough Rock Surfaces
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 zij (xi, yj) are the corner point coordinates of the mesoscopic plane, i and j are the labels of the corner points of the mesoscopic plane, Cn is a mutually independent random number with a standard normal distribution, n is the number of random numbers, D and λ are fractal variables, and An and Bn 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.
2.2. Fracture Closure
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.  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 . 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 :
where ΔVf is the fracture closure deformation; Vm is the maximum aperture closure deformation; σn is the normal stress; and Kni 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 × 1010 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.
2.3. Fluid Flow
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 is the flow velocity (LT−1); ρ is the fluid density (ML−3); and is the body force vector (MLT−2).
3. Results and Discussion
3.1. Fracture Morphology Distribution of a Mesoscopic Fracture During Normal Stress Loading
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 Ac (blue solid line), the projection area of the fracture Af (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 Ac and Af, which is subject to the size of the model and was equal to 100 mm2 in this study. Thus, Af decreased as the two fracture surfaces moved closer together because σn increased. c is defined as the ratio of Ac to Af , 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 . 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.
3.2. Fluid Flow in a Mesoscopic Fracture Dduring Normal Stress Loading
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 um (green line) and velocity difference (red line), where is defined as follows:
where is the velocity difference (LT−1), u1 is the outlet velocity (LT−1), and u2 is the inlet velocity (LT−1). For fluid flow in the X-direction, and σn are positively correlated, while um 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 um 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 and σn for fluid flow along the Y-direction in the mesoscopic rough-walled fracture models, and the difference in between the X- and Y-direction is obvious, suggesting that normal stress-induced flow anisotropy exists. When σn is 1.0607 MPa, = 1.388 m/s and = 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.
3.3. Interdependence between Flow and Stress in a Mesoscopic Fracture Under Normal Stress Loading
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  proposed that the relationship between the hydraulic aperture eh and closure Un is as follows:
where . 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 fc is the same in the X- and Y-direction. Thus, the flow direction has no impact on fc, and its effect on c is confirmed. The correlation coefficients (R2) 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 R2 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.
3.4. Verification of Interdependence between Stress and Flow in a Fracture under Normal Stress Loading
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.
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.
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.
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 fc 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.
Conflicts of 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.
The data that support the findings of this study are available from the corresponding author upon reasonable request.