Influence factors of unpropped fracture conductivity of shale

Hydraulic fracture is widely performed to enhance oil and gas production in the development of unconventional reservoirs. This paper proposes a numerical model for the calculation of the unpropped fracture conductivity of shale. This numerical model consists of three basic models, namely, the numerical reconstruction model of fracture morphology, fracture deformation model under closed stress, and fluid flow model. The numerical model can be used to analyze the influencing factors of the unpropped fracture conductivity of shale quantitatively. The results demonstrate that the fracture misalignment of shale is a prerequisite for unpropped fractures with conductivity at high closed stress. Furthermore, the unpropped fracture conductivity of shale decreases exponentially with the increase in closed stress. Moreover, closed stress is the factor having the most significant influence on the unpropped fracture conductivity of shale, followed by roughness and fracture morphology. The negative effect of closed stress can be reduced by increasing the fracture morphology during the fracturing operation.


| INTRODUCTION
"Horizontal well + segmented multi-cluster perforation + large-displacement hydraulic fracturing" has led to the commercial development of unconventional reservoirs, especially of shale gas. 23,33,35 The main and branch fractures formed during fracturing are the tensile fractures of the shale reservoir, whereas the secondary fractures are the shear misalignment fractures,all the fracture surfaces are rough. 30,39 When the fracturing is completed, the fractures with proppant cannot be closed, and the unsupported fractures are relatively dislocated owing to the geostress and cannot be completely closed, thereby resulting in unpropped fractures. 34,38 An unpropped fracture with conductivity is the basis of fluid flow in fracture networks (Appendix A-D).
Fracture conductivity refers to the ability of a fracture to fill a proppant to pass through a fluid under reservoir stress, which can be expressed by the fracture width multiplied by the fracture permeability. So far, the studies on unpropped fracture conductivity are primarily based on experiments, and the mathematical models are empirical formulas fitted according to the experimental results. For instance, Penny 20 determined through experiments that time and temperature are the main factors influencing the decline in long-term conductivity. Van and De 28 demonstrated experimentally that fractures without proppant still have conductivity even under closed stress. Numerous experiments that demonstrate that the fracture conductivity decreases exponentially with the increase in closed stress have been conducted. 9 Although experiments can simulate the formation state to obtain more accurate conductivity, there are several factors that influence the unpropped fracture conductivity of shale. Hence, it is difficult to conduct these experiments owing to the following reasons: (a) According to the American Petroleum Institute (API) standard, the shale sample is large in size (length: 178 mm, width: 38 mm), which makes processing difficult. (b) The experiment is time-consuming owing to the close correlation between the experimental steps. (c) There are several factors that affect the conductivity of shale (such as fluid properties, temperature, stress, initial fracture surface morphology, rock mechanical properties, and misalignment movement), making the analysis of the data difficult.
"Acid etching fracture conductivity" is another important method that can be used to study fracture conductivity. Nierode and Kruk 19 proposed the Nierode-Kruk model for the quantitative evaluation of acid etching fracture conductivity. Based on this model, researchers have developed more accurate conductivity models (Table 1) considering factors such as closed stress, rock mechanical properties, and fracture surface roughness. Table 1 indicates that the theoretical model of acid etching fracture conductivity is directly or indirectly related to the roughness and morphology of the fracture surface. Owing to the large difference in surface roughness between the acid etching fracture and unpropped fracture, the application of these models to calculate the unpropped fracture conductivity is limited. Therefore, the morphology and roughness of a fracture are vital in calculating the fracture conductivity of shale.
The fracture morphology comprises two parts: wall roughness and opened fracture surface. 1,24 Tse et al 26 and Turk et al,27 proposed to reconstruct the rough-walled fracture surface using an interpolation method. However, this method is prone to errors due to frequent interpolation iterations. 32 Lespinasse and Sausse 15 used an automatic surface profilometer to measure the fracture surface roughness and consequently proposed a set of methods to describe the distribution of a rough-walled fracture. Issa

