Microfractures with different properties in rocks interconnect to form complex microfracture networks^{[1,2,3,4,5,6]}. The microfracture network improves the seepage capacity of the rock to varying degrees, which makes the permeability of the rock under the same porosity vary greatly, and the permeability and porosity of the reservoir are poorly correlated, which increases the difficulty in the study on reserve mobility and seepage capacity^[7,8]. It is necessary to propose a new method to evaluate the overall property of microfracture network.

Microfractures in rocks are mainly identified by rock thin section, cast thin section and scanning electron microscopy^[2]. Scholars have characterized the basic features of single microfracture by using fracture length^[9], equivalent fracture aperture^[10], fracture surface roughness^[11,12], fracture tortuosity^[12,13] and dip angle^[9]. How-ever, it is difficult to get an overall understanding of microfracture network from microscopic observation results. With the development of micro-CT and digital core technology, three-dimensional digital core can intuitively display the complex morphology of microfracture networks in rocks and evaluate microfracture networks qualitatively. But evaluating microfracture network quantitatively is much harder. As different microfractures have different dip angles, it is difficult to use a unified dip angle describing the overall strike of a microfracture network^[4,5,6]. Meanwhile, the connectivity of a fracture network is often defined as the average number of intersections on each fracture^{[9, 14]}, or the proportion of connected fracture network density in the overall fracture network density^[15]. The permeability of porous media with a locally connected fracture network is not necessarily large^[16]. There still lack of effective means to quantitatively characterize the properties of microfracture networks^{[2, 9, 16-17]}.

The numerical simulation method based on finite difference is applicable to flow simulation at macro-scale^[18], but cannot meet the demand of flow simulation in microfracture network. Microscopic numerical simulation methods mainly include pore network simulation^[19], computational fluid dynamics (CFD)^[20] and lattice Boltzmann method (LBM)^[17]. Among them, LBM based on micro-kinetics, with clear background of microscopic particles, can deal with irregular geometry boundary and has great advantages in studying coupled flow in irregular micro-cracks and pores^{[17, 21-23]}. Gao et al.^[24] and Li et al.^[25] used LBM to simulate fluid flow in porous media and proved its reliability. The multi-relaxation-time lattice Boltzmann method (MRT-LBM) overcomes the deficiency that the simulation precision is affected by relaxation time of single relaxation time lattice Boltzmann method (SRT-LBM), and effectively improves the accuracy of simulation results^[26].

Based on the understanding of micro- fracture and pore structure in rocks, a method for quantitative evaluation of the connectivity and strike of microfracture network is proposed. A series of fractured porous media and fracture-pore porous media models are constructed by numerical algorithm. MRT-LBM is applied to simulate fluid flow in them, to reveal the influence mechanism of microfracture network on rock permeability.

Fig. 1 shows some thin section, scanning electron microscopy and micro-CT scanning images of tight sandstone and carbonate samples from Sichuan Basin. Fig. 1a and Fig. 1b are thin section images of tight sandstone samples. Fig. 1c and Fig. 1d are thin section images of carbonate samples. Fig. 1e-Fig. 1h are scanning electron microscopy images of tight sandstone samples. Fig. 1i and Fig. 1j are micro-CT images of tight sandstone samples at the resolution of 0.62 μm. Fig. 1k and Fig. 1l are micro-CT images of carbonate samples at the resolution of 14.10 μm. The porosity and permeability of tight sandstone samples are 4.1%-7.3% and 0.05×10^-3-0.14×10^-3 μm², respectively. The porosities of carbonate samples are 3.2%-7.5%, while their permeabilities are 0.03×10^-3-1.5×10^-3 μm². Fig. 1 shows that microfractures with complex morphology exist in both tight sandstone and carbonate rocks, including structural fracture, intergranular fracture, interlayer fracture, and dissolution fracture, etc.; some of the microfractures intersect to form microfracture network. Among them, microfractures in the tight sandstone samples are 0.002 mm to 0.05 mm in aperture, and microfractures in the carbonate samples are 0.002 mm to 0.100 mm in aperture, which are similar to the findings of Wang et al.^[27].

Micro-images and 3D digital cores can show the mor-phology of microfracture network in rocks, and the aperture and number of microfractures can be further obtained. However, due to the limitation of observation scale, some microfractures extend beyond the thin section, making it difficult to accurately obtain lengths of microfractures^[2]. The fracture length often affects the connectivity of microfracture network, which plays an important role in seepage capacity of rock. Moreover, 3D digital core fails to capture most of the matrix pores in rock due to the limited resolution. If simulation is carried out directly based on digital core, the connectivity of microfracture network in rock would be overestimated and the seepage capacity of matrix would be underestimated.

