An experimental investigation of the nonlinear gas flow and stress‐dependent permeability of shale fractures

The fluid flow characteristics in shale fractures are of great significance for shale gas reservoir evaluation and exploitation. In this study, artificial tension fractures in shale were used to simulate the hydraulic fractures formed by fracturing, and a gas flow test under different pressure gradients was conducted. The nonlinear gas flow and stress‐dependent permeability characteristics were analyzed. The experimental results show the following: (a) CO2 flow in shale fractures exhibits strong nonlinearity. Forchheimer's law, which considers gas compressibility, satisfactorily describes the nonlinear relationship between the flow rate and the pressure gradients in shale fractures. (b) The permeability sensitivity of shale fractures under stress is very strong, and the exponential relationship better describes the pressure dependency of the permeability for the tested shale samples. The permeability of the shale fractures is similar when measured parallel or perpendicular to bedding. Furthermore, the pressure dependence of fractures in shale obeys the Walsh permeability model. (c) As the effective stress increases, the nonlinear flow behavior appears earlier. Based on the Reynolds number and the nonlinear coefficient, a friction factor model is proposed. (d) The normalized transmissivity exhibits a strong correlation with the Reynolds number. CO2 flow through shale fractures is generally dominated by transitional flow. The critical Reynolds number ranges from 1.8 to 102.88 and decreases with increasing effective stress.


| INTRODUCTION
Shale gas is an important unconventional energy source, which is abundant and widely distributed in China. It is expected to become the most important substitute for conventional oil and gas resources. In recent years, substantial breakthroughs have been made in shale exploration and development of the Wufeng-Longmaxi Formation in the Sichuan Basin and its peripheral areas in southern China, indicating a good potential for resource development. For conventional reservoirs, interconnected pores contribute a lot to the reservoir permeability. However, the permeability of shale is very low in the range of nano-Darcy to micro-Darcy, 1 so reservoir stimulation must be carried out to achieve industrial production. At

| 2809
SHEN Et al. present, hydraulic fracturing is mainly used to stimulate shale reservoirs, producing complex fractures and improving the gas desorption and flow conductivity. 2 Since fractures are the main gas flow channels, fracture conductivity is the key to gas production in gas wells.
In recent years, numerous studies have been conducted to understand flow behavior in a single fracture generated in different rock types, including coal, 3 sandstone, 4 and granite. 5 Previous studies on shale fractures have mostly focused on the hydromechanical behavior of propped shale fractures, such as smooth shale plates, steel plates, artificial fractures, and natural fractures. 1,[6][7][8] However, there are unsupported natural and man-made fractures in shale reservoirs, and little research has been conducted on these unsupported fractures. Previous studies have shown that shale permeability is greatly affected by the effective stress. Fractures can greatly improve shale reservoir permeability, and the stress sensitivity of fractured shale is greater than that of intact shale. [9][10][11][12] Guo et al 13 studied the influence of fracture surface roughness, confining stress, proppant type, and distribution on fracture permeability using core fractures in shale from the Shengli Oilfield. They found that compared with the shale matrix, the permeability of the matching fractures and the supporting fractures increased by about one to three and two to four orders of magnitude, respectively, under 35 MPa stress. Zhou et al 14 studied the influence of filling minerals, fracture surface morphology, and acidizing treatment on the permeability of natural shale fractures, and found that the overall roughness and the local roughness of a fracture surface affect permeability, while the size and hardness of the cemented particles affect the stress sensitivity. Most previous studies conducted on fracture flow in shale are based on the cubic law, in which the flow rate is linearly proportional to the pressure gradient. 15,16 The test stress is generally low, 15,17 and the proppant embedment and stress sensitivity of the shale fractures are commonly taken into consideration. [18][19][20] In addition, most shale gas reservoirs in China have burial depths of more than 3500 m. 12 Thus, high stress generally exists, especially in the Sichuan Basin. The complex stress conditions have an important impact on fluid flow in shale fractures. In the process of shale gas production, the fluid pressure in the fractures decreases, so the effective stress on the fracture increases, which leads to the closure of reservoir fractures. Thus, the interaction between the stress and fluid pressure makes the flow in shale fracture more complex.
Shale fractures have a heterogeneous spatial geometry and flow channel. This makes it easy to cause friction and non-negligible inertial losses, so non-Darcy flow commonly occurs in fractures. Non-Darcy flow significantly reduces the ability of fluid conduction in fractures, thus affecting the productivity of hydraulic fracturing wells. 21 Schrauf et al 22 and Yang et al 19 used gas flow tests on single fractures in rocks to demonstrate that there is a nonlinear relationship between the gas pressure gradient and the fluid volumetric flow rate. Chen et al 23 studied the nonlinear flow behavior caused by the roughness of fractures, the contact condition, and the change in the flow passage through rock fractures. Su et al 24 further experimentally confirmed the nonlinear relationship between the flow rate and gas pressure gradient and found that the nonlinear relationship overestimates the permeability at high pressures. In recent studies, the friction factor involved in nonlinear flow through fractured rocks has been widely studied. 25 Li et al 26 studied the friction effect of gas flow in shale fractures and established a friction factor model based on the relative roughness and Reynolds number. Su et al 27 found that rock fracture permeability is affected by the contact area and void space and established a nonlinear flow model that takes into consideration the evolution of the fracture geometry under stress. Although the nonlinearity of fluids in rock fractures has been widely studied, only a few studies have investigated the nonlinear evolution and range of the critical Reynolds number in shale fractures. Understanding the nonlinear flow characteristics of gas in shale fractures is very important as it is meaningful for resource evaluation and gas productivity analysis.
In this study, shale samples with a single fracture were prepared in the direction parallel to the bedding and perpendicular to the bedding, and the nonlinear gas flow and stress-dependent permeability characteristics of the shale fractures were analyzed. CO 2 flow tests were carried out under different confining stresses , and the inlet gas pressure varied from 0.2 to 2.2 MPa.