Tsang and
Witherspoon 25 observed that fracture wall surfaces generated by natural faults exhibit fractal characteristics. Therefore, the fractal dimension can be used to describe the fracture wall roughness. Brandt and Prokopski 2 observed that the fracture roughness is proportional to the fractal dimension. Barton et al 1 and Barton 43 observed that the distribution of fracture width obeys the lognormal and negative exponential distributions. The unpropped fracture conductivity of shale is also affected by vertical stress. 13 Studies have shown that loading effective stresses on the fractures can significantly change the fracture width and the fluid velocity and pressure distribution within the fracture. 4,5,40 A rough-walled fracture can be deformed by vertical stress, thereby reducing the fluid flow ability of the fracture. 22 Subsequently, Louis 17 proposed that vertical stress can be more accurately expressed using "effective stress." He established an exponential relationship between effective stress and flow coefficient. 17 Scholars have always hoped to describe the relationship between stress and fluid flow in fractures quantitatively through a mathematical model,however, most models are empirical formulas obtained via fit and regression according to the experimental data.
Unpropped fractures of shale have characteristics of large size, high roughness, and random morphology. 21 Thus, the Navier-Stokes equation is often used to study the fluid flow in fractures. 8 To simplify the calculation, the inertial force is often neglected, and the nonlinear three-dimensional (3D) Navier-Stokes equation is reduced to a linear 3D model. 42 However, when the Reynolds number is too large, the inertial force cannot be neglected, and the nonlinear 3D Navier-Stokes equation cannot be simplified to ensure calculation accuracy. 14 The 3D model is approximated to a two-dimensional model by using the "regular discrete fracture morphology model" under uniform fracture width. 16 When the rough-walled fracture length is greater than five times the fracture width, the error between the simplified model and the Navier-Stokes equation is >10%. 41 The Navier-Stokes equation has considerable limitations in the study of rough-walled fracture morphology. The Lattice Boltzmann (LBM) discrete numerical simulation method can effectively match the geometry of rough-walled fractures and examine the law of fluid motion. 31,37 The LBM method discretizes fluid particle motion comprising two main steps: migration and collision. 36 Migration is the conversion of particles from one node to the next on the grid. Collisions preserve the momentum by redirecting "collisions" or particles that occupy the same node. Martys and Hagedorn 18 improved the LBM to make it easier to describe node velocity and fluid momentum through probability distribution functions.
In this paper, based on the fracture shear slip and characteristics of rough-walled fractures, we proposed a rough-walled fracture numerical restructure model. Additionally, we proposed a fracture deformation model under closed stress based on the mechanical characteristics of shale. We adopted the 3D LBM method to establish the fluid flow model in rough-walled fractures. We constructed a mathematical model to calculate the unpropped fracture conductivity of shale using the joint fracture numerical restructure model, fracture deformation model under closed stress, and fluid flow model. The shale used in this study was obtained from Longmaxi Formation in the south of the Sichuan Province to verify the correctness of the numerical model.

| MATHEMATICAL MODEL
The mathematical model used for the calculation of the unpropped fracture conductivity consists of a numerical reconstruction model of fracture (Appendix A), fracture deformation model under closed stress (Appendix B), and fluid flow model (Appendix C), as illustrated in Figure 1.

| Numerical reconstruction model of fracture
The height matrices of the upper and lower fracture surfaces after misalignment are expressed as Equations 1 and 2, respectively.
The height division for the upper surface of the fracture is expressed by Equation 1: The height division for the lower fracture surface is expressed by Equation 2: where A represents the initial height matrices of fracture obtained via laser scanning; Z is the height matrix of fracture after misalignment; a is the element of the matrix;A ′ 0R is the vertical-reverse height division of the upper fracture surface;A ′ 0S is the symmetric placement matrix; the subscripts 0, 1 and u, d represent the upper and lower surfaces of the fracture, respectively; SDIS indicates a fracture having misalignment movement. (1) LU et aL.

| Fracture deformation model under closed stress
The morphology of the unpropped fracture under closed stress can be expressed by Equations 3 and 4.
Here, Z 1 and Z 2 represent the microelement height matrices of the upper and lower surfaces of the fracture, respectively;Z ′ u represents the height division after the upper surface of the fracture moves down; Z c represents the height division after fracture deformation; ΔZ represents the compression deformation matrix.
By substituting Equation 4 into 5, the relationship ( = f (Z) ) between the closed stress and fracture morphology can be obtained via cyclic calculation as follows: where E, , and m represent Young's modulus, Poisson's ratio, and the compressive strength of shale, respectively.
If the closed stress on shale is known, then the fracture width can be expressed by Equation 6.
where W f is the fracture width, m.