The high pressure mercury injection experiment can accurately measure properties of micro pores in rock. Taking the reservoir rock of a gas reservoir in Sichuan Basin as an example, its matrix pore radii measured are mainly 0.1-0.7 μm and approximately meet the normal distribution (Fig. 2). Numerical algorithms were used to construct fractured porous media model and fracture-pore media model based on the results of digital core and high pressure mercury injection experiment. Among them, the fractured porous media model (called fractured model for short) contains penetrative fracture networks (ignoring the effect of matrix pores), and the fracture-pore media model (called fractured-pore model for short) contains pores and non-penetrative fracture networks. The construction process of the fracture-pore model is as follows: (1) Set the area of porous media and the whole area was set as solid phase. (2) Based on the results of high-pressure mercury injection experiment, the matrix pore radii were set to satisfy the normal distribution shown in Fig. 2. Then the node numbers of pores with different sizes were determined, and the radii of generated pores were randomly assigned. When the node numbers of pores generated reached the expected value, pore growth was stopped, and pore connectivity was ensured (Fig. 3a). (3) Based on the results of micro image and digital core, the microfractures were set with different apertures and lengths, which were randomly distributed in the porous media. (4) Different number of microfractures was set, steps (1)-(3) were repeated to generate fracture-pore models with different microfracture networks (Fig. 3b). Steps (3) and (4) were needed only to generate different fractured models.

The seepage capacity of rocks is only related to the interconnected pores and microfractures. The four-connected region labelling algorithm^[23] was adopted to further remove the isolated pores and isolated fractures in porous media. Based on morphology principle^[28], inserting straight lines in the porous media along x and y directions at certain line spacing, the ratio of average number of fractures encountered by unit length line was defined as strike factor of fracture network (A):

Fig. 4 shows six micro-crack networks composed of 10 microfractures with length of 400 grids and aperture of 20 grids. Set the line spacing at 5 grids and insert straight lines in the porous media along x and y directions to obtain n_x and n_y. According to Eq. (1), A values of 6 microfracture networks are obtained between 0.05 and 17.61, which means that microfracture networks with different strikes have different A values. When the number of microcracks dipping to y direction is more than that dipping to x direction, A is less than 1. When the number of microfractures dipping to x direction is more than that dipping to y direction, A is greater than 1. When the number of microfractures dipping to x direction equals to that dipping to y direction, A is 1. As the microfractures are all trending in y direction, A is the smallest (A=0.05), while A is the largest (A=17.61) when microfractures all trending in in x direction. With the increase of microfractures dipping to x direction, A increases. Therefore, A can quantitatively characterize the overall strike of microfracture network.

Microfracture networks will form new flow channels in porous media, and microfractures may be directly connected or indirectly connected through matrix pores with each other. Previous evaluation methods for connectivity of fracture network fail to take into account the relationship between local fracture network and whole porous media, leading to a poor applicability^[15]. In this paper, the ratio of the maximum horizontal length (L_f,max) of microfracture connected to matrix pores at a flow path to the porous media length (L) is defined as connectivity coefficient (f ) of the microfracture network, representing the maximum proportion of microfracture on the fluid flowing path. f is between 0 and 1. The larger the f, the better the overall connectivity of the microfracture network in porous media. When f is less than 1, the microfracture network is inside the porous media, and the porous media is fracture-pore model. When f equals 1, the microfracture network connects the inlet and outlet ends of the porous media, and the porous media is called fractured model. Applying numerical algorithm to screen fractures connected with pores in the porous media and distinguish the connectivity between fractures, L_f,max can be obtained (Fig. 5, red circles are the starting point and ending point selected). f of Fig. 5a and Fig. 5b calculated with Eq. (2) are 0.58 and 0.70, respectively.

MRT-LBM is directly derived from Boltzmann equation, and can be expressed as^[23]:

where Ω=M^-1SM

M is the transformation matrix of distribution function f_i^[29]. Eq. (3) can be written as:

where m=Mf_i m_eq=Mf_eq,i

In the two-dimension and nine-discrete-velocity-direction model (D2Q9), the equilibrium distribution function is:

where ${{c}_{s}}={c}/{\sqrt{3}}\;$

The discrete velocity c_i is defined as: c₀=(0, 0), c₁=-c₃=(c, 0), c₂=-c₄= (0, c), c₅=-c₇=(c, c), c₆=-c₈= (-c, c). The diagonal matrix S is defined as:

