Abstract

The progressive shear failure of a rock mass under hydromechanical coupling is a key aspect of the long-term stability of deeply buried, high fluid pressure diversion tunnels. In this study, we use experimental and numerical analysis to quantify the permeability variations that occur in an intact marble sample as it evolves from shear failure to shear slip under different confining pressures and fluid pressures. The experimental results reveal that at low effective normal stress, the fracture permeability is positively correlated with the shear displacement. The permeability is lower at higher effective normal stress and exhibits an episodic change with increasing shear displacement. After establishing a numerical model based on the point cloud data generated by the three-dimensional (3D) laser scanning of the fracture surfaces, we found that there are some contact areas that block the percolation channels under high effective stress conditions. This type of contact area plays a key role in determining the evolution of the fracture permeability in a given rock sample.

1. Introduction

In many buried underground construction projects such as deeply buried diversion tunnels, the surrounding rocks are subjected to high geostress and high hydraulic pressure. The complete rock mass may experience shear failure due to high geostress conditions; if the fractures generated by shear failure continue to slip due to ongoing tectonic stresses, the permeability of the rock mass will continue to evolve over time [1]. The high fluid pressure may further trigger shear slip along the fractures in the formation [2, 3], which results in shear compaction or dilation of the fractures or faults and the permeability reduction or enhancement [4, 5]. This physical process can be simply described by the Coulomb-Mohr criterion as follows [6]: where is the shear stress, is the cohesive strength, is the coefficient of friction (also known as the frictional strength), is the effective normal stress, is the normal stress applied to the fracture or fault plane, and is the fluid pressure acting on the fracture. Therefore, by exploring the changes in the rock mass permeability during the entire process of the shear failure of rocks at high confining pressures and fluid pressures, we hope to improve our understanding of the hydraulic characteristics of rock masses, which will provide valuable insight into preventative measures or design ideas that should be incorporated into deep underground engineering projects such as diversion tunnels.

In recent years, many studies have focused on the permeability characteristics of a single fracture [711]. These studies conclude that the permeability characteristics of a fracture can be linked to the physical parameters of the fracture surface, including the aperture, the roughness, the contact area, and the matedness [12]. In addition, the cubic law obtained from the smooth parallel plate test has been modified according to different definitions of roughness [8, 13, 14]. The permeability characteristics of a fracture are related to not only its physical parameters but also to the shear displacement [12, 15]; however, the evolution of the fracture permeability during the shear process has rarely been investigated.

Most of the existing experimental studies on the permeability evolution of fractures during shearing use three-dimensional (3D) printing or poured cement materials to create prefabricated rough fracture surfaces with varying degrees of roughness [6, 1618]. While these empirical studies seek to quantify the permeability and frictional sliding characteristics of fracture surfaces with different roughness during the shearing process, the fracture surface stiffnesses and elastic moduli of these materials are much lower than those of a natural hard rock sample. Other ways to create rough fracture surfaces include splitting rocks such as granite, sandstone, marble, and shale [1921] or preroughened surfaces using other methods to obtain rough fracture surfaces [22, 23]. Then, through shear flow experiments, it is possible to explore the effects of parameters such as normal stress, shear displacement, and surface roughness on the evolution of fracture permeability during the shear process. Their results show that at lower normal stress, as the shear displacement increases, the rough fracture surface exhibits obvious dilatancy characteristics during sliding, which leads to a distinct increase in permeability. As the normal stress increases, the dilatancy of the shear process may be restrained, and then, the increase in permeability would be reduced [24, 25]. If the initial roughness of the fracture surface is high, then the fracture permeability may be enhanced during shear slip [26, 27]. Conversely, when the initial fracture surface is relatively smooth, the fracture permeability decreases with increased shearing [5, 22, 28]. In addition, Fang et al. [6] pointed out that the permeability evolves in a fluctuating pattern for significantly rough fractures. Despite the wealth of information provided by these experiments, none of the prefabricated rough fracture surfaces obtained through the 3D printing or Brazilian splitting methods can truly reproduce the fracture characteristics created by the shear failure of intact rocks under high geostress and hydraulic pressure.

