1 Introduction

The fractures are usually of multi-scale heterogeneity and thus affect the nonlinear flow of fluids in geological formations. Fractures are the main fluid flow channels in geological formations for nuclear waste disposal (Wu and Wang 2020), geothermal resource development (Xin et al. 2020; Ye and Wang 2020), oil and gas exploitation (Zheng et al. 2020; Sharifi et al. 2021), and CO2 geological storage (Wang and Wang 2018). These fractures may naturally form fracture networks by the expansion and connection of microcracks, holes and joint fissures. Artificial stimulations such as hydraulic fracturing can help to form a fracture network in shale gas reservoirs (Dou et al. 2021). A fracture network has two features: rough surface of each fracture and fractal tree-like branching, as shown in Fig. 1. The surface morphology of natural fractures is always rough, with self-affine surfaces and non-uniform aperture. Natural micro-fractures are interconnected and can be simplified as a fractal tree-like fracture network. These two features heavily affect the gas transport ability or permeability of this fracture network through the fracture aperture, connectivity, fracture length and surface roughness (Wang and Cheng 2020; Zheng et al. 2022; Wan et al. 2023). However, an analytical permeability model for the gas flow in the fracture network with rough surfaces has not been available.

Fig. 1
figure 1

The rough surface and fractal tree-like branching characteristic of nature fractures

Full size image

The gas flow rate is usually calculated by the classical cubic law for fractures with two smooth parallel surfaces (Snow 1969; Witherspoon et al. 1980; Akbarzadeh and Chalaturnyk 2014; Feng et al. 2020). Numerous laboratory tests observe that the classical cubic law produces obvious errors in the calculation of the flow rate in rough fractures (Ju et al. 2013; Tzelepis et al. 2015; Liu et al. 2016; Yin et al. 2017). Thus, the cubic law for smooth fractures is not applicable to the rough fracture structures.

The cubic law has been modified for the gas flow in rough fractures. Jin et al. (2017) developed a permeability model for single fracture whose rough fracture profile is generated by Weierstrass–Mandelbrot (W–M) function. They analyzed the impacts of surface roughness on gas transport by considering the surface tortuosity, surface roughness and hydraulic tortuosity of a rough fracture. Zeng et al. (2018) predicted the permeability of the microchannels with self-affine rough surfaces through shale gas flow simulations. They found that the influence of surface roughness on shale gas transport is greater in the transition flow region than in the slip flow region. Wang et al. (2018) established a modified cubic law to estimate the flow rate of single-phase fluid in a rough fracture. Their model has a simple form compared to the standard cubic law and can more accurately describe the fluid flow properties in rough fractures. Zhang et al. (2019) simulated the water flow in the rough fracture with several rough elements of rectangular or triangular bulges regularly arranged on the fracture surface. They found that these rough elements increase flow resistance and cause the separation and backflow of flow streamlines. Ju et al. (2019) proposed a water permeability model for rough fractures based on an optimized segmentation algorithm and self-affine rough profiles. Their self-affine rough profiles are generated by the W–M function and the optimized segmentation algorithm can decompose the fracture profiles into smaller parts. Dong and Ju (2020) further extended their optimized segmentation algorithm to the fractures with variable aperture and modified the cubic law. They used the lattice Boltzmann method (LBM) to simulate the water flow through the fracture with variable aperture (Guo et al. 2002). Liu et al. (2020a) proposed a successive random addition algorithm to generate a 3D rough surface of fractures and simulated the fluid flow through the 3D rough fractures. These previous studies only investigated the fluid flow in single rough fracture. The permeability and gas flow mechanisms within a rough fractal tree-like fracture network have not been investigated so far.

Fracture network has complex morphology and plays a key role in fluid flow mechanisms in shale oil or gas reservoir (Jafari and Babadagli 2012; Miao et al. 2015; Xu et al. 2016; Hu et al. 2019; Lahiri 2021). Jin et al. (2019) constructed several cleat networks in coal reservoirs based on self-affine rough fracture profiles and carried out the LBM simulation of gas flow in these cleat networks. They established an empirical formula to predict the permeability of cleat networks in coal reservoirs. Hu et al. (2020) generated several fractal discrete fracture networks in a shale gas reservoir based on the fractal distribution of fracture length. They simulated the shale gas transport in different fractal discrete fracture networks and observed that the cumulative production of shale gas is increased by 113% when fractal dimension increases from 1.5 to 1.7. Previous studies demonstrated that the distributions of fracture apertures and lengths follow fractal scaling laws, and corresponding permeability models for the fracture network have been developed based on fractal theory (Miao et al. 2015; Li et al. 2016; Wang et al. 2019; Gao et al. 2021). However, these previous permeability models assume that the fracture network is composed of a series of parallel fractures. This assumption simplifies the structure of fracture network but ignores the connectivity between fractures. This is derivated from the fact that the natural fracture networks in a shale gas reservoir are usually composed of intersected fractures. Therefore, the accurate description of the distribution of the intersected fractures is a challenge.