SPECIMENS
The test sample was collected from an outcrop in Changning, Sichuan Province. The sampling location is shown in Figure 1A. Changning, which is located on the edge of the Sichuan Basin, is an important shale gas-producing area in China. The shale of the Longmaxi Formation and the Wufeng Formation is widely developed in this area. The mechanical and permeability parameters of the shale in this area are shown in Table 1. The shale cores were taken parallel and perpendicular to the bedding plane. 28 All of the samples were cut into 50 × 50 mm cylinders. The upper and lower ends of the sample were placed on the precision grinding machine to polish the surface and ensure that the surface fluctuation of the ends of the samples was within ±0.05 mm. Based on the study conducted by Deng et al, 29 most of the fractures created by hydraulic fracturing are tensile fractures, so we used a selfmade Brazilian splitting device to split the shale specimens, forming a fresh tensile fracture approximately parallel to the axial direction. The process of fracture generation is shown in Figure 1B, and the shale fractures are shown in Figure 1C.
The surface morphology of a rough fracture has an important influence on the fracture transmissivity. A high-precision, noncontact, 3D scanner was used to scan and measure the surface characteristics of each shale fracture. The scanning point distance was 0.1 mm, and the scanning accuracy was 0.02 mm. The three-dimensional coordinate data for the surface points of each scanned specimen were imported into the ParaView software for data analysis and visual interpretation. Figure 2 shows the fracture surface morphology of each specimen. The maximum asperity in the fracture parallel to the bedding was much smaller than that perpendicular to the bedding because the generation of the shale fracture perpendicular to the bedding was disturbed by the bedding. When the fracture is parallel to the bedding, the fracture directly passes through the bedding, and the surface of the fracture is relatively smooth.