| Fluid flow model
According to the height division of the lower fracture surface, Z d , under closed stress and deformation, and the corresponding fracture width distribution matrix, W f , a physical model for the LBM flow field calculation was constructed. The meshing of the LBM flow field consists of regular squares. The "0 and 1 flag" is used to set the solid-liquid boundary of the fracture space structure. In this study, the LBM calculation method in Bhatnagar-Gross-Krook (BGK) format is used to simulate the flow field in unpropped fractures. The D3Q19 single-speed model (d'Humières et al 44 ) is employed, and the equilibrium distribution function of this model is determined as follows: First, the discrete velocity is expressed by Equation 7: The equilibrium distribution function is expressed by Equation 8: where w is the weight coefficient, which is related to the velocity component and can be expressed by Equation 10.
The solid-liquid interface is set as the standard bounce format, and the fluid boundary is set as the gravity-driven periodic boundary. After calculating the equilibrium distribution at each particle density, the particle collision and migration are calculated by Equations 11 and 12, respectively, and a new density distribution function is generated.
The collision calculation function is expressed by Equation 11: The transfer calculation function is expressed as follows: where f i (x,t) and f i (x,t) are the density distribution functions of the fluid after collision and in the current state, respectively; f is the density distribution function of the fluid at equilibrium; is the dimensionless relaxation factor.  Figure 2 depicts the sequence of the process of solving the unpropped fracture conductivity numerical model, which can be described as follows:

| Numerical model solution
• The fracture surface is scanned using a laser to obtain data (A 0 and A 1 ). • The initial unpropped fracture surface morphology is determined based on the fracture surface data and misalignment. • According to the Young's modulus and Poisson's ratio of the shale, the initial unpropped fracture of the shale is subjected to the fracture deformation model, and the stress-misalignment chart is obtained. • The deformation is determined in the stress-misalignment chart, according to the closed stress. • The initial fracture is subjected to deformation calculation to obtain the fracture morphology. • The LBM method is used to mesh and simulate the fractures after deformation, the fracture permeability under closed stress can be subsequently obtained, and the fracture conductivity at any point can be calculated according to API RP 19D-2008.

| Experiment
We evaluated the conductivity of two shale samples (sample 1 and sample 2) selected from the Longmaxi Formation in the Sichuan Basin with different surface roughnesses to verify the accuracy of the numerical model of the unpropped fracture conductivity. The surface of sample 1 is smoother than that of sample 2, and the mechanical properties of the two samples are listed in Table 2.
A "self-supporting fracture test analysis device" 3 was used to test the unpropped fracture conductivity. The experimental steps are as follows: • A shale sample was prepared with a length of 176.0 mm and a width of 36.0 mm. • The sample was divided into two pieces along its natural fracture.
• The upper and lower rock plates of the sample were dislocated by 2.5 mm, and then, they were integrally fixed with tape and sealed with Vaseline. • The rock sample was placed in the test device. A leakage pressure test was performed by applying an injection pressure of 5 MPa and a back pressure of 3 MPa. • The conductivity test chamber was heated at 3°C/min, the closed stress was applied to the rock sample at 6 MPa/ min, and then, nitrogen was slowly injected at 100 mL/ min. • The experimental data were recorded when the closed stress, temperature, and flow rate of nitrogen reached a steady state.

| Numerical reconstruction of fracture morphology
The fracture surface was numerically reconstructed using the upper and lower fracture surface height matrices (Equations 1 and 2, respectively) and the surface data. The results are shown in Figure 3. Sample 1 and Sample 2 are two typical features of shale. It can be observed from the color distribution (Figure 3) that the height of the fracture surface of sample 1 is not as evident as that of sample 2. This indicates that the fracture surface of sample 1 is smoother than that of sample 2. If the fractures are misaligned to form an unpropped fracture, then the unpropped fracture of sample 2 would be more evident than that of sample 1. The distribution of fracture width directly reflects the internal space structure of the fracture. The numerical reconstruction model (Equations 1 and 2) is used to reconstruct the fracture surface with a misalignment of 2.5 mm. The distribution cloud of the misalignment fracture and the distribution frequency of the fracture width are shown in Figure 4. The fracture width distribution frequencies of the two samples demonstrate a single peak, and the distribution of fracture width is concentrated, demonstrating a normal distribution. The misalignment fracture width of sample 1 is mainly concentrated in the range of 1.019-2.719 mm, whereas the fracture width of sample 2 is mainly concentrated in the range of 1.921-5.766 mm. The misalignment fracture of sample 2 provides a more convenient passage for fluids to flow therein.

| Calculation of fracture deformation under closed stress
Based on the fracture deformation model under different closed stresses (Equation 5), we obtained the relationship between the deformation and closed stress of the rock samples, as illustrated in Figure 5. The stress-deformation curve of a fracture has three phases: (a) In the first stage, the deformation increases sharply with the increase in closed stress under lower stresses. (b) In the second stage, as the closed stress increases, the deformation tends to increase linearly. (c) In the third stage, the fracture deformation increases sharply with the increase in closed stress. Combined with the initial fracture width distribution frequency (Figure 4), owing to the narrow distribution of the fractures, the contact points are concentrated, and the microbody rapidly reaches the yield condition and then deforms as the closed stress increases. As the closed stress increases, more convex points come in contact with the rough fracture surface, and the relationship between the fracture deformation and closed stress demonstrates a linear change. More support points reach the yield state with the continuous increase in closed stress owing to the decrease in the number of wide parts, thereby sharpening the fracture width ( Figure 5). The distribution of the fracture under different closed stresses is shown in Figure 6.

| Comparison of calculation and experimental results
The numerical reconstruction model of fracture, fracture deformation model under closed stress, and 3D LBM model are used to calculate the unpropped fracture conductivity of shale with a 0.2-mm grid, and the comparison of the calculation and experimental results is illustrated in Figure 7. The fluid medium used for calculation is nitrogen gas with a viscosity of 17.85 × 10 −3 mPa·s and density of 1.25 kg/m 3 for low-speed flow. As illustrated in Figure 7, the misalignment of the unpropped fracture of shale is a prerequisite for its conductivity under high closed stress (≥4 MPa); the unpropped fracture conductivity decreases exponentially with the increase in the closed stress. This law of unpropped fracture conductivity is consistent with that of sandstone and the previous experimental results for shale. 6,9 The unpropped fracture conductivity obtained using the numerical reconstruction model of fracture, fracture deformation model under closed stress, and LBM model is consistent with the reducing trend of that obtained via experimental tests, as shown in Figure 7. However, under low closed stress, the calculation accuracy is higher than that under high closed stress; furthermore, the calculation accuracy for wide fractures is higher than that for narrow fractures. For example, when the closed pressure is 1 MPa, the errors between the experimental and calculated results of unpropped fracture conductivity of samples 1 and 2 are only 4.51% and 4.17%, respectively. However, when the closed stress increases to 45 MPa, the errors of samples 1 and 2 increase to 40.16% and 22.08%, respectively. This is because, under high closed stress, the width of the partial fracture is smaller than the mesh size (0.2 mm) used for the LBM calculation. When using a 0.2-mm cube mesh for calculation, the wide fracture has a larger number of compute nodes, and hence, the calculation accuracy is higher. When the physical model has a constant size, the calculation accuracy can be improved by reducing the mesh size.

CONDUCTIVIT Y
Several factors affect the unpropped fracture conductivity of shale. To consider the impact of various factors fully, we first use the orthogonal analysis method to design the calculation program; then, we use the conductivity numerical calculation model to calculate the unpropped fracture conductivity of shale; finally, the influence law of each factor is analyzed using the range summation method of the statistics (Appendix D). The higher the bit level of a certain influencing factor, the greater is the unpropped fracture conductivity. A grid of 0.2 mm is selected for low closed stresses (1-4 MPa) and a 0.1-mm grid is selected for high closed stresses (≥4 MPa) to ensure the accuracy and efficiency of the calculation results.

| Fracture surface roughness
The effect of fracture surface roughness on the unpropped fracture conductivity demonstrates an increasing trend first and then a decreasing trend, as shown in Figure 8. When the fracture surface roughness is 2.7-2.8, the range summation is at its maximum. This indicates that, when the fracture surface roughness is 2.7-2.8, the conductivity of fracture is at its highest. When the fracture surface roughness value is too small, the fracture is smooth, and it is not easy to form an unpropped fracture; thus, the range summation is small. However, when the fracture surface roughness value is too high, the fracture width decreases under closed stress, which leads to a decrease in the conductivity. Therefore, the range summation decreases when the roughness value is too high.

| Misalignment
It can be observed from Figure 9 that, when the fracture morphology is less than 4 mm, the range summation increases F I G U R E 4 Distribution of cloud pattern and fracture width after misalignment by 2.5 mm F I G U R E 5 Fracture morphological stress-strain curve with the fracture morphology. When the fracture morphology is 4 mm, the range summation decreases slightly. When the fracture morphology reaches 4 mm, the tops of the support points on the fracture surfaces are in contact with each other, and the unpropped fracture conductivity is at its highest, indicating that the bit-value summation is at its highest. If the fracture morphology is further increased, the range summation will decrease owing to the decrease in the fracture width.

| Closed stress
As can be observed from Figure 10, the range summation decreases with the increase in closed stress. When the closed stress increases from 1 MPa to 10 MPa, the range summation decreases by 307.5/MPa. When the closed stress increases from 10 MPa to 30 MPa, the range summation decreases by 51.6/ MPa. This is consistent with the exponential decrease in the unpropped fracture conductivity with the increase in closed stress.

| Rock mechanical properties
The influence of the mechanical properties of shale such as Young's modulus, Poisson's ratio, and compressive strength on the unpropped fracture summation is shown in Figure 11. It can be observed from the figure that the range summation demonstrates no regularity with the change in the rock mechanical properties. This indicates that the mechanical F I G U R E 1 1 Effect of Young's modulus, Poisson's ratio, and compressive strength of shale properties of shale have an irregular effect on the unpropped fracture conductivity of shale. Moreover, the mechanical properties of shale have a much weaker effect on the conductivity than that of roughness, misalignment, and closed stress.

| Order of influence factors
We use the range (maximum range summation of influencing factor minus the minimum range summation) to evaluate the effect degree of each influencing factor and to identify the key factors affecting the unpropped fracture conductivity of shale. The results are illustrated in Figure  12. The influence of various factors on the unpropped fracture conductivity of shale can be ranked as follows: closed stress > fracture surface roughness > misalignment > compressive strength > Young's modulus > Poisson's ratio. Closed stress has the most significant influence on the unpropped fracture conductivity of shale; thus, the influence of closed stress must be considered during the design of the fracturing process.

| CONCLUSIONS
In this study, a numerical model of unpropped fracture conductivity of shale was proposed, and it was verified through shale fracture conductivity tests. This numerical model was based on the numerical reconstruction model of fracture, fracture deformation model under closed stress, and 3D LBM model. Furthermore, the numerical model of unpropped fracture conductivity was used to analyze the influencing factors of the unpropped fracture conductivity of shale. The main findings are as follows: • Misalignment of shale is a prerequisite for the unpropped fracture with conductivity at high closed stresses (≥4 MPa).
• The influence of various factors on the unpropped fracture conductivity of shale can be ranked as follows: closed stress > fracture surface roughness > misalignment > compressive strength > Young's modulus > Poisson's ratio. • The unpropped fracture conductivity of shale decreases exponentially with the increase in closed stress. • The effect of closed stress can be reduced by increasing the fracture morphology during the fracturing operation.

Numerical reconstruction model building process
The data of the upper and lower surfaces of the fracture are obtained via laser scanning. Then, the data are collated, and the original height matrices of the upper and lower surfaces of the fracture are obtained, defined as A 0 and A 1 , respectively.
The numerical reconstruction of misalignment includes the following steps: removal of the substrate, symmetric placement, misalignment movement, and double-wall single-point contact.

Removal of the substrate
First, the substrate is selected as the lowest point of the fracture surface. Then, all the grid data of the fracture surface are subtracted from the height of the lowest point ( Figure A1).

Symmetric placement
As the upper fracture surface is in the reverse state when laser scanning is performed, it needs to be vertically reversed to obtain the vertical-reverse height division of the upper fracture surface, A ′
The height division of the upper fracture surface is added to the highest value of the lower fracture surface to obtain a plane symmetric placement matrix of the upper fracture surface, and the lower fracture surface remains unchanged, thereby obtaining a double-wall face-symmetric placement map ( Figure A3).

Misalignment movement
The misalignment movement is performed by fixing the lower fracture surface and moving the upper fracture surface. Therefore, the misalignment movement of the fracture has three modes: forward movement, horizontal movement, and inclined movement.
Considering the directions of tensile fracture and geostress, the forward movement mode is often used. The misalignment movement steps are as follows: • The fracture morphology is set to (δ/10) mm. • The fracture morphology δ needs to be added to the matrix ordinate of the upper fracture wall surface, and the matrix abscissa is unchanged. • The abscissa is maintained constant, and point A is defined in Figure A5 (c) as 0 on the ordinate and point B as the maximum value on the ordinate.
After the misalignment movement is completed, the height divisions of the upper and lower fracture surfaces are obtained as A 0SDIS and A 1SDIS , respectively.

Double-wall single-point contact
The fracture minimum width w min after misalignment is subtracted from the mesh data of the fracture upper surface, and the height matrices of the upper and lower fracture surfaces Z u and Z d , respectively, after the contact are obtained.
The original fracture surface is calculated by A1-A6, and the numerical simulation models (A7) and (A8) of the fracture surface are obtained.

Conductivity calculation results for or thogonal analysis
The conductivity of five rock samples was calculated using the established unpropped fracture conductivity model of shale. The results are presented in Table D1. Some actual core samples after the experiment are shown in Figure D1.