Quantitative Study of the Geometrical and Hydraulic Characteristics of a Single Rock Fracture

Abstract

Three-dimensional images of fractured rocks can be acquired by an X-ray micro-CT scanning technique, which allows researchers to investigate the ‘true’ inner void structure of a natural fracture without destroying the core. The 3D fractures in images can be characterised by measuring morphological properties on both fracture apertures and its trend surface, like the medial surface, that reveals the undulation of fracture planes. In a previous paper, we have proposed a novel method to generate fracture models stochastically. Based on a large number of such fracture models, in this work a modified factor was proposed for improving the performance of the cubic law by incorporating the flow-dominant characteristics, including two parameters (aperture roughness and spatial correlation length) for fracture apertures and two (surface undulation coefficient and spatial correlation length) for fracture trend-surface. We assess and validate the modified cubic law by applying it to natural fractures in images that have varying apertures and extremely bended trend-surfaces, with the permeabilities calculated by a Lattice Boltzmann Method as ‘ground truths’.

1. Introduction

Natural fractures are common features in many rocks, occurring at various length scales and at different intensities. Understanding the fluid flow characteristics in natural fractured rock is of importance for petroleum reservoirs [1], geothermal energy development [2,3], nuclear waste disposal [4] and the geologic storage of carbon dioxide [5]. Since intact rocks are often low in permeability, and the fracture network forms the main channels of flow, the effective permeability of a fractured rock is dictated by the permeability of the fracture network, which consists of numerous individual fractures. Therefore, the fluid flow through a single fracture could form a basic building block for understanding the mechanical-hydraulic interactions of natural fractured rock, which has been extensively studied in recent years [6,7,8,9,10,11].

For steady laminar flow between smooth parallel plates, separated by a constant aperture, fracture transmissivity is proportional to the cubic power of fracture aperture, which is well known as the ‘cubic law’ [12]. However, the fractures in natural rocks have complex geometry and widely deviate from the parallel plate model. The cubic law should, therefore, be tuned to incorporate more factors that influence fluid flow [7,13,14].

Earlier theoretical works [15,16,17,18,19,20] usually use two-dimensional aperture fields to represent rough-walled fractures, and many efforts [15,17,19] have been made to assess the effect of variable fracture apertures on fluid flow. For example, Zimmerman [19] modified the cubic law with additional fracture roughness, which is defined as the ratio of aperture standard deviation over mean aperture. Tsang [17] and Moreno [15] analyzed the effect of aperture anisotropy on fracture fluid flow by taking account of the ratio of the spatial correlation length in different directions of aperture field. In addition, the fractal characteristics of fracture void space, consisting of variable apertures at different locations, was considered by Pyrak-Nolte [14] and it was subsequently incorporated into the study of fluid flow through rough-walled fractures. Although these literatures have thoroughly investigated the influence of fracture variable apertures on fluid flow in rough-walled fracture, using aperture field only to characterize a rough-walled fracture over-simplifies the fracture void space to a big extent. It assumes the void space to be perfectly planar, ignoring the natural undulation [20] or waviness in the alignment of the fracture voids. Therefore, in order to analyze the influence of rough fracture characteristics on fluid flow, it is constructive to consider the morphological features of the trend-surface (i.e., medial surface) of fracture profile, which realistically describes the undulation of fracture planes. Similarly, the morphological analysis can also be applied to the pores of rock matrix [21,22].

Advances in imaging technology make it possible to obtain 3D images of fractured rocks. One example of the imaging technology is X-ray computed tomography (CT) [1,23,24,25,26] that can detect the inner structure of nontransparent objects without destroying rock samples. The principle of CT is that different rock components (e.g., minerals, pores, fracture voids, fluids) have different densities, resulting in different X-ray absorption coefficients, which are then converted into different gray values in resultant 2D images that can be used to distinguish the solid and pore space (fracture void) [25] in this work. Simply, they are both precise and direct methods for developing digital fractured rocks, providing a basis for quantitating the geometric features of the fracture void space [27,28].