| Gas flow tests
The gas flow tests conducted on each single shale fracture was conducted using a triaxial test system. A structural diagram of the triaxial loading test system is shown in Figure 3. During the test, the axial stress and confining stress were accurately controlled using a computer data acquisition unit, and the gas flow rate was measured and recorded by a digital electronic flowmeter. The sample was placed in a heatshrinkable tube and fixed in the triaxial chamber using the upper and lower pedestals. In addition, two O-rings were placed on the upper and lower loading pistons to ensure that the oil could not infiltrate into the sample. The gas flow tests conducted on all of the shale fractures were carried out at room temperature (25°C), at which the dynamic viscosity of CO 2 is 1.49 × 10 −5 Pa/s. Under the condition of constant temperature and low gas pressure, the dynamic viscosity of CO 2 is considered to be constant. During the test, the axial stress was maintained at 3 MPa, the inlet gas pressure and confining stress were simultaneously varied, ensuring that the effective stress in the fracture remained unchanged in each stage. The depth of the shale reservoir in the Sichuan Basin is about 1500-3000 m, and the shale is in a high-stress environment, 30 so the effective stress was gradually increased from 1 MPa to a maximum of 30 MPa; the inlet gas pressure was varied from 0.2 to 2.2 MPa. During the test, it was ensured that the gas pressure was less than the confining stress. The development of the loading process in the triaxial loading test is shown in Figure 4. The applied confining stress was calculated using the principle of effective stress (equation 1). 31 The Biot coefficient was assumed to be one in the fractured medium. The gas pressure in the fracture was assumed to be the average value of the inlet pressure and the outlet pressure, so the effective stress can be expressed by equation (2): where σ is the effective confining stress (MPa); ′ is the applied confining stress (MPa); is the Biot coefficient; P is the gas pressure within the fracture (MPa); P in is the gas pressure at the inlet (MPa); and P out is the gas pressure at the outlet (MPa).

| Test results
In this test, the inlet gas pressure varied from 0.2 to 2.2 MPa, and the CO 2 flow rate at the outlet was measured by the flowmeter. Figure 5 shows the relationship between the pressure gradient and the volumetric flow rate at the outlet under different effective stresses. The flow patterns of the CO 2 in the four shale fractures are basically similar. The gas volume flow at the outlet increases as the loading pressure gradient increases, but the growth rate decreases. In addition, the relationship between flow rate and pressure gradient is not strictly linear, which is mainly due to the fluid needs more energy to overcome the pressure loss caused by the inertial force. As the effective stress increases, the slope of the curve becomes steeper, especially in the initial stage. This is because the increase in effective stress leads to fracture closure and increased contact area, which leads to a decrease in the overflow flow of the fracture. In addition, when the effective stress is 30 MPa, the sample also has a certain flow capacity, indicating that the fracture is still the main migration channel for shale gas. Since the fracture surface is rough, the CO 2 flow between the two sides of the fracture will be subject to significant resistance along the way. Forchheimer introduced a quadratic equation, the quadratic coefficient of which can represent the nonlinear flow characteristics of the fluid in a fracture (equation 3). When CO 2 flows in a rough fracture, since it has a strong compressibility and the stress distribution and density of CO 2 vary in shale fractures, the compressibility must be considered. Based on the mass conservation theorem (equation 4) and the gas equation of state, Su et al 27 proposed a formula to describe the relationship between the pressure gradient and the flow rate that takes into consideration the gas compressibility (equation 5). The parameters A and can be obtained by fitting. The fitting results are shown in Table 2: where ∇P is the pressure gradient (MPa/m); A is the linear coefficient (10 11 kg s −1 m −5 ); Q is the volumetric flow rate (m 3 /s); is a nonlinear coefficient (10 16 /m −5 ); is the fluid density (kg/m 3 ); out is the fluid density at the outlet (kg/m 3 ); Q out is the volumetric flow rate at the outlet (m 3 /s); m out is the mass flow rate at the outlet (kg/s); and L is the test specimen length (m).
A good correlation was found for all of the correlation coefficients (R 2 ) with values greater than 0.97, which reveals that equation (5) can satisfactorily describe nonlinear gas flow in shale fractures. The coefficients A and reflect the linear and nonlinear parts, respectively. The linear coefficient A is related to the viscous force of the interaction between the gas and the fracture, which represents the flow capacity of the rough fracture. The hydraulic aperture e h and the intrinsic permeability of the rock fracture (k = e 2 h /12) can be calculated using equation (6). 32 The nonlinear coefficient is an important parameter for investigating the nonlinearity of the flow because it represents the inertial effect of the irreversible loss of kinetic energy caused by the change in the velocity of the fluid within the fracture. The inertial effect describes the deviation of the fluid flow from the cube law due to the surface roughness in the fracture. The changes in A and with effective stress are shown in Figure 6, which shows that A and increase by 1-3 and 2-4 orders of magnitude, respectively: where w is the width of the fracture surface perpendicular to the flow direction (mm), is the dynamic viscosity of CO 2 (Pa·s), and e h is the equivalent hydraulic aperture (µm).