A complex fracture network can be described by the fractal tree-like fracture network model composed of a main fracture and several secondary branching fractures (Lorente et al. 2002; Xu and Yu 2006; Zheng and Yu 2012; Alalaimi et al. 2015; Zhang et al. 2021). When high pressure fluid is injected into a shale gas reservoir, hydraulic fractures interact with primary fractures to create more small fractures, thus increasing effective gas flow paths and enhancing shale gas production efficiency (Zhang et al. 2018; Dou et al. 2019). Xu et al. (2016) believed that the fractal tree-like fracture network can capture the statistical characteristics and topology of the fracture system in fractured porous media. They derived analytical expressions of effective permeability, tortuosity and fracture density for fractured porous media. Peng et al. (2019) proposed a thermal conductivity and an effective permeability model for a fractal tree-like fracture network and investigated the effects of structure parameters, such as length ratio, branch angle, branch level and width ratio on the heat transfer and fluid flow in the fracture network. Liu et al. (2020b) and Huang et al. (2020) developed their theoretical permeability models for the fractal tree-like fracture network with constant and variable aperture width. They found that the permeability of a fracture network with variable aperture is smaller than that of fracture network with constant aperture. These fractal tree-like fracture network models have been successfully applied to the complex fracture system of fractured porous media, but they are based on the smooth fracture surfaces in each branching and ignore the effect of surface roughness on the gas transport in complex fracture networks. The rough surfaces in each branching may have significant influences on the shale gas transport in a complex fracture network and should be further investigated.

In this study, an analytical permeability model is proposed for the rough fractal tree-like fracture network in fractured shale gas reservoirs. It considers the rough surface and branching characteristic of fracture network, which can effectively predict the permeability of fracture network of shale gas reservoirs, and is of great significance for unconventional oil and gas extraction. This paper is organized as follows. In Sect. 2, a rough fractal tree-like fracture network is constructed by self-affine rough surface profiles and branching characteristics, and an analytical permeability model is proposed to consider both surface roughness and branching characteristic of a rough fractal tree-like fracture network. In Sect. 3, the shale gas flow in the rough fractal tree-like fracture network is numerically simulated by the computational fluid dynamics (CFD) module within COMSOL software. In Sect. 4, the analytical permeability model is verified with the CFD simulations. In Sect. 5, the gas velocity distribution in the rough fractal tree-like fracture network at both micron and millimeter scales and the effects of structure parameters such as fractal dimension, length ratio, aperture ratio, branch level and branch angle on permeability are analyzed. Main conclusions are drawn in Sect. 6.

2 An analytical permeability model for a rough fractal tree-like fracture network

2.1 Construction of rough fractal tree-like fracture network

The fracture network is of great significance for the gas transport and permeability of shale gas reservoirs. Following Weierstrass–Mandelbrot (W–M) random function is used to characterize the rough fracture surface with self-affine property (Berry and Lewis 1980; Jin et al. 2017; Ju et al. 2017):

$$Z(x) = G^{{D_{f} - 1}} \sum\limits_{{i = N_{l} }}^{{N_{h} }} {\gamma^{{i(D_{f} - 2)}} } \left[ {\cos (\phi_{i} ) - \cos \left( {\frac{{2\pi \gamma^{i} x}}{{L_{s} }} - \phi_{i} } \right)} \right]$$
(1)

where \(Z(x)\) is the surface height at the coordinate x; \(D_{f}\) is the fractal dimension of fracture surfaces; \(\gamma\) is the rescaling factor that controls the density of the spectrum; \(G\) is the scaling constant; \(\phi_{i}\) is the random number that ranges from 0 to 2π; \(L_{s}\) is the sample length; \(N_{h}\) and \(N_{l}\) control the maximum cosine frequency and minimum cosine frequency, respectively.