With 3D digital rocks, there are two major numerical approaches to simulate single- and multi-phase flow processes: direct simulation and pore-network modelling (PNM) [29,30,31,32,33,34,35,36]. Examples of the direct simulation method includes finite difference/element/volume methods [37,38,39] and the Lattice Boltzmann method (LBM) [40,41,42,43,44,45]. However, for the finite difference/element/volume methods, it is a challenging and time consuming issue to deal with the expensive mesh generation and irregular geometries. With regards to those limitations, the LBM has been widely applied to simulate the flow behaviors in 3D digital rock. The advantages of the LBM include simple algorithms, efficient parallelization, and flexibility in implementing complex flow geometries, which makes the LBM profoundly powerful in the simulation of single- and multi-phase fluid flow through porous media. As for the PNM, it is designed to simulate fluid flow on a simplified structure described by a network of pore bodies (nodes) connected by pore throats (bonds). Each network element (i.e., node or bond) is associated with a regular (e.g., a star, triangle, square, or circle) cross section and assigned with a set of geometrical properties (e.g., radii, volume, shape factor, and length), which makes the numerical simulation based on PNM more complex than LBM. Therefore, in this work, we use the LBM simulations of fluid transport with X-ray topographies to estimate fracture permeability.

The objective of this work is to quantitatively determine the collective effects of both fracture apertures and the trend-surfaces on fluid flow in single fractures. To address this, their characteristic parameters are firstly extracted from digital fractured rocks. Fracture structural models are then generated with various parameters about fracture apertures and trend-surfaces. Subsequently, the LBM is applied to simulate fracture flows for predicting fracture permeability, and further a modified version of the cubic law, which incorporates geometric fracture undulation, is proposed and validated by applying it to natural fractured cropped out from 3D μCT images.

2. Methodology

2.1. Fracture Morphology

2.1.1. Natural Fracture in CT Images