| Permeability characteristics
The matrix permeability of shale is very low (Table 1), so its influence on the outlet volumetric flow rate is relatively small compared with that of fracture permeability, and thus, it is assumed to be ignored in this study. In addition, the average free path of a CO 2 gas molecule is 0.98i, which is far less than the aperture of the shale fracture, so the slippage of gas in the shale fracture is considered to be negligible in our study. 12 Figure 7 shows the evolution curve of permeability with effective stress, which indicates that the fracture's permeability is very sensitive to changes in the effective stress. The permeability decreases with increasing effective stress, and permeability and stress exhibit a negative exponential relationship. In the initial stage of loading, the permeability decreases rapidly with increasing effective stress. After 8 MPa, the decrease in permeability slows down. When the effective stress increases from 1 to 8 MPa, the permeabilities of the four shale samples decrease by 70.95%, 76.89%, 56.39%, and 76.62%. The damage of the initial stress on the permeability is greater. After 8 MPa, as the contact area increases, the ability of the fracture to resist deformation increases, and the permeability decreases slowly. Hydraulic fracturing of a shale reservoir mainly extends along the bedding plane or perpendicular to the bedding plane. 28 The permeability in both directions is very important for gas well productivity. The average permeability k in the two directions is used as a reference to compare the permeability of the fractures in the direction parallel to the bedding and perpendicular to the bedding. In equation (7), k is the ratio of the average permeability in the direction perpendicular to the bedding plane (k ⊥ ) to that parallel to the bedding plane (k 11 ). The permeability of shale fractures is similar in both directions, and the ratio increases from 0.762 to 1.735 with increasing effective stress. Before the effective stress reaches 12 MPa, the permeability parallel to the bedding plane is greater than that perpendicular to the bedding plane, which may be related to the fact that a small amount of debris roughens the channels in the bedding plane in Figure 2. When the effective stress is higher than 12 MPa, the permeability perpendicular to the bedding plane is larger than that parallel to the bedding plane. This may be due to the fact that when the fractures are made perpendicular to the bedding plane, the asperity of the surface of the fracture is large. Under high effective stress, the normal rigidity of the shale fracture is large and the deformation is small, so the residual aperture is  large. In addition, the fracture's permeability is still very high under an effective stress of 30 MPa: where k ⊥ is the average permeability in the direction perpendicular to the bedding plane (1e −10 /m 2 ) and k 11 is the average permeability in the direction parallel to the bedding plane (1e −10 /m 2 ).

| Permeability sensitivity of shale fractures under effective stress
Fracture permeability is one of the key factors influencing gas production. It is of great theoretical and practical significance to studies of the stress sensitivity of fractures and for revealing the law of production change and determining a reasonable development plan for shale gas wells. Since there is no specific stress sensitivity standard for shale permeability, we refer to China's oil and gas industry standard (SY/ T5336,5358,6385) for analysis. The increase in the effective stress leads to compression of the fractures and a decrease in the permeability. The fracture compression coefficient is a parameter that characterizes the compression ability of the fracture. 1 By fitting the effective stress and permeability data with equation (8), the compressibility of the fracture can be obtained. The fitting results are shown in Table 3: where k is the permeability (1e −10 /m 2 ), k 0 is the initial permeability of the shale fracture (1e −10/ m 2 ), and C f is the fracture compression coefficient (MPa −1 ), which is a parameter that describes the stress sensitivity of the fracture.
The fitting correlation coefficient is >0.944. The compressibility coefficient of the shale fracture is between 0.034 and 0.07 MPa −1 . The reason why the coefficients differ from each other may be related to the surface morphology and contact state. The compressibility coefficient obtained in this paper is highly consistent with that of the Longmaxi deep shale (0.014-0.077 MPa −1 ) under 4 to 80 MPa of confining stress. 12 However, the shale compressibility coefficient obtained in this study is smaller than that obtained from supporting shale fractures (0.041-1.01 MPa −1 ), which are located in the same area of the Sichuan Basin. 1 The reason for this may be that proppant supporting shale fractures are easily deformed and embedded.
Shale fractures are very sensitive to changes in stress, so the change in stress must be quantified to characterize the degree of damage to the shale fracture's permeability. We define the permeability damage rate (DR) as equation (9), which reflects the proportion of the permeability damage under the effective stress: where k 1 is the permeability when the confining stress is 1 MPa, and k i is the permeability when the effective stress is i MPa. Figure 8 illustrates how the permeability damage rates change with effective stress. The permeability damage rate increases rapidly in the initial stage. The increase slows after 12 MPa, gradually becoming constant. This occurs because the two fractures contact each other and form a self-supporting flow channel in the initial stage of loading. Since the contact area is small in the initial stage, the rate of increase of the effective stress is higher, leading to a larger closure of the