Furthermore, due to the difficulty of the fluid tightness during the shearing process, the fluid pressures in the aforementioned shear flow experiments are typically less than 1 MPa, and the normal stress are relatively small. As such, the permeability results from these experiments are not necessarily applicable in shearing scenarios with high confining pressures and fluid pressures. Additionally, experiments that utilize prefabricated rough fracture surfaces (i.e., those that are split, prepolished, poured, or 3D printed) cannot reflect the permeability evolution during intact rocks from shear failure to slippage.

The purpose of this study is to explore the evolution of the permeability of intact rock samples from shear failure to shear slip under high confining pressures and fluid pressures. First, based off of the characteristics of deeply buried tunnels, we designed and created shear modules suitable for shear flow experiments under high confining pressures and high fluid pressures. We then conducted shear flow experiments on an intact marble core sample obtained from the surrounding rock of a deeply buried tunnel and recorded the permeability changes that occurred over time. Additionally, we developed a numerical model for fluid flow via 3D laser scanning of a fracture surface after shearing. Based on the results of our experiments and numerical models, we have gained significant insight into how the presence of high confining pressures and hydraulic pressures impacts the fracture permeability in intact rock masses during the shear process.

2. Experimental Method

Using a self-designed shear module, we conducted shear flow experiments with different fluid pressures and confining pressures. We used a 3D laser scanner to gather the point cloud data from the fracture surface after shearing. These data provide the experimental basis for our investigation of how the confining pressure and the fluid pressure affect the mechanical properties and flow evolution of a rock mass during rock shear failure and slippage.

2.1. Rock Sample Preparation

The marble cores were collected from the surrounding rock of a diversion tunnel in the Jinping Mountains in Sichuan Province in China. The core samples were processed into several cylindrical samples. In order to reduce short circuiting along the edge of the sample, we evenly coated the sides of the rock sample with a layer of silica gel and let it stand until the silica gel was firm. After placing the cylindrical samples in a vacuum saturation device for 24 hours to evacuate the air, we loaded the samples into the module for testing. During the application of the shear load, a certain amount of osmotic fluid pressure was applied to the inlet. When the shear load increased to form a water-conducting fracture, the water seeped along the evolution path of the shear plane crack. In addition, the fluid pressure provided the normal stress required for fracture expansion, forming the shear flow process.

2.2. Experimental Apparatus

The triaxial shear flow experiments were conducted using a self-adaptive, fully automatic triaxial test machine (located at SINOSTEEL in Anhui Province) and our self-designed shear module, which consists of a pair of shear blocks and a pair of semicircular silicon plugs (Figure 1(a)). The specimen and the shear module were encapsulated by a PVC casing (Figure 1(b)) to avoid the intrusion of the confining oil. The test system mainly consists of three automatic servo pumps, a confining pressure tank, and an automatic data acquisition system. The three pumps (P1, P2, and P3) control the axial pressure, the confining pressure, and the fluid pressure, respectively (Figure 1(b)). The normal stress is equal to the confining pressure applied by pump 2, and the fluid pressure is applied by pump 3. Owing to a large difference in elastic modulus between shear blocks and soft silicone plugs, the specimen experiences shear force during axial loading in the triaxial rig. There are two linear variable differential transformers (LVDTS) on the rigid indenter, and similarly, d1 and d2 are directly installed on the outside of the sample loading table to measure the axial deformation. In order to reduce short circuiting along the edge of the shear module, an O-type sealing ring was also installed on the lower shearing module. We also wrapped the sample in a heat shrinkable PVC casing to prevent the oil from entering the porous rock. The main technical parameters of the device are as follows: the maximum axial pressure of the test system is 100 MPa, the confining pressure is 60 MPa, the fluid pressure is 20 MPa, the maximum flow rate is 15 mL/min, and the loading rate of the displacement control is 0.005–100 mm/min.

2.3. Experimental Plan and Procedure

In this study, we conducted shear flow tests on intact rock samples under different confining pressures and fluid pressures with a fixed shear loading rate. The specific test schemes and sample numbers are shown in Table 1.