In this paper, s_ρ=s_j=1.0, s_e=1.19, s_ε=1.4 and s_q=8(2-s_v)/ (8-s_v)^[29]. Among them, s_v is related to the kinematic viscosity. The density, velocity and kinematic viscosity of macroscopic fluid can be expressed as Eq. (7).

The pore surface of porous media is coarse, and classical LBM boundary conditions can be directly applied^[28]. In this paper, pressure boundary is adopted at the inlet and outlet of the porous media^{[25, 30]}, and standard bound boundary is adopted at the upper and lower boundary of the porous media.

The simulation program of MRT-LBM was written with Fortran language, and its accuracy was verified in simulating fluid flow in porous media.

2.2.1. Simulation of fluid flow in the fractured model

Combined with the apertures, lengths and numbers of microfractures in reservoir rocks, 28 fractured models were constructed with numerical algorithm. MRT-LBM was used to simulate fluid (water) flow in porous media with microfracture apertures of 30, 60 and 120 μm, respectively. The construction parameters of the fractured models and detailed simulation conditions are shown in Table 1. The fluid flow in the simulation was laminar (with small Re). The simulation deviation was calculated by Eq. (8). When the simulation deviation was less than 10^-6, it is considered that simulation results wouldn’t change anymore and the iterative process was ended^[30,31].

Fig. 6 shows the simulation results of the fractured model with microfracture aperture of 60 μm, which indicates that fracture network significantly affects fluid flow path and velocity field distribution in the porous media. Based on the simulation results, the microfractures in the porous media not participating in the flow were further deleted to obtain connected porosity of the fractured model (i.e. the ratio of interconnected pores volume to the total volume of the porous media). The line spacing was set at 5 grids, and Eq. (1) was adopted to calculate strike factors of fracture networks in the fractured models. The strike factors of fracture networks in the fractured models calculated are between 0.99 and 2.70, and the connected porosities of the models are between 2.3% and 5.0%, showing there is no obvious correlation between the strike factor and connected porosity (Fig. 7).

2.2.2. Simulation of fluid flow in the fracture-pore models

Combined with the reservoir rock properties, the matrix porosity of the fracture-pore model was set at 5%, and the pore radii meet normal distribution (Fig. 2), while the porosity of microfractures was set at 5%. Considering that microfractures in rocks are generally larger than matrix pores, microfractures were set at 2 μm and 4 μm in aperture, respectively. The microfractures were randomly distributed within the fracture-pore models, and their strikes and lengths were randomly set. The starting and ending points of the microfractures were connected with matrix pores. Detailed construction parameters are listed in Table 2. To consider the effect of randomness, 20 fracture-pore models were generated at each microfracture aperture, and a total of 40 fracture-pore models were generated.

Connectivity coefficients of fracture networks in the fracture-pore models were calculated by Eq. (2). The D2Q9 MRT-LBM was used to simulate fluid (water) flow in the models, and the detailed simulation parameters are shown in Table 2. Fig. 8 shows the simulation results of the matrix models and fracture-pore models with fracture network connectivity coefficients of 0.36, 0.70 and 0.80 (fracture aperture of 2 μm). It shows that the fracture network significantly affects fluid flow path and velocity field distribution in fracture-pore model, the fluid in some flow channels increases significantly in velocity (Fig. 8c, 8d), and the dominant channel effect is significant. Then, microfractures not involved in the flow were deleted. The line spacing was set at 5 grids, and strike factors of fracture networks in the fracture- pore models were calculated with Eq. (1). Figs. 9 and 10 show the relationships of connected porosity versus strike factor and connectivity coefficient of fracture network. At the fracture aperture of 2 μm, the connected porosities of the fracture-pore models are between 8.4% and 9.5%, the strike factors between 0.99 and 1.58, and the connectivity coefficients between 0.11 and 0.80. At the fracture aperture of 4 μm, the connected porosities of the frac-ture-pore models are between 8.6% and 9.2%, the strike factors of fractured networks between 0.89 and 3.63, and the connectivity coefficients between 0.34 and 0.70. Clearly, there is no significant correlation between the connected porosity and fracture network strike factor and fracture network connectivity coefficient of fracture- pore model.

Permeability is often used to characterize seepage capacity of porous media, while tortuosity reflects the tortuous degree of flow paths in porous media which is an important parameter affecting permeability of porous media^[25]. Based on MRT-LBM simulation results, Eq. (9) and Eq. (10) were used to calculate the permeability and tortuosity of the fractured model and fracture-pore model^[28]. Figs. 11 and 12 show that there is no signficant correlation between permeability and connected porosity of two kind of models. The matrices of the fracture-pore models are same in porosity and pore radius distribution, which indicates that the permeability is mainly affected by the microfracture network. The relationship between permeability of porous media and fracture network properties needs to be further analyzed to find out the influence mechanism of microfracture network on the permeability of porous media.