The rough fracture surface profiles with different fractal dimensions \(D_{f}\) are shown in Fig. 2. These profiles show that the fractal dimension \(D_{f}\) has a great influence on the morphology of fracture surfaces. A larger fractal dimension has a larger height range of surface profile and more obvious frequency variation. Figure 3 shows the diagrams of single rough fracture surface profile, a single rough fracture and a rough fractal tree-like fracture network. A single rough surface profile moves vertically for a short distance \(e_{0}\), where \(e_{0}\) can be regarded as the fracture aperture. A single rough fracture is formed by two rough surface profiles with the same shapes as shown in Fig. 3b. This single rough fracture is divided into two identical sub-branches according to the constant branching angle \(\theta\). A binary rough fracture architecture is shown in Fig. 3c. The length ratio \(\alpha\) and aperture ratio \(\beta\) of the two adjacent branching levels are defined as

$$\alpha = \frac{{l_{k + 1} }}{{l_{k} }}$$
(2)
$$\beta = \frac{{e_{k + 1} }}{{e_{k} }}$$
(3)

where \(l\) and \(e\) are the fracture length and aperture; the subscript \(k\) or \(\left( {k + 1} \right)\) denotes the fracture branching level. Thus, the fracture length and aperture at the kth level are

$$l_{k} = l_{0} \alpha^{k}$$
(4)
$$e_{k} = e_{0} \beta^{k}$$
(5)
Fig. 2
figure 2

Rough fracture surface profiles generated by the W-M function with different fractal dimensions

Full size image
Fig. 3
figure 3

Construction of a rough fractal tree-like fracture network

Full size image

2.2 A permeability model for single rough fracture

The gas flow rate \(q\) in a single fracture with two parallel flat plates is calculated by the classical cubic law:

$$q = \frac{{we^{3} }}{12\mu }\frac{\Delta P}{l}$$
(6)

where \(w\) is the fracture width; \(\Delta P\) is the pressure difference between the inlet and the outlet; \(\mu\) is the fluid viscosity.

The rough surface morphology in Fig. 2 has a great influence on shale gas flow. The roughness has three effects (Jin et al. 2017): the local roughness effect \(f_{local}\), the surface tortuosity effect \(f_{s}\), and the hydraulic tortuosity effect \(f_{\tau }\). Their total roughness effect \(f_{total}\) is then expressed by

$$f_{total} = f_{local} f_{s} f_{\tau }$$
(7)

The local roughness effect \(f_{local}\) is semi-empirically expressed as (Witherspoon et al. 1980)

$$f_{local} = 1 + \alpha_{c} \left( {\frac{\sigma }{e}} \right)^{{\alpha_{e} }}$$
(8)

where \(\sigma\) is the root mean square of surface height; \(\alpha_{c}\) and \(\alpha_{e}\) are the two fitting parameters related to the surface geometry.

The tortuous surface geometry narrows fracture aperture which is modified as \(e_{m} = {e \mathord{\left/ {\vphantom {e {\tau_{s} }}} \right. \kern-0pt} {\tau_{s} }}\). Thus, the surface tortuosity effect \(f_{s}\) is (Turcotte 1997)

$$f_{s} = \frac{{e_{m}^{2} }}{{e^{2} }} = \frac{1}{{\tau_{s}^{2} }}$$
(9)

where \(\tau_{s} = {{A_{s} } \mathord{\left/ {\vphantom {{A_{s} } {A_{0} }}} \right. \kern-0pt} {A_{0} }}\) is the surface tortuosity; \(A_{s}\) and \(A_{0}\) are the total area of the rough surface and its corresponding projection area. The rough surface area is calculated through the fractal scaling law as (Sun and Koch 1998)

$$A_{s} = A_{0}^{{{{D_{s} } \mathord{\left/ {\vphantom {{D_{s} } 2}} \right. \kern-0pt} 2}}} e^{{2 - D_{s} }}$$
(10)

where \(D_{s}\) is the fractal dimension of surface tortuosity.

Thus, the surface tortuosity \(\tau_{s}\) is

$$\tau_{s} = \frac{{A_{0}^{{{{D_{s} } \mathord{\left/ {\vphantom {{D_{s} } 2}} \right. \kern-0pt} 2}}} e^{{2 - D_{s} }} }}{{A_{0} }}{ = }\frac{{e^{{2 - D_{s} }} }}{{A_{0}^{{1 - {{D_{s} } \mathord{\left/ {\vphantom {{D_{s} } 2}} \right. \kern-0pt} 2}}} }}{ = }\left( {\frac{e}{{\sqrt {A_{0} } }}} \right)^{{2 - D_{s} }}$$
(11)