This test system, which provides precise control over the shear displacement rate, accurately measures the changes in the permeability and the shear force of an intact rock sample during shear flow. After placing the test specimen in a confined pressure chamber (Figure 1(c)), we then connected the well LVDT displacement gauge and sealed the confined pressure chamber. Using the pumps, the confining pressure was set to normal pressures () of 6 MPa, 10 MPa, or 15 MPa and was kept constant during each test. A constant fluid pressure was set up at the water inlet using a constant pressure pump. The fluid pressure () was set to 0.5 MPa, 2 MPa, or 4 MPa. The permeability of the marble in this area is about 10-20 m2 [29]; therefore, the intact marble sample is initially impervious to infiltration. The axial pressure was loaded at a constant rate of 0.1 mm/min via the axial compression system until the axial displacement reached 4.5–5.5 mm. The measurements were taken continuously and were averaged at recording rates of 1 Hz. When the flow velocity is low, the local cubic law can be used to describe the flow in fractures [30, 31]. The permeability of the fracture was calculated according to the cubic law [7]: where [m2] is the fracture permeability, [m] is the equivalent hydraulic aperture, is the viscosity coefficient of the fluid ( Pa·s at 20°C for distilled water), [m] is the length of the percolation channel, [m] is the fracture width, [m3/s] is the measured flow rate, and [Pa] is the differential pressure between the upstream and downstream areas.

2.4. Surface Profiling

In order to obtain the 3D morphological parameters of the fracture surface after shear failure until shear slip of the intact rock under different confining pressures and fluid pressures, we used a 3D blue light scanning instrument (Wiiboox) (Figure 2(b)) to obtain the 3D point cloud data for the fracture surface after shearing (Figure 2(a)). This 3D optical scanner uses blue light raster scanning technology with a measurement accuracy of 5–15 μm, a single-sided measurement range of  mm2 to  mm2, an average sampling point distance of 0.04–0.16 mm, and a scanning speed of less than 1.5 s. In order to better visualize the surface, we applied brightener to the reflective areas.

Using the surface morphology test data for the fracture specimens, we calculated the fractal dimension and the average aperture between the two separated fracture surfaces of the rock fracture. In addition, we gridded the surface topography measurement data, which made it convenient to import the data into the finite element numerical software. Then, a numerical model of the fluid flow was established to analyze the evolution of the fracture permeability characteristics.

3. Experimental Results and Analysis

In this section, the effects of normal stress and fluid pressure on the variation in the permeability of intact rocks during shearing are discussed using the results of the shear flow tests.

3.1. Evolution of Permeability during Shearing

As shown in Figures 35, the permeability of the sample decreases with increasing normal stress (increased from 6 MPa to 15 MPa), which is due to the normal closure of the fracture aperture caused by the elastic deformation [23, 32]. However, when the confining pressure is 15 MPa, the water pressure is 2 MPa, the shear displacement is 3~4 mm (Figure 6(b)), and the permeability increases greatly. The reason may be that the aperture distribution changes in the original contact area perpendicular to the flow direction, and new flow channels are formed during the shear process. This is also reflected in the numerical simulation in Section 4.2.

Under the same normal stress, the fluid pressure will also affect the fracture permeability. When the fluid pressure reaches 4 MPa, the effective normal stress on the fracture surface is significantly reduced. When the confining pressure is 6 MPa or 10 MPa (Figures 3(c) and 4(c)), while the confining pressure is 6 MPa and the displacement is 3.5 ~ 5.5 mm, the fracture permeability decreases. But the general trend is that the permeability of the sample increases with increasing shear displacement. It is possible that when the shear displacement is 3~4 mm, the new flow channels are generated, which result in the increase of fracture permeability, even more than that caused by shear dilatancy. However, when the confining pressure is 15 MPa, even under 4 MPa of fluid pressure, there is no obvious pattern in the permeability variation with the shear displacement (Figure 5(c)). Ji et al. [33] mentioned that for fractures under high effective normal stress, the fluid pressure will be concentrated at the entrance, resulting in heterogeneity of the fluid pressure distribution in fractures. Therefore, there is no obvious regularity of permeability evolution. When the confining pressure is only 6 MPa (Figure 3(a)), under 2 MPa of fluid pressure, the permeability also increases with increasing shear displacement. However, with a low fluid pressure of only 0.5 MPa, regardless of whether the normal stress is 6 MPa, 10 MPa, or 15 MPa, there is no distinct pattern in the permeability with the shear displacement, and the permeability evolves in a fluctuating pattern.