The permeability and strike factor of fracture network of the fractured model are positively correlated (Fig. 13). When the strike factor of fracture network increases from 1.0 to 2.7, all the fractured models with different fracture apertures increase by more than 160% in permeability. The tortuosity is negatively correlated with strike factor of fracture network (Fig. 14), and the polynomials of them have the best fitting, with a correlation coefficient square of 0.953 (Table 3). In contrast, there are not significant correlations between the permeability and tortuosity and fracture network strike factor of the fracture-pore model (Figs. 15 and 16). Further analysis shows that the greater the strike factor of microfracture network, the greater the fluid velocity at the outlet of porous media will be (Fig. 17). From the microscopic point of view, the larger the strike factor of fracture network, the smaller the tortuosity of the fractured model, the flatter the flow path, and the greater the permeability of the fractured model will be. However, for the fracture-pore model (f<1), the matrix pores and microfractures generally have magnitude difference in size, so the matrix pores are a key factor affecting the flow capacity of fracture-pore model^[32]. Microfracture network strike has much stronger impact on the flow capacity of fractured model than that of fracture-pore model.

The permeability of fracture-pore model is positively correlated with the connectivity coefficient of fracture network (Fig. 18). At the fracture aperture of 2 μm, when the connectivity coefficient of fracture network increases from 0.11 to 0.80, the permeability of fracture-pore model increases by 61.83%. At the fracture aperture of 4 μm, when the connectivity coefficient of fracture network increases from 0.34 to 0.70, the permeability of fracture- pore model increases by 48.43%. At the same time, the exponential relation between permeability and connec- tivity coefficient of fracture network of fracture-pore model has the best fitting effect (taking the fracture aperture of 2 μm as an example), with a fitted correlation coefficient square of 0.832 (Table 3).

Subsequently, the relationships between connectivity coefficient of fracture network and the sum of fluid velocities of flow channels at the outlet and the flow proportion of different channels (that is, the ratio of fluid flow through one channel to the total fluid flow through porous media) were analyzed (Fig. 19). With the increase of connectivity coefficient of fracture network, the sum of fluid velocities of all flow channels at the outlet of fracture-pore model generally increases, and the local difference in fluid velocity is caused by the randomness of structure of the porous media. The flow proportion of each channel changes with connectivity coefficient of fracture network, but there is no obvious correlation between them. This means the greater the connectivity coefficient of fracture network, the greater the proportion of fractures in the fluid flow path will be, the less energy the fluid will consume when flowing through the model, and the greater the permeability of fracture-pore model will be. Therefore, the connectivity coefficient of fracture network has a stronger influence on the seepage capacity of fracture-pore model than the strike of fracture network.

In addition, the dominant channel effect is significant in some fracture-pore models. For example, at the connectivity coefficient of fracture network of 0.70 and 0.80, the flow proportion of channel 2 is over 50%, and at the connectivity coefficient of fracture network of 0.36, the flow proportion of channel 4 is over 50% (Fig. 19). In general, the greater the connectivity coefficient of fracture network, the more likely the dominant channel effect will occur in porous media (especially for f≥0.7). When dominant channel effect occur, the flow proportion of the remaining channels in the porous media decreases, while the flow rates through the remaining channels mostly increase, and the flow rates through individual remaining channels decrease in a few cases.

The permeability of the fractured models with fracture aperture of 60 μm and 120 μm are about 4 times and 16 times than that of the fractured model with fracture aperture of 30 μm, respectively. The correlations between strike factor of fracture network and connected porosity, permeability and tortuosity of the models with fracture aperture of 60 μm and 120 μm are in consistent with those of the fractured model with fracture aperture of 30 μm (Figs. 13 and 14). Namely, the seepage capacity of fractured model mainly depends on the strike of fracture network and fracture aperture, and the effect of fracture aperture is more significant. The characteristic length of matrix pores in fracture-pore model is 0.344 μm^[31]. When the ratio of fracture aperture to characteristic length of matrix increases from 5.814 to 11.628 (that is the fracture aperture increase from 2 μm to 4 μm), the permeability of fracture-pore model changes less than 10% (Fig. 18). The fracture aperture has much smaller impact on seepage capacity of fracture-pore model than fracture network connectivity.

