1. Introduction
With the increasing depletion of hydrocarbon resources, unconventional resources such as shale, tight sandstone, and other reservoirs, which were previously considered to be undevelopable, have become hot issues in recent decades [
1,
2,
3]. Compared to conventional low permeability reservoirs, shale reservoirs possess extra-low permeability, usually among nanoscale range, making the development methods for conventional low permeability reservoirs unsuitable for shale reservoirs development [
4,
5,
6]. In recent years, a variety of methods were attempted for the development of shale reservoirs [
7,
8,
9,
10]; however, there is still no final conclusion about which method is more suitable. However, there is a consensus that although the shale reservoirs hav extremely low permeability, they posess the common characteristics of developed bedding and natural fracture as well as a high content of brittle minerals. Thus, large-scale volume fracturing technology is usually recommended to transform shale reservoirs [
11], forming a volumetric fracture network structure in three-dimensional space, which could communicate the matrix of shale reservoirs and the wellbore, realizing the effective migration of hydrocarbon from the matrix to the wellbore. Therefore, the volumetric fracture network structure is the key to the effective development of shale reservoirs. The volumetric fracture network structure is composed of crisscrossed natural fractures, induced fractures, and hydraulic fractures. According to the fracture morphology, the aforementioned three kinds of fractures can be divided into tensile fractures and shear-slipping fractures [
11]. Regardless of the nature and morphology of fractures, the volumetric fracture network structure is essentially composed of fractures, although the fractures are of different shapes.
As the only path for fluid migration from the reservoir matrix to the wellbore, the fracture that directly contacts with reservoir fluid has a binding effect on the migration of fluid [
12,
13,
14]. The characteristics of fracture surface morphology have a direct impact on fluid migration, which requires an accurate and reasonable evaluation of fracture surface morphology. However, in the early studies of fluid migration in fractures, the influence of fracture surface morphology was generally ignored, and the fracture was usually equivalent to a smooth parallel plate model [
15,
16,
17,
18]. The equivalent parallel plate model was used for a long time since its’ proposal. However, with the deepening of understanding, it was found that the equivalent model had a great influence on the research results in specific applications. Therefore, in order to understand the characteristics of fluid migration more accurately, the characteristics of the fracture surface morphology cannot be ignored. As for the fluid migration in the fracture [
19,
20], proppant placement in the fracture [
21,
22], groundwater seepage [
23], and other related situations [
24], the actual results that considered the characteristics of the fracture surface morphology were indeed quite different from those of ignoring the morphology characteristics [
19,
20,
21,
22,
23,
24], which also verified that the characteristics of the fracture surface morphology had a significant influence on fluid migration in the fracture. Therefore, various parameters were proposed to describe the characteristics of fracture surface morphology. However, only one single parameter was usually selected for the characteristics description of fracture surface morphology in specific applications [
19,
20,
21,
22,
23,
24].
At present, few reports on the systematic characteristics description of fracture surface morphology are available in the literature, while no reports about the characteristics description of shale fracture surface morphology were found. In order to better serve the development of shale reservoirs, the shale samples from Barnett Shale was taken as the example. The characteristics of shale fracture surface morphology are described by parameters such as roughness, joint roughness coefficient, fractal dimension, tortuosity, and dip angle, so as to comprehensively describe the characteristics of shale fracture surface morphology and provide a solid foundation for accurately describing fluid migration in shale fractures. In addition, this research will also provide a reference for fracture description in other unconventional reservoirs.
2. Acquisition of Shale Fracture Surface Height Distribution
Three shale samples collected from Barnett Shale were fractured using the Brazilian test, and the fractures of the three shale samples are shown in
Figure 1. The height distribution of the shale fracture surface was obtained by a three-dimensional optical profilometer, and both the lateral and longitudinal interval for data collection was 0.01 mm. The principle of data acquisition by a three-dimensional optical profilometer was listed as follows: First, a lateral coordinate was fixed, and the longitudinal height distribution under the lateral coordinate was collected successively with the data sampling accuracy taken as the interval until the data acquisition under this lateral coordinate was completed. Then, the lateral coordinate was moved by one longitudinal interval and fixed. The longitudinal height distribution under this lateral coordinate was also collected successively until all the information along this lateral coordinate was collected. The aforementioned steps were repeated until all the height distribution information of the shale fracture surface was collected.
According to its scanning principle, the three-dimensional optical profilometer transformed the morphology of the shale fracture surface into the point set of the profile. However, it could not obtain all the information of the shale fracture surface, as the scanned part of the adjacent profile was the blind area where the information could not be collected. Therefore, the blind area was excluded when calculating the relevant description parameters; in other words, the description parameters of shale fracture surface morphology were calculated according to the information from the profile obtained using a three-dimensional optical profilometer. The specific data acquisition process is shown in
Figure 2.
3. Description Parameters of Shale Fracture Surface Morphology and Their Calculation
3.1. Height Distribution Pre-Processing of Shale Fracture Surface Morphology
As the fracture was split by the Brazilian test, in order to prevent the influence of end discontinuity and data acquisition, the central position of the shale fracture surface was firstly located. The rectangular region with 35 mm length and 21 mm width that lay in the central area of shale fracture surface was selected as the object for characteristics description, and the central position was used as the center of the rectangular region.
Prior to calculating the description parameters of shale fracture surface morphology, the height distribution of the shale fracture surface was processed as follows: for the data point set to be studied, the initial point of both lateral direction and longitudinal direction were reset to zero; then, the minimum height of shale fracture surface was regarded as the zero point of height distribution, and the minimum height was subtracted from each data point to obtain a new set of height distribution values for the shale fracture surface. The study focused on the characteristics of shale fracture surface morphology, which was only related to the fracture morphology and was independent of the location and the absolute height of the point. Therefore, the data processing had no effect on the results.
The fracture surface morphologies of three shale samples were reconstructed after the aforementioned pre-processing, and the reconstructed images are shown in
Figure 3.
3.2. Characteristics Description of Shale Fracture Surface Morphology
So far, there are no systematic reports available regarding the characteristics of fracture surface morphology, and the characteristics description methods for fracture surface morphology are randomly selected. Generally speaking, roughness and fractal dimension are widely accepted to describe the characteristics of fracture surface morphology in specific studies that involve fluid migration, while joint roughness coefficient is a classical parameter recommended in the field of rock mechanics. However, tortuosity is often adopted when the study is related to both the characteristics description of fracture surface morphology and fluid migration. Besides, dip angle is a parameter used to describe the overall trend of fracture surface morphology. There are few reports that utilize dip angle alone to describe the characteristics of fracture surface morphology.
One single parameter was mostly used to describe the characteristics of fracture surface morphology, but the selection basis was not expounded [
19,
20,
21,
22,
23,
24]. Most of the selection was based on experience; however, it was doubtful whether one single parameter could comprehensively describe the characteristics of fracture surface morphology. Considering the disadvantage of the single description parameter used currently, the characteristics description parameters of fracture surface morphology were summarized, and their calculation methods were reasonably selected according to the scanning principle of the three-dimensional optical profilometer to systematically and comprehensively calculate the characteristics description parameters of shale fracture surface morphology.
3.2.1. Roughness
Roughness is not only the most basic parameter to describe the characteristics of fracture surface morphology but also the most widely verified parameter [
25]. During the study of nonlinear flow through rock fracture networks, it was shown that the surface roughness that contributes to complex void geometries and streamlines structures could significantly reduce the critical
Re for the flow regime transition [
25]. Roughness belongs to microscopic geometrical error, and it refers to the asperity of the surface morphology with both small spacing and small peaks/valleys. The smaller the roughness, the smoother the surface morphology will be. The commonly recommended parameters for calculating roughness include standard deviation, skewness, kurtosis, arithmetic average height, the maximum height of the peak, etc. [
26].
As for the description of fracture surface morphology, standard deviation is usually employed to calculate fracture roughness, which can be divided into overall standard deviation and sample standard deviation [
27]. Taking the principle of data acquisition into account, the roughness of profiles along both lateral and longitudinal directions were calculated, and the arithmetic square root of the two was regarded as the roughness of shale fracture surface morphology. Sample standard deviation was introduced to characterize the roughness of shale fracture surface morphology, and the specific methods are shown in Equations (1)–(3):
The roughness of shale fracture surface morphology is shown in
Table 1.
As shown in
Table 1, the roughness of the shale fracture morphology was different, among 0.0834 mm–0.2427 mm. For the same shale fracture surface morphology, roughness varied greatly along different directions, with a maximum difference of seven times. Therefore, there was heterogeneity of roughness along different strikes.
3.2.2. Joint Roughness Coefficient
The joint roughness coefficient is a parameter recommended for characterizing the joint morphology in the field of rock mechanics [
28]. Joint refers to those with fractures and no obvious relative displacement on either side of the fracture surface, e.g., with an essential class of fracture structure. However, joints and fractures are not easily distinguished in the specific application of joint roughness coefficient [
29]. In the experimental study of fluid flow through fractured granite, a joint roughness coefficient was employed to qualitatively describe the joint, and specimens with fairly smooth joints subjected to low confining pressures showed a linear laminar flow region against fluid pressures [
30]. Barton firstly proposed the concept of joint roughness coefficient (
JRC) through the summarization of rough fracture profiles [
31], and ten standard profiles (
JRC values among 0–20) were proposed. The actual value of the joint roughness coefficient could be acquired by comparing the actual profile of fracture surface to that of the standard joint roughness coefficient. This method has been widely applied in related descriptions of fracture surface morphology since it was proposed and recommended, and many empirical equations used for the calculation of fracture surface parameters were developed according to their relationship to fit the joint roughness coefficient [
32,
33,
34,
35].
The most widely used statistical parameter
Z2 was adopted to calculate the joint roughness coefficient of fracture profiles [
33], while the
Z2 along the lateral direction was listed as follows,
The equations of the relationship between
Z2 and joint roughness coefficient developed by Tse and Cruden were chosen to calculate the joint roughness coefficient of the profile along the lateral direction [
32], and it was formulated as,
Both the calculations of the statistical parameter
Z2y and the joint roughness coefficient
JRC3Dy along the longitudinal direction were similar to that of the lateral direction. The joint roughness coefficient of shale fracture surface morphology is shown in Equation (6).
The
JRC3D of shale fracture surface morphology is shown in
Table 2.
As shown in
Table 2, the joint roughness coefficient of different samples varied greatly, among 2.5715–10.9368. The difference between the joint roughness coefficient along the lateral direction and that along the longitudinal direction was significant, up to two to three times.
3.2.3. Fractal Dimension
The concept of fractal dimension was proposed based on the self-similarity of geometric features, which is an important metric parameter to determine the disorder of geometrical morphology and the irregularity of complex geometries [
36]. Compared to the integer dimension reflecting the static characteristics of the object, the fractal dimension represents the dynamic change process of the object. If it is extended to the dynamic behaviors and phenomena of nature, the fractal dimension is a correlated characterization of the whole system that is composed of small local characteristics in natural phenomena. That is, for an object, it can accurately reflect its irregularity and complexity using the dimension scale of the measured non-integer value; as a result, the dimension of the non-integer value is regarded as the fractal dimension [
37]. During the experimental analysis of single-phase flow through rough fracture replicas, the compensation of turbulence effect developed earlier for rough fractures was compared to the smooth parallel plate model and was longer for the fractures with higher fractal dimensions. Transmissivity values showed a decrease with increasing fractal dimensions and normal loading and also exhibited directional-dependent behavior. The percentages of water-invaded wet planar areas showed a tendency to decrease with increasing fractal dimension [
38].
For irregular profiles of complex fracture surfaces, the yard stick method and covering method are commonly used [
39]. However, for fracture surfaces, it is impossible to directly cover the rough surface with two-dimensional Euclidean geometry with a certain scale.
Currently, the reported methods for calculating the fractal dimension of fracture surface include the triangular prism surface area method, projection covering method, and cube covering method [
39]. The first two methods possess the defect that the calculation of rough surface area is approximate, which leads to the deviation of calculation results. In contrast, the cube covering method can effectively avoid the calculation deviation caused by approximation. Besides, the fractal dimension estimated by the cube covering method is the pure geometric fractal dimension. Each calculation step has an accurate method, and there is no approximate process during the calculation; therefore, the calculated fractal dimension is close to the true value.
The specific process of the cube covering method was provided as follows [
40]: assuming there is a square grid on the plane
XOY, and the length of its basic element unit is
δ, the four corners of the square grid correspond to the four heights
h(
i,
j),
h(
i,
j + 1),
h(
i + 1,
j), and
h(
i + 1,
j + 1) (
i,
j is among (1,
n − 1), and
n is the number of measuring points each side), respectively. Cubes with side length
δ were employed to cover the fracture surface, and the number of cubes in the covering area
δ ×
δ was calculated, that is, the number of cubes
Ni,j covering the fracture surface in the grid (
i,
j) was provided as follows,
Therefore, the total number of cubes required to cover the entire fracture surface was
Then, the whole fracture surface was covered again by changing the measuring dimension, and the total number of cubes for covering the whole surface was calculated using Equation (8). According to the fractal theory, the following relationship could be obtained,
Moreover, the logarithms of N and δ at both ends of Equation (9) were, respectively, calculated. The relationship between them was fitted and plotted, and the slope was fractal dimension D.
The fractal dimension of shale fracture surface morphology is shown in
Table 3.
3.2.4. Tortuosity
The concept of tortuosity was primarily proposed in the study of fluid seepage in porous media, and it is defined as the ratio of the fluid actual seepage path to its apparent length. Tortuosity represents the complexity of the fluid seepage path, so it essentially describes the complexity of the fluid seepage channel, which is also known as the channel tortuosity. Tortuosity has a negative effect on the permeability of the porous media, which means that the fluid seepage resistance increases with the increase of tortuosity [
41]. Additionally, tortuosity is closely related to the physical properties of porous media, such as porosity and permeability [
41].
The methods to study tortuosity mainly include casting thin sections, a capillary pressure curve, and empirical formula [
42,
43]. Currently, there are two main methods used to calculate the tortuosity of fracture surface morphology, which are, respectively, based on the profile and area of fracture surface morphology. Based on the acquisition principle of fracture surface height distribution, the calculation based on the profile of fracture surface morphology was selected to calculate the tortuosity of shale fracture surface morphology. For the characteristics of shale fracture surface morphology, three-dimensional distribution, both the lateral and longitudinal tortuosity along fluid seepage direction, were calculated, respectively, and the root mean square of the two tortuosities was taken as the tortuosity of shale fracture surface morphology. The specific equations are provided as follows:
The tortuosity of shale fracture surface morphology is shown in
Table 4.
3.2.5. Dip Angle
The concept of dip angle was proposed during the description of fracture structure, and then dip angle was gradually adopted for the description of fracture structure. Dip angle is a parameter used to describe the trend of the fracture surface; however, the application of dip angle in fracture description is less than that of other parameters. There is almost no report of taking dip angle alone to describe the characteristics of fracture surface morphology, and it is generally used in combination with other parameters. During the study of gas flow in shale micro-fracture, dip angle combined with other parameters were employed to fit shale gas flow, and the results showed that dip angle had an inhibitory effect on gas flow in shale micro-fracture [
44].
The angle of the profile can be divided into the angle with an upward trend and the angle with a downward trend along the profile; the influence of the two on fluid seepage is significantly different. The angle of the upward trend has an inhibitory effect on fluid seepage, while the angle of the downward trend has a positive effect. Therefore, the dip angle is divided into absolute dip angle and overall trend dip angle. The absolute dip angle refers to the absolute value accumulation of the angle along the profile, and upward and downward trends are not distinguished. The absolute angle represents the complexity of the profile. The overall trend dip angle refers to the superposition of the angles along the profile, distinguishing the angle of the upward trend from the angle of the downward trend. The overall trend dip angle represents the overall trend of the profile.
Based on the scanning principle of a three-dimensional optical profilometer, the dip angle of each profile was calculated successively, and the dip angle of the whole shale fracture surface morphology could be obtained by weighted averaging the dip angles of all profiles.
The absolute dip angle was calculated by Equation (13), and the arithmetic square root of the absolute dip angle product in both directions was taken as the absolute dip angle of shale fracture surface morphology. The details are shown in
Table 5.
As can be seen in
Table 5, the absolute dip angle of shale fracture surface morphology showed different characteristics along both the lateral and longitudinal directions. The absolute dip angle of shale sample 1 was basically the same along both directions, which were all around 24°. The absolute dip angle of shale sample 2 along the lateral direction was about 3° higher than that along the longitudinal direction, which was 19.406° and 16.235°, respectively. The absolute dip angle of shale sample 3 along the lateral direction was about 10° smaller than that along the longitudinal direction, which was 19.539° and 29.861°, respectively. In general, the absolute dip angle along the lateral direction was 19.406°–24.256°, and the absolute angle along the longitudinal direction was 16.235°–29.861°. The difference between shale samples along the lateral direction was about 5°, while it was about 13° along the longitudinal direction.
Meanwhile, the overall trend dip angle was calculated using Equation (14). The arithmetic square root of the overall trend dip angle product in both directions was taken as the overall trend dip angle of shale fracture surface morphology. The details are shown in
Table 6.
The overall trend dip angle of shale fracture surface morphology is shown in
Table 6. The overall trend dip angle of different shale samples varied greatly in both the lateral and longitudinal direction, and that of shale sample 3 was minimal, about 2.4°, while that of shale sample 1 was the maximum, about 12°. The difference between the overall trend dip angle of shale samples along the lateral direction was smaller, about 2.5°, while that along the longitudinal direction was larger, about 22.5°, roughly nine times larger than that along the lateral direction.