The gas flow rate \(q\) is replaced with \({q \mathord{\left/ {\vphantom {q \tau }} \right. \kern-0pt} \tau }\) to consider the hydraulic tortuous effect. The hydraulic tortuosity \(\tau\) of a rough fracture is (Duda et al. 2011)

$$\tau = \frac{U}{{U_{x} }}$$
(12)

where \(U\) is the average flow velocity in fractures and \(U_{x}\) is the component of the average flow velocity \(U\) in the x direction. The ratio of \({U \mathord{\left/ {\vphantom {U {U_{x} }}} \right. \kern-0pt} {U_{x} }} = {{({q \mathord{\left/ {\vphantom {q \tau }} \right. \kern-0pt} \tau })} \mathord{\left/ {\vphantom {{({q \mathord{\left/ {\vphantom {q \tau }} \right. \kern-0pt} \tau })} {q_{x} }}} \right. \kern-0pt} {q_{x} }}\) is satisfied in the single fracture. The hydraulic tortuosity effect \(f_{\tau }\) is expressed as

$$f_{\tau } = \frac{{q_{x} }}{q} = \frac{1}{{\tau^{2} }}$$
(13)

The tortuous fracture length \(l_{t} (e)\) is related to the fracture aperture \(e\) and the tortuosity fractal dimension \(D_{\tau }\) as (Yu and Cheng 2002)

$$l_{t} (e) = l^{{D_{\tau } }} e^{{1 - D_{\tau } }}$$
(14)

The hydraulic tortuosity \(\tau\) is the ratio of tortuous length \(l_{t} (e)\) to the straight length \(l\) (Cai and Yu 2011):

$$\tau { = }\frac{{l_{t} \left( e \right)}}{l} = \left( \frac{e}{l} \right)^{{1 - D_{\tau } }}$$
(15)

For a rough fracture, the \(D_{\tau }\) and \(D_{s}\) are approximately equal to \(D_{f}\) and \(D_{f} + 1\), respectively (Mandelbrot 1983). In order to eliminate the limitation that the fractures are composed of the parallel and smooth plates as described in Eq. (6), the total roughness effect should include the local roughness effect, the surface tortuosity effect and hydraulic tortuosity effect. Thus, the gas flow rate in a single rough fracture is

$$q = \frac{{we^{3} }}{{12\mu f_{total} }}\frac{\Delta P}{l} = \frac{{we^{3} }}{{12\mu f_{local} }}\frac{1}{{\tau^{2} \tau_{s}^{2} }}\frac{\Delta P}{l}{ = }\frac{{we^{3} }}{{12\mu f_{local} }}\left( \frac{e}{l} \right)^{{4(D_{f} - 1)}} \frac{\Delta P}{l}$$
(16)

2.3 A permeability model for a rough fractal tree-like fracture network

Each fracture branch at any branching level of a fractal tree-like fracture network can be regarded as a single fracture. Based on Eq. (16), the gas flow rate \(q_{k}\) of each fracture branch at the kth level is

$$q_{k} = \frac{{we_{k}^{3} }}{{12\mu f_{local} }}\left( {\frac{{e_{k} }}{{l_{k} }}} \right)^{{4(D_{f} - 1)}} \frac{{\Delta P_{k} }}{{l_{k} }}$$
(17)

Then, the total flow rate \(Q\) at the kth level in all fracture branches is

$$Q = n^{k} q_{k} = \frac{{we_{0}^{{4D_{f} - 1}} }}{{12\mu f_{local} l_{0}^{{4D_{f} - 3}} }}\left( {\frac{{n\beta^{{4D_{f} - 1}} }}{{\alpha^{{4D_{f} - 3}} }}} \right)^{k} \Delta P_{k}$$
(18)

where \(\Delta P_{k}\) is the pressure difference of the kth fracture branch.

The total pressure difference \(\Delta P\) of a fractal tree-like network is the sum of all \(\Delta P_{k}\):