At present, the evaluation methods of reserve mobility mainly include experimental method and numerical simulation method. The method of calculating the lower limit of porosity for movable reserves by considering economic factors, experimental results and field test data has got good results. But it can be seen from the above analysis that there is no significant correlation between seepage capacity and porosity of porous media with microfracture network. Thus, the method of defining movable lower limit of reservoir reserves by porosity does not consider the influence of microfracture network on rock seepage capacity. Taking the exponential relationship between porous media permeability and microfracture network connectivity obtained in this work as an example, when the matrix porosity and fracture aperture are consistant, as the connectivity coefficient of fracture network increases from 0 to 1, the permeability of porous media increases by 120.98%, and the reserve mobility changes significantly. Therefore, it is necessary to analyze the microfracture properties (such as fracture aperture, number and connectivity, etc.) of rocks through microscopic experiments to evaluate the influence of microfracture network on rock permeability for different gas reservoirs. For gas reservoirs with weak influence of microfracture network, the porosity lower limit obtained by experimental method can be directly used to evaluate reserve mobility. For gas reservoirs with significant influence of microfracture network and weak water energy, the porosity lower limit obtained by the experimental method should be lowered according to different influence degrees of the fracture network to improve the accuracy of reservoir evaluation. For gas reservoirs with significant influence of microfracture network and strong water energy, the existence of dominant channel effect will lead to a decrease in sweep efficiency during waterflooding in rocks^[33], aggravating the water locking of gas, and reducing the mobility of reservoir reserves.

Accurate reservoir permeability is the key to the evaluation of gas well productivity. The reservoir permeability in near-well area and far-well area obtained from well test analysis often differ widely, showing strong reservoir heterogeneity. It can be seen from this work that the difference of reservoir flow capacity in near-well area and far-well area may be caused by different fracture network strikes. At the same microfracture aperture and microfracture number, when the strike factor of fracture network increases from 1.00 to 2.64, the fractured model decreases by more than 25% in tortuosity and increases by more than 160% in permeability. Considering the fracture network strike can further improve the prediction of gas well productivity.

The strike factor and connectivity coefficient of microfracture network can quantitatively characterize the overall strike and connectivity of microfracture network in porous media. With the increase of fracture network strike factor, the fractured porous media decreases significantly in tortuosity and increases gradually in permeability, and the tortuosity and strike factor of fracture network satisfies polynomial relation. With the increase of fracture network connectivity coefficient, the fracture- pore model is more likely to have dominant channel effect and increases in permeability, and the permeability and fracture network connectivity coefficient of the model satisfy the exponential relation. This research results can be used to calculate permeability tensor of rock samples with different fracture development conditions, to improve the evaluation accuracy of reservoir reserve mobility and gas well productivity.

A—strike factor of fracture network, dimensionless;

c—lattice speed, m/s;

c_i—discrete speed, m/s;

c_s—sound speed, m/s;

e_i—lattice velocity, m/s;

f—connectivity coefficient of fracture network, dimensionless;

f_i—discrete distribution function, kg/m³;

f_eq,i—equilibrium distribution function, kg/m³;

i—discrete direction;

K₀—model permeability, 10^-3 μm²;

K_e—matrix permeability of fracture-pore model, 10^-3 μm²;

L—length of porous media, m;

L_fn—horizontal length of a microfracture on the flow path, m;

L_f,max—the maximum value of the sum of horizontal lengths of microfracture networks on the flow path, m;

m_i—function of the moment space;

m_eq,i—equilibrium function of the moment space;

M—transformation matrix;

n—total number of fractures on flow channels;

n_x, n_y—average number of fractures encountered by a line per unit length as the line crosses the porous media in the x or y direction;

R—correlation coefficient, dimensionless;

Re—Reynolds number, dimensionless;

s_e, s_j, s_q, s_v, s_ε, s_ρ—energy moment, momentum moment, energy flux moment, stress tensor moment, energy square moment, and density moment, dimensionless;

S—diagonal matrix;

t—time, s;

u—fluid velocity, m/s;

u(x, y)—fluid velocity at point (x, y) in porous media, m/s;

u(x, y, t)—fluid velocity at point (x, y) in porous media at time t, m/s;

u_x(x, y)—fluid velocity at point (x, y) in x direction, m/s;

${{\bar{u}}_{x}}$—average fluid velocity in the x direction, m/s;

w_i—weight coefficient, dimensionless;

x—position vector of porous media;

x, y—two directions in the rectangular coordinate;

δ_t—time step, s;

δ_E—relative deviation of simulation result, dimensionless;

p—pressure gradient, Pa/m;

μ—dynamic viscosity, Pa·s;

v—kinematic viscosity, m²/s;

ρ—fluid density, kg/m³;

τ—tortuosity, dimensionless;

Ω—collision matrix;

ϕ—total porosity of porous media, %;

ϕ_eff—connected porosity of porous media, %.