A 3D digitized rock, which reveals the pore structure in the rock matrix and fracture opening/void, is made up of a series of 2D images which are obtained by X-ray computed tomography (CT). Figure 1 shows an example of a stack of cross-sectional slice CT images of a shale core, which includes fracture structure (e.g., red box in Figure 1a) and were provided by iRock Technologies (http://www.irocktech.cn/). These 1024 successive 2D shale images in Figure 1a were aligned and stacked to form a 3D image (see Figure 1b), with 1024 × 1024 × 1024 voxels and 12.9 μm/voxel.

2.1.2. Fracture Characterization

Prior to discussing how to represent the fracture geometries, we must first introduce the approaches in measuring fracture morphological properties regarding fracture apertures and fracture trend-surface.

Figure 2a shows a 2D slice of a cropped grayscale sub-image of 300 × 300 × 100 voxels from the shale CT-image (red box) shown in Figure 1a. A binary image in Figure 2b, consisting of only solid and fracture void, can be obtained by segmenting the grayscale image using the threshold segmentation method [46]. Thus the fracture volume fraction (or porosity) of the 3D shale sub-image can be calculated as 9%, which is the ratio between the fracture void and rock volume. Meanwhile, based on the binary image, three structural components of a fracture, as shown in Figure 2c, can be extracted by the edge detection operator [46]. Defining the plane of height Z = 0 as a reference, the spatial relative positions of the upper and lower fracture-wall surfaces at a coordinate ζ = (x, y) on the plane can be determined, denoted by z^U(ζ) and z^L(ζ), or z^U and z^L for short. The fracture local aperture b(ζ) at a point ζ, which is described as the length of the two-way arrows in Figure 2c, can be simply calculated by:

$$b=z^U-z^L$$

(1)

It is worth noting that if the fracture profile in the rock extends horizontally at the XY plane, the fracture aperture can be approximated by the height difference between the upper and lower fracture surfaces, as referred to in [8]. If the fracture is placed non-horizontally at the XY plane, the fracture aperture, i.e., the shortest distance between fracture two surfaces, can be calculated with the help of the fracture medial surface, which the calculation details can be referred to [47,48]. In our work, all fractures are assumed horizontal extending at the XY plane, and the height difference between the upper and lower fracture surfaces is choose as the fracture aperture (see Equation (1)).

The red curve in Figure 2c was obtained by smoothing the mid-surface, z^M (see Equation (2)), using the image filter operator [46], and denoted as z^T (see Equation (3)). The heights of the trend-surface quantifies the spatial fluctuation of a fracture.

$$z^M=\left(z^U+z^L\right)/2$$

(2)

$$z^T=g\left(z^M\right)$$

(3)

where g(.) is a moving average smoothing operator.

Therefore, fracture void space can be equivalently characterized by two elements: apertures and trend-surface.

With respect to fracture apertures that vary in a range, mean aperture, μ_b, and standard deviation, σ_b, can be estimated from the aperture density distribution, and then we have the ratio of the standard deviation and the mean aperture, σ_b/μ_b, to quantify fracture roughness [19]. In addition, the spatial correlation between apertures at different locations can be described by a scalar, i.e., aperture spatial correlation length, which can be described by the semi-variogram function [49]. Let φ be a stochastic variable defined on a discrete domain Ω ⊂ R and h = (x, y) be a (direction) vector, the semi-variogram γ_φ(h) between variables φ(ζ) and φ(ζ + h) is defined as:

$${\gamma }_{\phi }\left(\mathit{h}\right)=\frac{1}{2}E\left[{\left(\phi \left(\mathsf{\zeta }\right)-\phi \left(\mathsf{\zeta }+\mathit{h}\right)\right)}^2\right]$$

(4)

where E[] denotes the mathematical expectations, ζ is evaluated over all positions in Ω, and the variable φ can be the fracture apertures (Equation (1)) or surface heights (Equation (3)). Note that the closer two locations (i.e., ζ and ζ + h) are, the smaller the semi-variogram value ${\gamma }_{\phi }$ becomes.

Similarly, for the fracture trend surface (Figure 2d), z^T, the corresponding standard deviation, σ_T, and spatial correlation length, λ_T, can be estimated using the height density distribution of z^T and the surface semi-variogram, respectively. It is noticed that σ_T can also be considered as a coefficient to characterize fracture undulation, which reveals the undulation strength on fracture plane.

In summary, in this work the key properties influencing fracture morphology are μ_b, σ_b, λ_b for fracture apertures (fields) and σ_T, λ_T for fracture trend surface.

2.2. Modified Cubic Law

A single fracture is usually modelled by two smooth parallel plates separated by a constant width (i.e., aperture), $\overline{b}$, which is approximated by the arithmetic average of apertures for a rough-walled fracture [50], i.e., $\overline{b}$ ≈ μ_b. Many researchers [13,20,43] have shown that for the steady state, isothermal, laminar flow between parallel plates, the absolute permeability, K_CL, of a fracture is a function of fracture spacing (i.e., aperture), and estimated by the following equation (i.e., the cubic law):

$$K_{CL}=\frac{{\overline{b}}^2}{12}$$

(5)

However, neglecting the morphological complexity of the void/opening space of natural fractures results in the failure using the cubic law to predict the permeability accurately due to –that $\overline{b}$ in Equation (5) is not able to sufficiently represent aperture variation, thus commonly a modified factor [43], f is introduced in Equation (5) to take into account of the fracture complex morphology as:

$$K=\frac{{\overline{b}}^2}{12\left(1+f\right)}$$

(6)

To determine f we can use the following formula:

$$f=\frac{{\overline{b}}^2}{12K}-1$$

(7)

where K is assumed to be true permeability, which is actually approximated by the LBM calculation in this work.

2.3. Fracture Modelling

In this work, each model of a rough-walled fracture is generated by combining two random two-dimensional fields into a fracture aperture map b(x, y) and a trend surface z^T(x, y). Assuming that these two fields are defined on the XY plane, and both the aperture field and the trend surface are suggested to be normally distributed and are spatially correlated. In our previous work [51], an image mean-filter method has been proposed to generate such fields independently and then to assembly into a 3D fracture, as shown in Figure 3, subject to five pre-measured morphological parameters, μ_b, σ_b², λ_b, σ_T² and λ_T. Note that the actual position of the trend surface has no impact on the aperture distribution of a resultant fracture, so μ_T can be randomly chosen in modelling.

The method includes three major steps: (1) sampling data from a uniform distribution and then filtering the data with a pre-defined mean mask to obtain a trend surface field that has the required distribution and the pre-specified correlated feature, associated with two input parameters: σ_T² and λ_T. The data field is used as fracture trend surface, illustrated in Figure 3, denoted as z^T; (2) Generating two more fields independently in the exact same means as the first step for the two fracture walls with three specified aperture parameters: μ_b/2, σ_b²/2 and λ_b. The resulting two data fields, denoted as d^U and d^L, are the relative heights in relation to the trend surface in two opposite directions, as demonstrated in Figure 3; (3) Positing the two wall fields above or below the trend surface in horizontal direction into an upper surface, z^T + d^U, and a lower surface, z^T − d^L, respectively. Let b be d^U + d^L, i.e., a local aperture at a location on XY plane. Over the whole XY plane the resultant fracture is determined by five statistical parameters: mean aperture μ_b, aperture variance σ_b², aperture correlation length λ_b, trend surface variance σ_T² and trend surface correlation length λ_T.

As an example, Figure 4 shows three 3D fracture models generated using our mean-filter method with five input morphological parameters. Note that the fracture shapes are different from each other even though with the same input morphological parameters, which gives a certain amount of freedom in fracture shape but capture the controlling morphology.

2.4. Fracture Permeability Calculation

On 3D fracture models, single-phase can be simulated using Palabos [52], Lattice-Botlzmann method (LBM) software, which is an open-source code for computational fluid dynamics. As the LBM is highly flexible and parallel implemented, it can handle a wide variety of fluid dynamics problems [40,41,42,43]. In this study, the Palabos LBM software is utilized to compute absolute permeabilities only for a large number of fracture models under an assumption of the bounce-back boundary condition for fracture walls. The fracture walls, as shown in Figure 5, parallel to the flow direction, are padded to eliminate the need for periodic boundaries. The standard Bhatnagar-Gross-Krook (BGK) collision operator is applied with a D3Q19 model. Initially, the fluid macroscopic density in the inlet is set to one and the fluid velocity is set to zero everywhere and fluid movement is induced by maintaining a constant pressure gradient between the inlet and outlet, which is set as 5 × 10⁻⁵ Pa in addition to a lattice viscosity of 0.167 for all models. Steady state will be reached when the standard deviation of the average energy, measured over a fixed number of time steps, fall below a given threshold value (e.g., 3000 time steps, or a threshold value of 10⁻⁴). Steady-state results of the model are then analyzed to simulate absolute permeability. Given the velocity distribution throughout a 3D domain, the “true” permeability K in Equations (6) and (7) of a fractured rock model is calculated by the Darcy’s law:

$$-\frac{dP}{dx}=\frac{\mu }{K}U$$

(8)

where dP/dx is the pressure gradient along X-direction, μ is the fluid viscosity, and U is the average velocity through a cross-section.

Actually, the fractured rock permeability, calculated by Palabos code, is in dimensionless lattice (or image voxel) units and marked as K_l hereafter. This value is easily converted to physical units K_r by Equation (9) [45]:

$$K_r=K_l·l^2$$

(9)

where l is the physical side-length of a voxel, as shown in Figure 5.

Considering that the rock models studied in our work contain a single fracture only, the fracture permeability K_f is computed instead by [53]:

$$K_f=\frac{K_r}{\varnothing }$$

(10)

where φ is the volume fraction of fracture void in the rock models.

Since the computation of permeability using the LBM has been evaluated as a reliable approximation to the real permeability, which have been validated by comparing simulated results with experimental data in reference [45]. In our work we also denote the permeability calculated by the LBM as K_LB, i.e., K_f = K_LB.

3. Results and Discussion

3.1. Robustness of Fracture Modelling

Regarding the fracture models generated by our mean-filter method, the fracture shapes demonstrate a certain amount of randomness between individual models for the same input morphological parameters, i.e., μ_b, σ_b, λ_b, σ_T², λ_T, see Figure 4. Further, it is worthwhile investigating the robustness of these models in terms of their permeabilities before we carry out studying the effect of fracture morphology on single-phase flow. In other words, we need to assess how robust the fracture models are that have the same input parameters on fracture permeability rather than fracture shapes. For instance, μ_b = 8.7 voxels, σ_b² = 3, λ_b = 6 voxels, σ_T² = 61, and λ_T = 200 voxels are the five morphological properties measured from the fracture illustrated in Figure 2, and can be used as input parameters to generate any number of fracture models (realizations). Figure 6 shows the relative errors in the LBM permeability between such 15 realizations and a nature fracture in Figure 2, in which the relative error ε, is defined as |k − k_o|/k_o, where k_o and k are two permeabilities of the natural fracture and one of the realizations, respectively. It is evident that different fracture realizations for the same input parameters differ only within a narrow range (ε < 15%), which indicates that our modelling method captures the essential morphology that dominates single-phase flow. Although the numerical fractures are not identical entirely to the actual fracture in shape, their corresponding permeability is sufficiently close, therefore the method is robust in this sense.

3.2. Influence of Fracture Trend Surface

The morphological parameters (i.e., σ_T², λ_T) are introduced in Section 2 to characterize the fracture trend surface. In order to study the impact of the trend surface on petro-physical properties (especially the permeability) associated with fractures, three groups of single fracture models of 300 × 300 × 100 voxels are generated using the algorithm briefly introduced in Section 2.3, within each group every fracture has the same mean aperture: 5, 9, 13 (in voxels) for group 1, 2, 3, respectively. In each group, there are 25 fracture models which generated by pairing these two sets of trend surface parameters: σ_T² = 10, 50, 90, 130, 170; λ_T = 120, 160, 200, 240, 280 (in voxels), meanwhile, to eliminate the influence of fracture aperture field, we fix σ_b² = 0 and λ_b = ∞ (i.e., constant apertures at everywhere). Only six fracture models are shown in Figure 7 to demonstrate the different morphology with respect to σ_T² and λ_T. Using the LBM, 25 fracture permeabilities can be calculated for each group, followed by calculating 25 modified factors (f_T) according to Equation (7).

With the correlation lengths (λ_T) and the undulation coefficients (σ_T²) about the fracture trend surfaces in each group, we relate the corresponding modified factors, as shown in Figure 8a,d, Figure 8b,e, and Figure 8c,f for group 1, 2 and 3 of mean apertures of 5, 9 and 13 voxels, respectively. It is observed that with fixed undulation coefficient and mean apertures, the modified factor is a power-law function of the correlation length as the dotted curves shown in Figure 8a–c. As the longer the correlation length is, the gentler the bending of the trend surface is, consequently the fracture wall is smoother and therefore the modified factor is smaller. This means that the flow path in the fracture is shorter, thereby increasing the fracture permeability. However, when the undulation coefficient decreases, the fracture models are more plate-like. In this circumstance, for instance σ_T² = 10 (see Figure 8a–c), changing the correlation length of the fracture trend surface hardly affects the fluid flow path in the fracture, thus the modified factor tends to be zero.

With fixed correlation lengths and mean apertures, Figure 8d–f reveals the relationship between the modified factor and the undulation coefficient for the three groups. It can be seen that when correlation length remains constant, the modified factor increases linearly as the undulation coefficient increases. The undulation coefficient reflects fracture tortuosity, and so the greater the tortuosity is, the longer the flow path is, therefore the fracture permeability decreases. In addition, in Figure 8d, we found that there are some crosses existing. In Figure 8d, the curves with different colors represent different correlation lengths (λ_T) of trend surface, and the distance between curves decreases with the increase of correlation length (λ_T). When the correlation length is 240 (blue curve) and 280 (purple curve) in Figure 8d, the two curves are closest to each other. Therefore, the crossing of curves in Figure 8d may be caused by numerical errors.

Depicting a triple of f_T, σ_T² and (dimensionless) λ_T/L as a dot in Figure 9 for every fracture in the three groups, a 3D surface can be fitted by the least square fitting method. The fitted function (Equation (11), R² = 0.86) for the surface can be obtained as:

$$f_T=0.064·{\left(\frac{{\lambda }_T}{L}\right)}^{-1.39}+0.0012·{\sigma }_T^2-0.084$$

(11)

This reflects that the modified factor increases linearly with the undulation coefficient, and decreases with the correlation length as a power function. In addition, we found that the power in Equation (11) is −1.39 and the power of the dotted lines in Figure 8b,c are closed to −1.39, while the value in Figure 8a is −0.99 and quite larger than −1.39. The reason of this phenomenon is that the dotted lines in Figure 8a–c are the optimal fitting of the power formula for the black solid lines. Although, the power in Figure 8a differs greatly from that in Figure 8b,c, the coefficient of power function in Figure 8a is much smaller than that in Figure 8b,c, which make up for the difference between the formulas in Figure 8a–c.

3.3. Influence of Fracture Aperture Variation

Regarding to the morphological parameters (i.e., μ_b, σ_b, λ_b) for fracture apertures, similarly, three groups of fracture models of 300 × 300 × 100 voxels are also generated, having the mean aperture (μ_b) of 5, 9, 13 voxels for Group 1, 2 and 3, respectively. Differently the two trend surface parameters are fixed as σ_T² = 0 and λ_T = ∞ to eliminate the impact of trend surfaces, while only the aperture smoothness and aperture correlation length are allowed to pair between σ_b/μ_b = 0.1, 0.2, 0.3, 0.4, 0.5, and λ_b = 6, 10, 14, 18, 22. As a result we have 25 fracture models in each group and 75 models in total, of which six fracture models are shown in Figure 10.

Figure 11a–c presents the relationship between the modified factor (Equation (7)) and aperture correlation length based on the generated fracture models in Group 1, 2 and 3, respectively. It is evident that the modified factor seems to be a power function of the correlation length when σ_b/μ_b is constant. This is because the longer the correlation length the stronger the continuity between local apertures. Consequently, the obstruction of fluid flow in the fracture is reduced and increasing the degree of flow movability.

Similarly, Figure 11d–f describes the relationship between the modified factor and aperture roughness when λ_b/L is constant. It can be concluded that under the same correlation length of aperture fields, the modified factor is linearly related with the aperture roughness. This is because the roughness determines the range of variation between local apertures. The larger the roughness the greater the undulation of apertures, which hinder fluid flow and so reduces the fracture permeability.

Depicting a triple of f_b, σ_b/μ_b and λ_b/L as a dot in Figure 12 for every fracture in the three groups, a 3D surface can be fitted in the same way as above. The fitted function (R² = 0.83) for the surfaces can be obtained as:

$$f_b=1.18·{\left(100·\frac{{\lambda }_b}{L}\right)}^{-0.7}+2.73·\frac{{\sigma }_b}{{\mu }_b}-0.49$$

(12)

It shows that the modified factor is a linear function of aperture roughness, but a power function of the aperture correlation length.

3.4. Accumulated Influence of Fracture Morphology

Quantitative studies on how the variation in both fracture trend surface and fracture apertures effects the modified factor in the cubic law (Equation (6)) have been discussed above separately. The accumulated influence may be complicated and has no systematical way to study, however it is necessary to investigate the function relationship between the two factors, f_T and f_b, and the total modified factor f in Equation (7).

In the same way as described in above sections, three groups of fracture models of 300 × 300 × 100 voxels can be generated, within each group every fracture has the same mean aperture: 5, 9, or 13 (in voxels), respectively. In each group there are 50 fracture models, of which 25 fracture models were generated by altering the trend surface parameters: σ_T² = 10, 50, 90, 130, 170, λ_T = 120, 160, 200, 240, 280, where the fracture aperture parameters were fixed as σ_b/μ_b = 0.5 and λ_b = 22; meanwhile the other 25 fracture models were generated by changing the aperture parameters: σ_b/μ_b = 0.1, 0.2, 0.3, 0.4, 0.5; λ_b = 6, 10, 14, 18, 22, but the fracture trend surface parameters were fixed as σ_T² = 100 and λ_T = 120. As a result, we have 150 models in total, of which eight fracture cross-sections are shown in Figure 13 to demonstrate the different morphology with respect to (σ_T², λ_T) or (σ_b/μ_b, λ_b). Using the LBM, we have 150 fracture permeabilities and the corresponding modified factors f, furthermore we have the corresponding f_T and f_b for each fracture model according to Equation (11) and Equation (12).

Depicting a triple of f, f_T and f_b as a dot in Figure 14 for every fracture in the three groups, a 3D surface can be fitted for the data by the least square fitting method. The corresponding function (R² = 0.85) for the fitted surface can be estimated by:

$$f=0.77·f_T+0.93·f_b$$

(13)

It is shown that the slope of the modified factor caused by the aperture parameters (f_b) is larger than the slope of the trend surface parameters (f_T) because more weight is assigned, which indicates that the morphological characteristics of fractures described by (μ_b, σ_b and λ_b) are the predominant factor that effect the fluid flow more significantly. Therefore, the modified factor in Equation (13) is used to refine the cubic law as given in Equation (6). Next, we will use natural fractures in CT images to validate the calculation of the modified factor by Equation (13).

4. Verification with Natural Fractures

Figure 15 shows seven slice images of a CT image acquired by scanning seven fractured rock cores, all the images have the same dimensions of 1024³ voxels, and the voxel resolution is assumed to be 1 μm/voxel without loss of generality. The first five are shales and the last two are sandstones. Instead of using the whole original images, a small sub-image (see Figure 15 is cropped out from each CT image to reveal a natural fracture. All fractures in Figure 15 can be considered to be horizontally placed. Using the measurement method described in Section 2.1.2, we can measure all the geometrical parameters of the trend surface and fracture aperture fields in the CT images, as listed in

Table 1. The characteristic parameters of eight different fracture samples presented in the green box of Figure 15 and the sample permeabilities calculated from LBM, cubic law and modified cubic law.

Samples	μb (μm)	σb/μb	λb (μm)	σT2	λT (μm)	KLB (mD)	Kcl (mD)	Kmcl (mD)	εcl (%)	εmcl (%)
SA1A	5.03	0.26	6	160.75	150	1110.49	2108.41	982.76	89.86	−11.5
SA1B	5.28	0.19	6	90.60	150	1443.05	2323.20	1232.62	60.99	−14.58
SA1C	5.52	0.25	10	231.2	200	1157.60	2539.20	1297.74	119.35	12.11
SA2	10.45	0.20	20	14.40	150	6166.89	9100.21	6372.08	47.57	3.33
SA3	11.75	0.16	10	161.14	200	6705.99	11,505.21	7434.62	71.57	10.87
SA4	6.74	0.18	10	45.99	150	2665.43	3785.63	2389.23	42.03	−10.36
SA5	11.45	0.14	10	18.6	150	8102.42	10,925.21	8551.35	34.84	5.54
SB1	5.25	0.18	7	94.25	150	1468.85	2296.88	1271.76	56.37	−13.4
SB2	10.81	0.30	20	70.5	140	5689.30	9738.01	5599.80	71.16	−1.57

. Three types of fracture permeabilities for each fracture sample were calculated using LBM, cubic law, and our modified cubic law, K_LB, K_cl, K_mcl, respectively. Let K_LB be the permeability calculated by the LBM, and K_cl and Kmcl be the permeabilities estimated by the original cubic law of Equation (6) and our modified cubic law (Equations (7) and (13)), respectively. The relative errors are ε_cl = (K_cl − K_LB)/K_LB and ε_mcl = (K_mcl − K_LB)/K_LB. Note that from the table that the relative errors calculated by the modified cubic law are far less than that of the original cubic law: Most of them, remain below 20%, and the average relative error is ~8.57%, which are acceptable in practical applications. Therefore, based on the combination of fracture trend surface and variable apertures, the modified cubic law ((Equations (7) and (13)) proposed in this work effectively calculates the permeability of rough fractures, greatly reducing the computational complexity.

5. Conclusions

Advances in 3D imaging technology allow the fracture structure in fractured rocks to be revealed truly and faithfully without destroying the core. In this study, the rough-walled fracture was investigated by combining fracture variable apertures and the fracture trend surface topography, which quantifies undulation about the fracture undulation. Using the computer laboratory, the geometric characteristics of variable aperture and fracture trend surface were quantitatively represented, in which the aperture roughness and aperture correlation length used to describe variable apertures, and the trend surface fluctuation coefficient and the trend surface correlation length used to represent the fracture trend surface. In addition, the influences of these representation parameters on the fracture permeability of fluid flow through rough-walled fracture were investigated. Lastly, a corrector factor is proposed for improving the performance of the cubic law by incorporating the characteristic parameters of fracture trend surface and variable apertures. For all the 3D digital fracture samples present in this work, the modified cubic law was determined to be an acceptable approximation to the numerical simulation (LBM).

Because single fracture is the basic element of fracture network, the methodology has the potential application in characterizing the geometric structure of fracture network in which the interaction between fractures can be take into account. Meanwhile, the flow properties of the fracture network can be evaluated by the modified cubic law proposed in this work. These tasks will be processed in our next work. Apart from this extension, the applicability of our methods proposed in this work may be limited within a narrow range of length-scale (fracture opening in μm, mm or cm, plus the model dimensions), we will explore both in the near future, however the inertial effect and surface contact are neglected in our work and will be investigated in the future.