$$\begin{aligned} \Delta P = & \sum\limits_{k = 0}^{m} {\Delta P_{k} } = \sum\limits_{k = 0}^{m} {\frac{{12\mu f_{local} l_{0}^{{4D_{f} - 3}} Q}}{{we_{0}^{{4D_{f} - 1}} }}\left( {\frac{{\alpha^{{4D_{f} - 3}} }}{{n\beta^{{4D_{f} - 1}} }}} \right)^{k} } \\ { = } & \frac{{12\mu f_{local} l_{0}^{{4D_{f} - 3}} }}{{we_{0}^{{4D_{f} - 1}} }}\frac{{1 - \left( {\frac{{\alpha^{{4D_{f} - 3}} }}{{n\beta^{{4D_{f} - 1}} }}} \right)^{m + 1} }}{{1 - \frac{{\alpha^{{4D_{f} - 3}} }}{{n\beta^{{4D_{f} - 1}} }}}}Q \\ \end{aligned}$$
(19)

where m is the branch level.

This total gas flow rate Q is expressed by the Darcy’s law as

$$Q{ = }A\frac{{K_{F} }}{\mu }\frac{\Delta P}{L}{ = }wE\frac{{K_{F} }}{\mu }\frac{\Delta P}{L}$$
(20)

where \(K_{F}\) is the permeability; \(A = wE\) is the cross-sectional area of rough fractal tree-like network. The height \(E\) of the rough fractal tree-like network is

$$\begin{aligned} E = & e_{0} + 2\left( {l_{1} \sin \theta + l_{2} \sin \theta + \cdots + l_{m} \sin \theta } \right) \\ { = } & \, e_{0} + 2l_{0} \alpha \sin \theta \frac{{1 - \alpha^{m} }}{1 - \alpha } \\ \approx & 2l_{0} \alpha \sin \theta \frac{{1 - \alpha^{m} }}{1 - \alpha } \\ \end{aligned}$$
(21)

The length \(L\) of the rough fractal tree-like network is

$$\begin{aligned} L = & l_{0} + l_{1} \cos \theta + l_{2} \cos \theta + \cdots + l_{m} \cos \theta \\ { = } & l_{0} \left( {1 + \alpha \cos \theta \frac{{1 - \alpha^{m} }}{1 - \alpha }} \right) \\ \end{aligned}$$
(22)

Substituting Eqs. (19), (21) and (22) into Eq. (20) obtains the permeability of the rough fractal tree-like network as

$$K_{F} = \frac{{e_{0}^{{4D_{f} - 1}} }}{{12f_{local} l_{0}^{{4D_{f} - 3}} }}\frac{{\left( {1 + \alpha \cos \theta \frac{{1 - \alpha^{m} }}{1 - \alpha }} \right)}}{{2\alpha \sin \theta \frac{{1 - \alpha^{m} }}{1 - \alpha }}}\frac{{1 - \frac{{\alpha^{{4D_{f} - 3}} }}{{n\beta^{{4D_{f} - 1}} }}}}{{1 - \left( {\frac{{\alpha^{{4D_{f} - 3}} }}{{n\beta^{{4D_{f} - 1}} }}} \right)^{m + 1} }}$$
(23)

Equation (23) is our analytical permeability model for a rough fractal tree-like fracture network. When \(f_{local} = 1\) and \(D_{f} = 1\), which ignores the surface roughness, the permeability becomes

$$K_{F} = \frac{{e_{0}^{3} }}{{12l_{0} }}\frac{{\left( {1 + \alpha \cos \theta \frac{{1 - \alpha^{m} }}{1 - \alpha }} \right)}}{{2\alpha \sin \theta \frac{{1 - \alpha^{m} }}{1 - \alpha }}}\frac{{1 - \frac{\alpha }{{n\beta^{3} }}}}{{1 - \left( {\frac{\alpha }{{n\beta^{3} }}} \right)^{m + 1} }}$$
(24)

When \(m = 0\), which ignores the branching characteristic, the permeability model of Eq. (23) is simplified into

$$K_{F} = \frac{{e_{0}^{{4D_{f} - 1}} }}{{12f_{local} l_{0}^{{4D_{f} - 3}} }}$$
(25)

This is the permeability of single rough fracture and has the same expression as that proposed by Jin et al. (2017).

3 Numerical simulation on gas flow in a rough fractal tree-like fracture network

The gas flow in fractures is usually described by the Navier–Stokes equation, which consists of the conservation laws of mass and momentum as (Navier 1827; Stokes 1845):

$$\frac{\partial \rho }{{\partial t}} + \nabla {\text{u}} = 0$$
(26)
$$\rho \frac{{\partial {\text{u}}}}{\partial t} + \rho \left( {{\text{u}} \cdot \nabla } \right){\text{u}} = - \nabla p + \mu \nabla \cdot \left( {\nabla {\text{u}} + \left( {\nabla {\text{u}}} \right)^{T} } \right)$$
(27)