In general, as shown in Figures 3(b), 3(c), and 4(c), a dilatancy effect occurs during shearing under low effective normal stress, which results in an increase in the permeability with increasing shear displacement. This phenomenon has been demonstrated in previous studies [26, 27]. Under low effective normal stress, the protrusions of the fracture surface may not be cut or worn down. At this time, it is possible to generate an obvious dilatancy effect, which increases the average mechanical aperture, leading to an increase in permeability. Under high effective normal stress conditions (Figures 4(a), 4(b), and 5), the dilatancy characteristics of the fracture are limited, but the high effective stress also enhances the asperities’ shear wear mechanism. At this time, the fracture dilatancy, compression, and shear mechanism of the asperities jointly affect the evolution of the fracture permeability [1], so the evolution of the permeability is more complicated. In addition, under high effective normal stress (Figures 4(a), 4(b), and 5), the relatively rapid decrease in shear stress is generally accompanied by an abrupt increase in permeability. This may be caused by the fact that under high effective stress, some of the large asperities are cut off. This results in the formation of new flow channels and increases the equivalent hydraulic aperture of the fracture surface, leading to an increase in the permeability [23].

3.2. Shear Stress Changes during Shearing

The shear failure of the intact rock samples occurs in three stages. The first stage is the elastic deformation stage, in which the shear stress increases linearly. The second stage is the fracture development stage, in which the shear stress gradually reaches its peak, the penetrating fracture forms, and the rock fails. The third stage is the slip after failure stage, in which the shear stress decreases rapidly from the peak and then decreases slowly with increasing shear displacement. As shown in Figures 35, as the normal stress increases (from 6 MPa to 15 MPa), the shear stress and the residual stress in the slip phase also increase significantly. In addition, because the fluid pressure acts as a surface force on the rock surface, it does not have an effect until a transfixion fracture forms in the intact sample. Therefore, the fluid pressure has little effect on the peak shear stress value. However, after shear failure occurs, the fluid pressure causes the effective normal stress at the fracture surface to decline, and an increase in the fluid pressure reduces the residual stress during the slip phase (Figures 3 and 4). However, when the normal stress reaches 15 MPa (Figure 5), the effect of the fluid pressure on the residual stress decreases. Moreover, when the fluid pressure is 4 MPa (Figures 3(c), 4(c), and 5(c)), the shear stress instantaneously drops after achieving its peak value. At low effective normal stresses, the marble core quickly loses stability after shear failure.

4. Interpretation and Discussion

We evaluated the evolution of the mechanical and flow characteristics of the marble sample during the shearing process. Our test results indicate that under low effective normal stress conditions, the permeability of the shear fracture increases with the shear displacement, while the permeability evolves more erratically under high effective normal stress conditions.

For the results shown in Figures 3(a), 4(a), 4(b), and 5, at high effective stress, there is no obvious regular change in the permeability with the shear displacement; this conclusion is supported by similar work in previous studies [6]. One explanation for this lack of correlation between the permeability and the shear displacement at high effective stresses is that for rough joint surfaces, due to the combined effects of dilation, compaction, and clogging during shearing, the permeability increases overall, but tends to alternate between lower and higher values along that general trend. In addition, Zhang et al. [1] also reported that the change in the fracture permeability during shearing is complex because it is affected by the dilatancy, shrinkage, and the shear wear of the fracture asperities.

In contrast to the previous techniques used to create fracture surfaces, in this study, we used an undisturbed fracture surface obtained through complete core shearing. The matedness of the hanging wall and the footwall of the resulting fracture should be greater than that of the prefabricated fracture surfaces. After the experiment, we found that the void space at the inlets of some of the samples was much smaller than that observed in the middle of the sample (Figure 6). It is questionable whether this phenomenon affects the average permeability of the fracture surface. Therefore, in this paper, we used the point cloud data for the fracture surface obtained through 3D laser scanning to establish a numerical model of the fractures using the COMSOL Multiphysics finite element software and propose a new explanation for the evolution of permeability during the shearing process.

4.1. Numerical Simulation of the Fracture Flow

The fracture surface was three-dimensionally scanned after shear failure, and the fractal dimensions of the hanging wall and footwall were calculated, as shown in Table 2. The physical parameters of the hanging wall and the footwall of the fracture are slightly different from one another after shear failure and shear slippage.