F I G U R E 6 Variation in coefficients
A and β with effective stress fracture. In addition, when the effective stress was 30 MPa, the permeabilities of the four shale samples decreased by more than 90%, indicating that the permeability of the shale fractures is greatly affected by the effective stress. The permeability sensitivity coefficient e under stress is introduced to quantitatively explore the influence of effective stress on shale fracture permeability. It is defined as the change in the shale's permeability divided by the initial permeability due to every unit change in effective stress. A larger e value indicates a more sensitive shale permeability. In contrast, a smaller e means a less sensitive permeability. The stress sensitivity of a shale's permeability can be calculated using equation (10): where e is the permeability sensitivity coefficient; Δk is the change in permeability (µm 2 ); and Δ e is the change in the effective stress (MPa). Figure 9 shows the evolution of the permeability sensitivity coefficients of four samples under different effective stresses. The permeability sensitivity coefficient ranges from 0.0015 to 0.337 MPa −1 . The permeability sensitivity coefficient rapidly decreases before 12 MPa and changes less after 12 MPa because of the fracture compaction effect. In addition, the permeability sensitivity coefficient fluctuates with changing effective stress, which may be related to the internal contact state and the spatial structure of the fractures.
The dependence of the shale fracture permeability on stress controls the final gas production performance. Walsh et al 33 established a model of permeability and stress, which is based on the stress dependence of permeability (k) generated when two random rough fractures come into contact. If the surface roughness of the fracture does not change with the effective stress, the model can be simplified as follows: 20,34 where h represents the scalar of the roughness height distribution corresponding to half of the fracture aperture value under the initial stress (µm), i represents the stress under i MPa, 1 represents the stress under 1 MPa of initial stress, and is the fitting parameter. Figure 10 shows that the relationship between (k i /k 1 ) (1/3) and ln(σ i /σ 1 ) is approximately linear. The fitting parameters are shown in Table 4. The fitting accuracy is high, which indicates that the relationship between the permeability of the shale fracture and stress obeys the Walsh theory. Since the initial aperture is considered to be a constant, the change in the effective stress will change the slope of the curve. There Relationships between the permeability damage rates and the effective stress are some inflection points in the curves, for example, an inflection point appears at 12 MPa in a1, 8 MPa in a2, and 6 MPa in b1, while b2 almost has no inflection point. The appearance of inflection points indicates that part of the rough asperity has been destroyed, which changes the permeability of the fracture.