where \({\text{u}}\) is the gas velocity; \(\rho\) is the gas density.

This Navier–Stokes equation is numerically solved by the computational fluid dynamics (CFD) module in COMSOL software. The geometrical structure of a rough fractal tree-like fracture network is shown in Fig. 4. For shale gas flow, the left boundary is the inlet and the right boundaries are the outlets. The pressure gradient is set between 1.4 and 12 kPa/m (Yin et al. 2020). The gas pressure of the inlet boundary is 120 Pa and the gas pressure of the outlet boundary is 0.01 Pa. The fracture surfaces are highly irregular, and the net movement of gas on their surfaces is almost negligible. Thus, the gas flow at the fracture surfaces can be regarded as no-slip boundary (Jin et al. 2017; Ju et al. 2019).

Fig. 4
figure 4

Geometry of a rough fractal tree-like fracture network and gas pressure distribution

Full size image

4 Model validation

4.1 Validation of simulation model

The experimental data of the single rough fracture from Yin et al. (2020) are used here to verify our simulation model. In their experiments, two single rough fractures with aperture of 1.0 mm and 1.5 mm are produced by the 3D printing technique. The fractures have the length of 100 mm and the fractal dimension of 1.2. In our simulation, the gas flow rate Q is obtained by integrating the flow velocity at the outlets. Figure 5 is the comparison of gas flow rate Q in the rough fracture between experimental data and simulation results. A clear linear relationship between gas flow rate and pressure gradient is observed when the fracture aperture is 1.0 mm. However, when the fracture aperture increases to 1.5 mm, they exhibit nonlinear characteristics, especially under high pressure gradient conditions. This may be because the gas flow rate increases with the increase of fracture aperture, resulting in a transition from linear to nonlinear gas flow. These simulation results match well with the experimental data of the two fractures. Hence, this simulation model can accurately describe the gas flow behavior in these rough fractures.

Fig. 5
figure 5

Comparison of flow rate between experimental data and simulations

Full size image

4.2 Validation of the analytical permeability model

The shale gas flow simulations in the rough fractal tree-like fracture network with e0 = 1 − 4 mm are implemented to verify the analytical permeability model of Eq. (23) for the rough fractal tree-like fracture network. Other parameters are set as: l0 = 30 mm, Df = 1.5, \(f_{local}\) = 1.05, θ = 30°, α = 0.8, β = 0.5, m = 2, and n = 2. Figure 4 also shows the gas pressure distribution in the rough fractal tree-like fracture network and shale gas flows from the inlet to several outlets at the end of the last fracture branches. Further, the permeability can be calculated by the gas flow rate Q according to Eq. (20). The permeability is compared in Fig. 6 between numerical simulation and the analytical permeability model of Eq. (23). They are generally in agreement, although some deviations are observed. The derivations are small and within an acceptable range of error. Thus, this analytical permeability model can be used to predict the permeability of the rough fractal tree-like fracture network.

Fig. 6
figure 6

Comparison of permeability between simulations and analytical model for the rough fractal tree-like network

Full size image

5 Results and discussions

5.1 Velocity distribution

The fracture apertures are set at millimeter scale in some experiments and simulations of fluid flow within rough fractures (Ju et al. 2013, 2019; Yin et al. 2017), but SEM images observe that the apertures of most fractures in shale gas reservoirs are at micron scale (Yu et al. 2022). However, the gas flow simulation in the fractures at micron scale has not well investigated so far. This study carried out a series of gas flow simulations in the rough fractal tree-like fracture networks at millimeter scale and micron scale. The aperture e0 of primary branch in the rough fractal tree-like fracture networks at millimeter scale is set as 1 mm, 2 mm, 3 mm and 4 mm, respectively, and length l0 of primary branch is 30 mm. Correspondingly, the aperture e0 of primary branch in the rough fractal tree-like fracture networks at micron scale is set as 1 μm, 2 μm, 3 μm and 4 μm, respectively, and length l0 of primary branch is 30 μm. The fractal dimension of rough surface is set as 1.2, 1.3, 1.4 and 1.5, respectively.