The distributions of contacting asperities/void space (hereinafter called “aperture distribution”) for the rock fracture consist of a pair of fracture surfaces. In this study, the fracture surfaces are assumed to be well matched in the initial state (Figure 7(a)). After implementing shear displacement, the two fracture walls shift horizontally (in the -direction) due to the shear stress and separate vertically (in the -direction) due to the shear dilation (Figure 7(b)). In addition, the asperities of the hanging wall and footwall overlap, the overlapping asperities are assumed to be the contact area (the aperture is zero), and the deformation of the surface is not considered [34, 35]. For each shear displacement interval, the opposing points on the two fracture walls are changed, and the local aperture is recalculated according to the following equation [36]: where is the shear displacement, is the aperture asperity height of the hanging wall, is the aperture asperity height of the footwall, and the value of indicates the contact point between the two fracture surfaces. is the normal displacement generated at , which is the given initial aperture in this paper. Because of the gradually decreasing nominal coincidence area between the hanging wall and the footwall during shear, the aperture field was extracted from the central part of the original model to ensure that the analyzed area was constant (Figure 7(c)) [36]. In addition, we treated the contact-area ratio, which is defined as the ratio of the contact area to the total analyzed area. The fracture is composed of a hanging wall and a footwall (Figure 8(a)), and Figure 8(b) shows the aperture distributions of the fracture after shearing and the boundary conditions of the fracture flow simulation. The distribution of the aperture and the permeability was calculated by simulating the laminar flow of a viscous, incompressible fluid using the Reynolds equation [3740]: where is the aperture, and is the pressure of the fluid.

In addition, due to the presence of microparticles, in order to consider the influence of the contact area perpendicular to the flow direction on the permeability, it is impossible to make full contact between the overlapping asperities of the fracture surface. Thus, the contact area was set to an ultrasmall nonzero aperture of 0.01 μm in the calculation [19], and thus the equivalent hydraulic aperture and permeability must be determined using Equations (3) and (4). The fracture fluid flow was numerically simulated using the finite element software COMSOL Multiphysics. It should be noted that this model does not strictly consider the interaction between the mechanical changes and the flow changes, so it is not a hydromechanical coupling model.

Based on the measurements of the fracture surface topography and the flow rate, the 2-D aperture distributions of the fractures under the confining stresses were numerically determined by the computer through the permeability matching method, in which the pairs of fracture surfaces are in contact with each other, so the aperture distributions have the experimentally determined fracture permeabilities [19]. Figures 9(a)9(i) show the aperture distributions of samples M1-M9, respectively.

These results, when combined with the data shown in Figures 3(b), 3(c), and 4(c), indicate that for samples M2, M3, and M6, where the permeability increases with the shear displacement, the flow rate obtained from the simulation is much larger than the value obtained from the shear tests, despite an initial aperture value of 0. A possible reason for this is that the simulation does not take into consideration the head loss caused by the rough surface or the presence of fine worn particles that could fill the fissure voids, which would decrease in the flow. In the flow direction, there is a contact area that blocks the flow channel in the middle of samples M1, M4, and M7 and at the entrances and exits of samples M5, M8, and M9. The equivalent hydraulic aperture, contact-area ratio, average aperture, given initial aperture, and other parameters of samples M1-M9 are shown in Table 2. Based on a comparison of samples M1, M4, and M7, under the same hydraulic gradient, there is no obvious relationship between the permeability and the contact-area ratio of the average aperture and the other parameters, which contradicts results obtained in previous studies [16, 27, 41]. This contradiction may arise because the matedness of the hanging wall and footwall of a fracture obtained from the shear failure of an intact core differs from that of a disturbed prefabricated fracture surface. It is possible that the effect of such contact areas on permeability was not taken into account in previously published studies.

Figures 10(a)10(i) show the velocity distribution and the main flow channels of the fracture. For samples with contact areas perpendicular to the flow direction, the main flow channels are all oriented along the flow direction where the width of the contact area is the shortest.

4.2. Influence of the Special Contact Area on the Permeability Evolution

As discussed in Section 4.1, the contact areas oriented perpendicular to the flow direction have a significant influence on the permeability. Thus, this type of contact area is probably the reason why the permeability of the test results changes without obvious regularity. Based on this, we changed the distribution of this type of contact area by setting different initial apertures and different shear displacements, and then, the following numerical simulations were constructed to investigate the variation in the permeabilities of the fractures in the different samples.