| Friction factor of the gas flow in the fracture
The aperture between two fractures changes, and the gas flow follows a tortuous path within the fracture. The friction between the gas and the fracture causes energy loss, which is called frictional resistance. The frictional resistance can be used to characterize the head loss in a rock fracture, and it is usually expressed by the Darcy Weisbach formula: 35 Because of the strong compressibility of gas, the friction factor cannot be directly calculated using equation (12), but when combined with equation (3), we can derive equation (14) for calculating the friction factor of the gas in a shale fracture: where f is the friction factor; v is the average flow velocity through a cross section perpendicular to the flow direction (m/s); and D h is the hydraulic diameter (twice the hydraulic aperture) (µm). CO 2 gas has a low viscosity and can diffuse through fractured rock more easily, so when CO 2 gas flows within a fracture, it is more likely to produce nonlinear flow. The Reynolds number (Re) is the ratio of inertial force to the viscous force (equation 15), making it an important parameter for judging the state of the fluid flow. It is widely used to distinguish between linear (laminar) and nonlinear characteristics of fluid flow. The fluid flow in shale fractures is very complex. When the fracture aperture is small and the velocity is low, the viscous force is dominant, and the inertial force can be ignored. When the velocity is higher in relatively large apertures, the inertial force is greater and the fluid flows irregularly or turbulently within the fracture. At this time, Darcy's law is no longer applicable, and the nonlinear phenomenon of fluid flow in the fracture is obvious: Non-Darcy flow is complex because it involves both viscous resistance and inertial resistance. Generally, the relationship between the Reynolds number, Re, and the friction factor, f , is used to determine the flow state. Figure 11 shows the variation in the friction factor with Reynolds number under different effective stresses. Their variation law is similar to that of the classical Moody's diagram. 36 All of the points fall on the upper side of the parallel plate model curve because as the effective stress increases, the fracture becomes relatively rough and the flow path becomes more tortuous, so the friction factor increases. The friction factor decreases with increasing Reynolds number, and the nonlinearity of the gas flow in rough fractures becomes more apparent at higher Reynolds numbers. When the Reynolds number is small, the Re-f curve is linear, and the state of the fluid flow is less affected by the effective stress. Therefore, the frictional resistance is mainly affected by the Reynolds number. In this case, the curve is close to the parallel plate model line, which indicates that the viscous resistance has a greater impact on the fluid flow, while the inertial force has a smaller impact. As the Reynolds number increases, the data deviate from the parallel plate model line, the curve becomes horizontal, and the fluid gradually changes from laminar flow to transitional flow, which indicates that the inertial force begins to play a major role in the fluid flow. At this time, the state of the fluid flow is not only affected by the Reynolds number, but it is also controlled by the effective stress. As the effective stress increases, the friction resistance coefficient increases, and the nonlinear flow appears earlier. The reason for is that with increasing effective stress, the contact area increases, and the fracture's relative roughness begins to play a leading role. Thus, the flow velocity and direction are affected by the fracture contact, which is more likely to lead to the occurrence of nonlinear flow. This phenomenon is similar to the result of Zhou et al 25 that with increasing relative roughness, the fluid approaches a state of Forchheimer flow and complete turbulence. There is no obvious boundary between the laminar and transition flows. The experimental curve does not reflect the entire process from laminar flow to partial turbulence and then to complete turbulence because of the limitation of the test data points. The relative roughness of the fracture surface and the frictional resistance of the fluid in the fracture change with increasing effective stress. According to our dimensional analysis, 36 the frictional resistance is related to the Reynolds number and the relative roughness. Many scholars have proposed friction resistance models, but these models have some limitations. For example, in the laminar flow state with small Reynolds numbers, the friction resistance coefficient is overestimated. 26 The fluid flow path is very tortuous because it is influenced by the irregular fracture geometry. The nonlinear coefficient, , in the Forchheimer equation is closely related to the fracture geometry and changes with effective stress. Based on the above analysis, we propose a friction model based on the Reynolds number and the nonlinear coefficient (equation 16). The first term on the right side of the equation   describes the linear flow, and the second term describes the nonlinear effect. Figure 12 shows the data and the fitted surfaces of the nonlinear friction model. The correlation coefficient is about 0.989: where a f and b f are the fitting parameters.