The velocity distribution and streamlines in the fractal tree-like fracture network with smooth surface and rough surfaces are shown in Fig. 7 at millimeter scale and in Fig. 8 at micron scale. For the millimeter scale, Fig. 7 shows that the eddy flow occurs in the local rough areas where the channel height changes greatly. These local eddy flows increase the internal frictional resistance and produce a local energy loss. However, for the micron scale, no eddy flow occurs near the rough surfaces as shown in Fig. 8. This is because the gas flow velocity is very low in the fracture networks at micron scale and cannot generate the eddy flow. The rough surface hinders the shale gas transport in fractures and results in drastic changes in gas flow velocity in the narrow area of the fractures. The gas flow velocity is much larger at millimeter scale than at micron scale.

Fig. 7
figure 7

Streamlines at different fractal dimension in rough fractal tree-like network at millimeter scale

Full size image
Fig. 8
figure 8

Streamlines at different fractal dimension in rough fractal tree-like network at micron scale

Full size image

Figure 9a is the comparison of the average gas flow velocity of the outlets at millimeter scale when the fracture surface is smooth or rough. Generally, the gas flow velocity at outlets of fracture network with smooth surface is the highest and decreases with the increase of fractal dimension. If the smooth surface case is used as the base, the growth rate of gas flow velocity with different fractal dimensions is shown in Fig. 9b for the fracture network at millimeter scale. This gas flow velocity decreases by about 50% when Df = 1.5. This phenomenon is mainly induced by two causes: One is that the rough surface increases the tortuosity of fractures and the streamline lengths of gas flow in fractures. The other is that the eddy flow more easily occurs near fracture surfaces for a higher fractal dimension. The increase of streamline lengths and eddy flow can increase the gas flow resistance and energy loss of shale gas, resulting in the decrease of gas flow velocity at the outlets. Similarly, the gas flow velocity and its growth rate at micron scale are presented in Fig. 10 for different fractal dimensions. Being different from millimeter scale, no eddy flow is observed near the rough surfaces at micron scale. The decrease of gas flow velocity at micron scale is only caused by the increase of gas flow path in fractures. As shown in Fig. 10b, the gas flow velocity at micron scale decreases by about 60% at Df = 1.5, which is larger than that at millimeter scale. The above analyses indicate that the reduction of gas flow velocity is more obvious for smaller fracture aperture and rougher fracture surface. This is because a smaller aperture will reduce the gas flow space, and a rougher surface will more easily hinder the gas flow in fractures. Thus, the gas flow resistance will be more significant for a fracture network with smaller aperture and rougher surface.

Fig. 9
figure 9

Variation of outlet velocity and its growth rate with fractal dimension at millimeter scale

Full size image
Fig. 10
figure 10

Variation of outlet velocity and its growth rate with fractal dimension at micron scale

Full size image

5.2 Impact of fractal dimension D f on permeability

Fractal dimension Df denotes the fracture surface roughness of a rough fractal tree-like fracture network. Figure 11 shows the effects of fractal dimension Df on the permeability KF of rough fractal tree-like fracture network under different e0. It can be seen that the permeability KF decreases with the increase of fractal dimension Df. The increase of fractal dimension Df means that the fractal tree-like fracture network has rougher surfaces. The local rough surface of fractures interferes with the gas flow and reduces the effective space for gas flow. Besides, higher fractal dimension of fractures also means higher tortuosity of fractures, which increases the length of shale gas flow path in fractures. The increase of gas flow length and the decrease of effective space increase the gas flow resistance, thus leading to the decrease of permeability of rough fractal tree-like fracture network.

Fig. 11
figure 11

Effects of fractal dimension Df on the permeability KF of the rough fractal tree-like fracture network

Full size image

5.3 Impacts of α and β on permeability

Length ratio α and aperture ratio β are two main scaling factors to structurally characterize the adjacent branches in a rough fractal tree-like fracture network. Figure 12 shows the effects of length ratio α on the permeability KF with different fractal dimension Df. It is found that the permeability KF logarithmically decreases with the increase of the length ratio α. Taking Df = 1.5 as an example, the permeability of the rough fractal tree-like fracture network decreases from 3.83 × 10–17 to 1.96 × 10–20 m2 when the length ratio α increases from 0.1 to 1.0. The reason is that the increase of length ratio α indicates an increase in the fracture length at the next branching levels and the total length L of rough fractal tree-like fracture network also increases. According to Eq. (22), the relationship between the total length L and length ratio α is shown in Fig. 13. The increase of length ratio α and the fracture length of the 0th level branch l0 can increase the total length L of the rough fractal tree-like fracture network. Longer total length L means bigger gas flow resistance, thus the permeability of the rough fractal tree-like fracture network decreases.