Figures 11(a)11(i) show the changes in the flow rate per second for the different initial apertures of samples M1-M9, respectively. The permeability values of samples M2, M3, and M6, which do not have contact areas blocking the flow path, increase with a given initial aperture, and the relationship between the aperture and the flow rate roughly satisfies the cubic law. As the initial aperture increases, the contact-area ratio of the sample decreases, the average aperture increases, and the equivalent hydraulic aperture increases. However, for the other six samples, when the given initial aperture is small, there are some contact areas that block the flow path, and the increase in the initial aperture has little effect on the permeability. The permeability values of these samples change by an order of magnitude. Once the given initial aperture reaches a certain critical value, as the initial aperture increases, the evolution of the permeability is the same as that in samples M2, M3, and M6, and it roughly satisfies the cubic law.

We established a numerical model using the surface morphology and a shear displacement of 5 mm to replace the initial surface morphology in order to simulate the evolution of the fracture permeability during the shearing process. The wear on the fracture surface during the shearing process was ignored. Samples M6 and M8 represent the low effective normal stress and high effective normal stress cases, respectively.

According to the curve of the normal deformation versus shear displacement obtained from our experiments (Figure 12), the initial aperture should be given in the numerical simulation when the shear displacements are 1 mm, 2 mm, 3 mm, 4 mm, and 5 mm during the shearing process.

As can be seen from Figures 13 and 14 when the shear displacement of sample M8 is 1 mm and 2 mm, there is no contact area blocking the flow channel. At this time, the average flow rate is larger, as shown in Figure 15(b). Although compared with the shear displacement of 1 mm, the initial aperture for the shear displacement of 2 mm is larger, and the contact-area ratio is lower, which means that there are more void spaces. However, when the shear displacement is 2 mm, the average flow rate is still smaller than that when the shear displacement is 1 mm. This is because a contact area that blocks the flow channel begins to gradually form at the inlet. As the shear displacement increases, the contact area at the inlet is more fully developed, and the average flow rate decreases by four to eight orders of magnitude compared to the flow rates recorded for shear displacements of 1 mm and 2 mm. However, sample M6, which represents the low effective stress case, never has a contact area blocking the percolation channel during shearing (Figures 16 and 17). Compared to sample M8, sample M6 has a smaller aperture, a larger contact-area ratio, and a much higher average flow rate (Tables 3 and 4). Furthermore, with sample M6, the average flow rate increases with the shear displacement (Figure 15(a)).

Based on the results of our numerical simulations, we propose a new explanation for the variation in the flow characteristics during shearing. We conclude that the presence or absence of contact areas that block the flow channels determines the flow characteristics of the fracture. When these contact areas are present, the evolution in permeability is not necessarily related to the shear displacement, the contact-area ratio, or even the hydraulic gradient. When the initial aperture reaches a threshold value, the contact area blocking the flow path disappears, and the relationship between the flow rate and the given initial aperture is roughly cubic.

Furthermore, we believe that during the shearing process, under high effective normal stress, there may be contact areas on the fracture surface that block the flow channels; these contact areas control the average permeability of the fracture. Significant changes in the fracture permeability during the shearing process depend on the generation and disappearance of these contact areas.

5. Conclusions

In this study, we explored the evolution of the rock mass fracture permeability during shearing under high confining pressures and fluid pressures using a new shear module to conduct complete shear flow experiments on intact marble core samples. The experimental results indicate that under low effective normal stress, the permeability increases with the shear displacement (0–5.5 mm). Under high effective normal stress, there is no obvious regularity for the effect of shear displacement on permeability, and the evolution of the fracture permeability fluctuates.

To investigate the reason for this phenomenon, we combined experimental observations with the digital rock fracture modeling of the aperture distribution and fluid flow. Consequently, we determined that under high effective normal stress conditions, the fracture contains some contact areas that block the flow channel. We then investigated the influence of such a contact area distribution on the permeability evolution by varying the initial aperture and the shear displacement. Our results indicate that this type of contact area controls the fracture permeability when a contact area blocking the flow passage exists. Furthermore, the shearing process is accompanied by the generation and disappearance of this type of contact area, which affects the evolution of the permeability, especially in the case of high effective normal stress.

Data Availability

All data included in this study are available upon request by contact with the corresponding author.

Conflicts of Interest

The authors declare that they have no conflicts of interest.

Acknowledgments

The authors acknowledge the support provided by the National Natural Science Foundation of China (U1765204, 41772340).