| Flow regime and critical Reynolds number analysis
In Darcy's law, the hydraulic transmissivity T is the product of the fluid permeability and the flow area of the fractured rock, and it is an important parameter used to characterize the seepage capacity. 37 In order to explain the nonlinear flow characteristics in rock fractures, some researchers 5,38 have used the normalized transmissivity (T/T 0 ) to reflect the flow capacity of rough rock fractures. T/T 0 (equation 17) is defined as the ratio of the apparent hydraulic transmissivity (T) to the inherent hydraulic transmissivity (T 0 ). T 0 describes the hydraulic transmissivity when the fluid velocity is very low and the inertial force can be ignored. It is worth noting that T/T 0 = 0.9, indicating that the nonlinear term accounts for 10% of the total pressure drop, and the nonlinear behavior of the fluid cannot be ignored. When T/T 0 < 0.1, which means the head loss caused by the inertial term accounts for more than 90% of the total head loss, the fluid enters the stage of complete turbulence: 25 Most of the T/T 0 values calculated in this study are >0.1. When 0.1 and 0.9 of the normalized transmissivity T/T 0 are used to judge the fluid flow state, the gas flow in the shale  fractures is generally in the transitional flow region, laminar flow seldom occurs (Figure 13), and complete turbulence cannot be observed, which is very consistent with the results obtained from the Re-f diagram.
To determine the start of nonlinear fluid flow (critical Reynolds number Re c ), we plotted the normalized transmissivity T/T 0 as a function of the Reynolds number Re (equation 18). 5,39 The fitting result is shown in Table 5. All of the experimental data points fall on the fitting curve, and the correlation coefficient is 100%. Zimmerman et al 32 carried out a water flow test through a rough fracture, thereby obtaining the relationship between the Reynolds number and the normalized transmissivity (T/T 0 ). They determined that the nonlinearity increased rapidly when the Reynolds number is >10. In this study, by fitting the experimental data, the range of all of the T/T 0 values was plotted in Figure 13: where is the fitting parameter. Figure 13 shows that the CO 2 flow within the fracture mainly undergoes three stages: the viscous stage, the weak inertia stage, and the strong inertia stage. For the normalized hydraulic transmissivity coefficient, when Re is small, the T/T 0 -Re curve is a horizontal straight line (slope close to 0) because the viscosity effect is dominant. As Re increases, the T/T 0 -Re curve begins to bend downward because the flow of CO 2 in the fracture is affected by both viscosity and inertia. When Re increases to a certain value, the T/T 0 -Re curve begins to decrease linearly. At this time, the inertial effect is dominant. The changes in the T/T 0 -Re curve clearly illustrate the flow process within the fracture as it evolves from the linear flow state to the nonlinear state. This trend is in good agreement with the results of Zimmerman et al. 39 When T/T 0 = 0.9 is adopted to determine the starting point of the transition from laminar flow to nonlinear flow, we can use equation (18) to calculate the critical Reynolds number Re c. The calculation results are shown in Figure 14. The critical Reynolds numbers of the four shale samples were found to be 3.07-78.25, 7.98-102.88, 1.8-34.29, and 2.88-46.1. The critical Reynolds number decreases with increasing effective stress, and the trend of the critical Reynolds number is very close to that of the Re-f curve in Figure 11. When the effective stress increases from 1 to 30 MPa, the critical Reynolds numbers of the four shale samples decrease by 96.08%, 92.24%, 94.75%, and 93.75%. The hydraulic transmissivity of the fracture is closely related to the hydraulic aperture, and small changes in the fracture aperture lead to large changes in the flow velocity and Re. Under low effective stress, the hydraulic aperture of the fracture is larger, and the surface roughness has less influence on the fracture flow, nonlinear flow does not occur easily, and the critical Reynolds number is large. As the effective stress increases, the aperture becomes smaller, the flow path becomes more tortuous, and thus, the nonlinear flow is more likely to occur, and the critical Reynolds number becomes smaller. In addition, the Re c ranges of the four shale samples are different because the nonlinear range in the rock fracture is related to factors such as surface roughness and effective stress. 40

| CONCLUSIONS
In this study, the nonlinear flow of CO 2 through rough shale fractures under different effective stress conditions was studied, and the change in permeability with effective stress was analyzed. The following conclusions are drawn: 1. There is a nonlinear relationship between the volumetric gas flow rate and the pressure gradient in shale fractures. The accuracy of the fitting of the experimental data based on the Forchheimer equation, which takes into consideration the high compressibility of gases, provides a new method for the analysis of gas flow in fractured rocks. 2. Shale fractures are very sensitive to effective stress. Their permeability decreases exponentially with increasing effective stress, and the initial stress greatly damages their permeability. The permeability measured parallel or perpendicular to the bedding plane is similar, and it is still