Fig. 12
figure 12

Effects of length ratio α on the permeability KF of the rough fractal tree-like fracture network

Full size image
Fig. 13
figure 13

Effects of the length ratio α on the total length L of the rough fractal tree-like fracture network

Full size image

The significant effect of aperture ratio β on permeability KF is shown in Fig. 14. The permeability KF increases with the increase of aperture ratio β under different fractal dimensions. When the aperture ratio β increases from 0.5 to 0.95, the permeability of the rough fractal tree-like fracture network with Df = 1.1 increases from 2.24 × 10–15 m2 to 2.73 × 10–14 m2, and the permeability increases from 1.46 × 10–17 m2 to 8.43 × 10–16 m2 when Df = 1.5. The permeability is much larger at Df = 1.1 than at Df = 1.5. This means that larger aperture ratio β and smoother fracture surfaces can provide larger gas flow spaces and lower gas flow resistance in sub-branching levels of the rough fractal tree-like fracture network. The effects of length ratio α and aperture ratio β on the permeability are consistent with the results from Liu et al. (2020b), but they ignore the rough surface of fracture branches in the fractal tree-like fracture network.

Fig. 14
figure 14

Effects of the aperture ratio β on permeability KF of the rough fractal tree-like fracture network

Full size image

5.4 Impacts of m and θ on permeability

The distribution of branches in the rough fractal tree-like fracture network is determined by two structural parameters: branch level m and branch angle θ. Figure 15 shows the effects of the branch level m on permeability KF at different fractal dimensions Df. The logarithm of permeability KF has an obvious negative linear relationship with the branch level m. When Df = 1.5, the permeability KF of the rough fractal tree-like fracture network decreases more significantly with the increase of branch level m than when Df = 1.1. For Df = 1.5, when the branch level m increases from 1 to 5, the permeability KF decreases from 1.64 × 10–16 to 2.10 × 10–20 m2 by more than 4 orders of magnitude. For Df = 1.1, the permeability KF decreases from 1.23 × 10–14 to 3.04 × 10–17 m2, the permeability KF decreases by approximately 3 orders of magnitude. It is because the increase of branch level m leads to the increase of the total number of smaller sub-branches. These smaller sub-branches reduce the permeability KF. Figure 16 shows the impacts of branch angle θ on the permeability KF with different fractal dimensions Df. The permeability KF decreases with the increase of branch angle θ. At Df = 1.1, the permeability KF decreases from 1.39 × 10–14 to 5.64 × 10–16 m2 by approximately two orders of magnitude when the branch angle θ increases from 5° to 85°. This is because the energy loss at the branch node increases with the increase of branch angle θ. The variations of permeability with branch angle θ and branch level m are consistent with previous studies (Wang and Cheng 2020; Hu et al. 2021), but the surface roughness of each sub-branch is not considered in their models.

Fig. 15
figure 15

Effects of the branch level m on permeability KF of the rough fractal tree-like fracture network

Full size image
Fig. 16
figure 16

Effects of the breach angle θ on permeability KF of the rough fractal tree-like fracture network

Full size image

6 Conclusions

In this study, an analytical permeability model was proposed for a rough fractal tree-like fracture network in shale gas reservoirs. This fracture network was constructed by rough surface profiles and branching characteristics. Then, the analytical permeability model was verified with CFD simulations. Finally, the gas velocity distribution and the effects of structural parameters of rough fractal tree-like fracture network on permeability were analyzed. Based on these analyses, the main conclusions are drawn as follows:

Firstly, this analytical permeability model can well describe the branching characteristics and self-affine surface roughness of branch fractures in the fractal tree-like fracture network. It can be used to estimate the permeability of complex fracture network in shale gas reservoirs.

Secondly, the mechanisms for the energy loss and gas flow resistance in fractures vary with fracture surface roughness and aperture scale. At millimeter scale, the gas flow velocity in fractures is much larger than that at micron scale. Eddy flow and backflow are observed in the rough fractal tree-like fracture network at millimeter scale, while no eddy flow occurs at micron scale. Higher gas flow velocity is more likely disturbed by local rough surfaces.

Thirdly, the permeability KF of rough fractal tree-like fracture network decreases with the increase of the branch level m, branch angle θ and length ratio α, but increases with the increases of aperture ratio β. Larger aperture ratio β and smaller length ratio α will produce lower gas flow resistance and larger gas flow channel.