A Review of Flow Mechanism and Inversion Methods of Fracture Network in Shale Gas Reservoirs

The pore structure of shale gas reservoirs has strong heterogeneity, and the flow mechanism in multiscale media is complex. The fracture network of hydraulic fracturing is significantly affected by reservoir in situ stress, rock mechanical properties, and natural fracture distribution. At present, there is no efficient and accurate inversion method for fracture networks. Accurately describing fracture network morphology and flow capacity distribution of induced fracture is an important basis for production analysis, fracturing evaluation, and production plan. This article focuses on the hot issues of shale gas development, from three aspects: flow parameter characterization method of organic/inorganic matter, multiscale mass transfer simulation of shale gas reservoir, and inversion method of fracture network morphology and flow capacity, to introduce relevant research progress in detail. At the same time, the advantages and shortcomings of current related researches are compared and analyzed. Based on this, the key scientific problems existing on flow mechanisms and inversion method of fracture network in shale gas reservoirs are proposed, which can provide guidance for further research.


Introduction
Stable economic growth has led to an increase in demand for oil and gas, and conventional oil and gas resources have been unable to meet the large energy demand. The shale gas resources are widely distributed and rich in reserves, which will become an important area for strategic succession and development of oil and gas resources, and the development prospects are optimistic. There are large differences between shale gas reservoirs and conventional gas reservoirs. Shale reservoirs have ultralow porosity/ultralow permeability, the matrix permeability is generally in the order of nano-Darcy, and the pore size is much smaller than that of sandstone and limestone [1]. The rock pore types are diverse, and the spatial structure is complex [2]. The gas occurrence is com-plex, mainly including dissolved gas, adsorbed gas, and free gas. Therefore, shale gas reservoirs are difficult to obtain economic production using conventional development strategy. In recent years, the development of horizontal wells and hydraulic fracturing technology has made shale gas reservoirs have development value under the current oil price and development technology [3]. Hydraulic fracturing not only produces primary fractures with high conductivity but also connects natural fractures to form a complex network of induced fractures [4,5] (Figure 1), which greatly enhances the original reservoir permeability and has a decisive impact on production. Therefore, how to accurately calculate the flows parameters of shale gas reservoirs, characterize the fracture networks morphology and flow capacity, and achieve an accurate prediction of development performance is an important prerequisite for the efficient development of shale gas reservoirs.
However, the distribution of fracture networks of shale gas reservoirs is extremely complex. The existing fracture propagation simulation and fracture morphology inversion methods have problems, such as low efficiency, poor accuracy, and inability to accurately match the actual fracture network morphology [6]. At the same time, due to the unusually complex spatial structure of matrix/natural fractures in shale gas reservoir and the large difference in fluid flow between different-scale media, it is more difficult to obtain the flow capacity of various induced fractures through dynamic fitting. Based on the investigation and analysis of inversion methods of induced fracture network, this paper proposes a new idea of fracture characterization based on the combination of fracture propagation calculation and multiscale flow simulation. That is, on the basis of revealing the influence of the micropore structure of shale gas reservoirs on gas storage and migration capabilities, this paper proposes a new research idea for the flow characterization and morphology inversion method of induced fracture network. The proposed idea can clarify the fracture network morphology, flow ability distribution of heterogeneous fracture network and matrix. The new approach can provide guidance of realizing accurate prediction and optimization of shale gas reservoir. Therefore, this paper discusses the current research progress from three aspects: Firstly, flow parameters characterization of porous media in shale gas reservoirs are discussed, and the current common methods for calculation porosity/permeability of inorganic/organic matter in shale reservoirs are introduced. Secondly, we will focus on the multimedia mass transfer simulation of shale gas reservoirs and summarize its current main research methods and progress. Finally, we will explore the research progress of inversion methods of induced fracture network of shale gas reservoirs.

Flow Parameters Calculation of Porous Media in Shale Reservoirs
The core displacement experiment is a common method of intuitively measuring matrix porosity/permeability and other flow parameters [8], including mercury intrusion method, nuclear magnetic resonance method, pressure pulse attenuation tests, pressure recovery method, and desorption flow method [9,10]. Based on laboratory experiments, the gas storage spaces and flow capacity under the influence of real pore space structure of large-scale cores can be considered. However, on the one hand, hydrocarbons in organic-rich shale have various occurrence types (dissolved gas, adsorbed gas, and free gas), and the gas density of different pore spaces varies greatly. The gas migration mechanism is complex (slip, adsorption and desorption, Knudsen diffusion, surface diffusion, etc.), and it is greatly affected by the temperature, pressure, and properties of pore surface [11,12]. Therefore, the permeability of shales is about three orders for magnitude higher than Darcy's permeability. Conventional porosity/permeability cannot fully reflect the occurrence and migration capacity of shale gas [13] ( Figure 2). On the other hand, the shale has nanoscale pores and poor fluid flow capacity. It takes a long time to measure with displacement experiments. The experimental instruments are not sensitive to small pressure drops or flow differences in the measurement process. It is greatly affected by external pressure, temperature, and other factors.
In view of the complex occurrence and migration mechanisms of fluids in shale gas reservoirs, Javadpour (2009) proposed an apparent permeability model for the first time to describe the multiple transports of shale gas in micronanopores. On this basis, Sheng et al. (2018) established apparent porosity model, which comprehensively considers the difference of adsorbed gas/free gases concentration and then couples the different gas occurrence of shale pores. The apparent porosity/permeability model combined with Darcy's equation can apply conventional numerical simulation methods into dynamic analysis of shale reservoirs, which has been widely used in flow simulation. Based on the microscale single-tube apparent permeability proposed by Javadpour, Akkutlu et al. [14], Civan et al. [15], Darabi et al. [16], Wu et al. [17], Mi et al. [18], Zhang et al. [19], and Sheng et al. [20] consider the complex structure of porous media and gas migration mechanism and proposed the analytical apparent permeability model for shale cores [21][22][23] ( Table 1). However, the current analytical apparent porosity/permeability models cannot consider the real pore structure, especially the three-dimensional spatial structure of pores.
At the same time, scholars proposed the approach for calculation apparent permeability, in which complex pores spaces modeling by digital core and numerical simulations are used [27]. The CT scanning or theoretical simulation is used to build a real three-dimensional pore network model and solved through numerical simulation methods [24,[26][27][28] (Figure 3). This method can consider the complex pore space distribution of shale, but limited by computational resources. It has the disadvantages of complex modeling, long time-consuming, and small simulation size. At the same time, the dynamic change of pore size resulting from stress sensitivity, organic matter desorption shrinkage, and clay minerals swelling with water during depressurization production has a large impact on fluid migration [8]. It is difficult to use numerical simulation methods to characterize these characteristics.
In 2020, Sheng et al. propose a new approach to describe the complex migration mechanism, in which the digital cores are used to describe the pore structure, and analytical models were used to simulate gas flow mechanism ( Figure 4) [30]. Combining the planar pore structure obtained by CT scanning, the pore structure of porous media is comprehensively considered, and the dynamic change of pore size (stress sensitivity, desorption shrinkage), the apparent permeability model of organic/inorganic porous media coupled with microscale migration mechanism, and the apparent porosity model coupled gases with different occurrence states are proposed. However, the existing analytical apparent porosity/permeability model based on CT scanning pore network cannot well reflect the three-dimensional spatial distribution of shales. Therefore, it is necessary to establish a new method to couple the actual complex pore space structure and the analytical method of gas occurrence and migration capabilities, to realize the analytical characterization of the threedimensional apparent porosity/permeability of shale gas reservoirs.

Mass Transfer Simulation in Multiscale Media
At present, commercial numerical simulators mainly use multicontinuous media models and discrete fracture models to couple fluid flow in the matrix and fracture network. Multicontinuous media model is a common method for flow simulation in matrix and fracture networks [31,32]. According to the fluid exchange between media in different-scale media, scholars proposed pseudosteady/transient cross-flow dual porosity models and dual permeability model [33,34]. The pore size of organic matter and inorganic matter in shale gas reservoirs is quite different. Due to the influence of factors such as slip effect, Knudsen diffusion, and surface diffusion, the permeability of shale is about three orders of  3 Geofluids magnitude higher than the Darcy permeability, which requires the dual permeability model to describe cross-flow between inorganic and organic matter in shale reservoirs [35]. Considering the multicontinuous media distributions in shale reservoirs, including organic matter, inorganic matter, and fracture network, scholars proposed the multicontinuous model to describe gas migration in shale gas reservoirs [35,36]. Since fracture network morphology and flow capacity have a greater impact on flow behavior, based on the conventional diffusion equation (uniform induced fracture spacing and uniform fracture network porosity/permeability, CDE), Cinco-Ley et al. [37], Reis [38], Ranjbar et al. [39], and other scholars proposed the modified conventional diffusion equations (MCDE) to simulate fluid flow between fractures and matrix. This equation takes into account the heteroge-neously distribution of induced fracture spacing and the uniformly distribution of fracture network porosity/permeability. Normal distribution and linear and exponential functions are usually used in MCDE to describe the distribution of induced fracture spacing. The results show that the fracture spacing distribution is an important parameter that affects fluid migration in matrix blocks during production.
At the same time, scholars also proposed the fractal diffusion equation (FDE), which takes into account the uniformly distribution of induced fracture spacing and the fractal distribution of fracture network porosity/permeability [40][41][42]. Cossio et al. [42] combined FDE with the trilinear flow model and proposed a semianalytical solution for vertical fracture flow with limited conductivity. Wang et al. [43] proposed a Simple, only suitable for straight pipes, ideal gas, without considering adsorption and desorption Civan [15] Second-order slip model, consider slippage effect Higher-order slip flow, contains empirical parameters Darabi et al. [16] Based on the hypothesis of slippage effect, Maxwell's theory is adopted, considering the surface roughness and Knudsen diffusion in porous media Consider the tortuosity and surface roughness of porous media, ideal gas, without considering adsorption and desorption Akkutlu et al. [12] Based on matrix-fracture dual media model, the matrix contains organic and inorganic matter, considering surface diffusion and adsorption Dual media system; complex numerical model Shabro et al. [24] Construct porous structure based on finite difference numerical model and geometric parameters Characterization of spatial structure of porous media, complex numerical models, ignoring adsorption Mehmani et al. [25] Apply Javadpour's model to pore network model Characterization of spatial structure of porous media, complex numerical models, ideal gases, ignoring adsorption Chen et al. [26] Based on SEM to construct the pore structure, considering Knudsen diffusion and slippage effects Characterization of spatial structure of porous media, complex numerical models, ignoring adsorption, and containing empirical parameters Wu et al. [17] Consider influence of real gas properties and stress on permeability and pore size Consider complex structure of pore network and the dynamic change of pore size under stress and desorption shrinkage; fail to consider the complex spatial structure Zhang et al. [19] Divide flow patterns based on Knudsen number, consider multiple migration mechanisms Considering multiple transport mechanisms of microscale fluids, failing to consider complex spatial structures Mi et al. [18] Considering matrix porous media and natural fractures, the series-parallel model is used for coupling permeability Considering microscale gas slippage effect, diffusion effect, and natural fracture flow characteristics, failing to consider complex spatial structure Sheng et al. [20] Based on molecular collision theory, considering the influence of adsorption layer migration on slip coefficient of free gas Analytical methods are used to realize the combination of multiple migration mechanisms, complex pore structure features and pore size dynamic changes, fail to consider complex spatial structures  4 Geofluids semianalytical flow model for multistage fractured horizontal wells by applying FDE in SRV and analyzed the influence of fractal fracture network porosity/permeability on pressure response. Relevant research results show that the heterogeneous distribution of fracture network porosity/permeability has a great influence on bottom hole pressure response. In the multicontinuous model, the fracture network morphology and flow capacity distribution are simple, while shale gas reservoirs have strong heterogeneity in fracture network morphology and flow capacity. Therefore, using the multicontinuous model to equivalently consider the fracture network heterogeneity will lead to larger errors. On the basis of quantitative understanding of fracture distribution, the discrete fracture model (DFM) uses unstructured grids to characterize the fracture network morphology and obtains the flow behavior by numerical methods [45][46][47][48]. DFM is widely used in the numerical simulator for heterogeneous reservoirs [49][50][51]. Compared with multicontinuous models, it can more accurately describe the fluid flow in complex fracture network. The embedded discrete fracture model (EDFM) was developed based on the DFM model, making the DFM simple for reservoir simulator and obtaining higher computational efficiency [52][53][54][55]. Feng et al. [48] established the apparent permeability model of shale gas reservoirs based on EDFM and studied the influence of natural fracture distribution on gas flow. AlTwaijri et al. [56] used EDFM to simulate complex fracture geometry and proposed a 3D flow numerical model for shale gas reservoirs. The mesh refinement method based on EDFM or DFM can describe fracture network morphology accurately; however, it cannot reflect the flow difference between organic and inorganic shale. On this basis, Sheng et al. proposed the dual fractal dif-fusion equation [44] ( Figure 5) and a heterogeneous discrete fracture coupling multicontinuous model (Figure 6), which comprehensively considers the heterogeneity distribution of induced fracture spacing and fracture network porosity/permeability. The model was verified with production data in Barnett reservoirs, and the dynamic matching accuracy was significantly improved. However, the complex pore space structure of shale has a great influence on fluid flow, and the relevant research of numerical simulation has not been fully considered. Therefore, it is necessary to comprehensively consider the discrete fracture and the complex flow mechanism between dual media of organic and inorganic, to realize the numerical simulation of mass transfer in multimedia of shale gas reservoir.

Inversion Methods for Fracture Network in Shale Gas Reservoirs
Accurate description of fracture network morphology and flow capacity is an important basis for production performance analysis, fracturing evaluation, and production plan formulation. The distribution of natural fractures in shale gas reservoirs is complex, and reservoir heterogeneity and in situ stress distribution make it more difficult to describe induced fractures. Laboratory experiments restore the real physical properties of rock, and it is also the main method to study fracture morphology [57]. Zhang et al. [58] carried out a series of laboratory experiments about shale reservoir fracture propagation, using the acoustic emission monitoring system to monitor the generation and expansion process of shale fracturing fractures and to observe the shape of hydraulic fracturing fractures. According to Li et al. [59], the 3D

Geofluids
fracture network is reconstructed based on serial-sectioned digital images of a shale cube hydraulically fractured under true-triaxial conditions. However, experimental simulation can construct a real reservoir environment to simulate the fracture morphology, but limited by the experimental environment, the simulation size is different from the actual size, so the combination of computer simulation can better solve the problem. [60]. The inversion of fracture network from    Figure 6: Heterogeneous discrete fractures coupled with dual media of organic and inorganic matter in shale gas reservoirs. 6 Geofluids microseismical data can avoid the shortcomings of existing fracture prediction methods, such as physical experiments and fracture propagation simulation, which are inefficient in calculation and difficult to be applied to actual reservoirs. Fisher et al. [61] used the microseismic fracture monitoring technology to evaluate the stimulated reservoir volume of the Barnett shale horizontal wells. Based on the microseismic image, Warpinski et al. [62] found that the fractures developed along the orthogonal direction, forming a complex fracture network. Scholars also found that the effective stimulated reservoir volume is often only 15%~25% of the microseismic monitoring volume, and the fracture network shape has a decisive impact on the effective stimulated reservoir volume [63]. At the same time, hydraulic fracturing experiments show that the fractal fracture propagation model can better match the fracture network morphology from fracturing experiment [64]. The fracture in different scales can be characterized by a fractal model [65,66]. Yu [67] proposed the self-similar fractal network model and discussed the influence of fractal bifurcation network on fluid flow. Cai et al. [68] quantitatively characterized the fracture distribution by fractal model and analyzed the dynamic mechanism of hydraulic fracture propagation. The fractal fracture morphology reflects the mechanical mechanism of fracture propagation to some extent. Therefore, scholars have proposed the inversion method of fracture network by describing fracture growth rules by fractal theory and matching with microseismic data [65,66]. Zhang et al. [69] used the principle of the shortest propagation path of seismic wave, and the equivalent ray tracing method based on Snell law to analyze microseismic data, inversed the complex fracture network morphology, and evaluated the fractal parameters of fracture network by combining with history matching. Zhou et al. [70] proposed a new induced fracture inversion method by matching microseismic data with fracture propagation rules based on fractal geometry theory and further applied the fractal bifurcation theory to fracture morphology inversion. On this basis, Sheng et al. established a random fractal fracture inversion method for fractured horizontal wells based on microseismic data and fractal bifurcation fractures. However, the fractal fracture propagation algorithm based on microseismic data does not take into account the effects of reservoir stress, rock mechanics parameters, and natural fracture heterogeneity on the fracture propagation. The accuracy of the inversion fracture network morphology needs further verification. Historical matching of production data can further to constrain fracture network morphology and flow capacity.
(1) In the aspect of analytical semi-analytical model, scholars used the dual media model to equivalently describe fracture network, then established the flow models for fractured horizontal wells. Combined with unstable pressure response of actual wells, the fracture parameters such as fracture length, fracture conductivity, and fracture network width are inversed [72,73]. To further improved the reliability of parameter inversion, Chen et al. [74] propose a dynamic inversion method by comprehensive analysis microseismic data, production data, and well test data, which further reduced the multisolution of fracture parameter inversion.
(2) In terms of numerical simulation, based on microseismic data and production, Mayerhofer et al. [4] used the historical matching method to evaluate fracture network volume, fracture spacing, and fracture length. Based on the numerical fracture model with automatic gridding technology and production dynamic data, Cipolla et al. [75] establish a fracture parameter inversion method which comprehensively considered microseismic data and production data. However, most of the above work is based on the assumption of a single fracture or orthogonal fractures (Figure 7), without considering the actual fracture network morphology. The inversion of complex fracture networks cannot be achieved. In recent years, automatic history matching technology has made important breakthroughs and developments and achieved large-scale parameter inversion. Zhao and others scholars have made progress in automatic inversion. Based on largescale parameter, dimensionality reduction algorithms, gradient-free algorithms, and other methods are established for automatic history matching. However, most of the methods currently do not combine actual monitoring data, and the accuracy of inversion fracture parameters needs to be further improved. Therefore, the complex morphology of fracture network in shale gas reservoirs brings serious ambiguity to the inversion of fracture flow capacity. It is urgent to establish a comprehensive inversion method by   7 Geofluids combining fracture growth algorithms, microseismic data, and production data history matching. Using the existing fracture network simulation method, the initial fracture morphology can be calculated through the actual geological parameters, and optimization algorithms such as SPSA optimization algorithm can be used to match the fracture morphology and microseismical data. The actual production data can also be combined to calculate the most realistic fracture morphology. The integrated fracture inversion method can reduce the ambiguity of the inversion and improves the inversion efficiency.

Conclusion
This paper introduced the research advance of flow mechanism and inversion methods of fracture network in shale gas reservoirs. The advantages and shortcomings of current researches about flow parameter characterization method of organic/inorganic matter, multiscale mass transfer simulation of shale gas reservoir, and inversion method of fracture network morphology and flow capacity were fully discussed. The following conclusion is obtained: (1) The laboratory experiment for shale porosity/permeability is time-consuming and not sensitive to small pressure drops or flow differences in the measurement process. Existing shale gas reservoir matrix apparent porosity/permeability models cannot fully consider the complex pore space distribution characteristics of shale (2) Traditional multicontinuous models cannot consider the flow difference between organic matter-inorganic matter-induced fractures in shale gas reservoirs and the influence of complex pore spaces structure on fluid flow (3) Accurate inversion of fracture networks morphology and flow capacity require comprehensive consideration of fracture growth algorithm, microseismic data, and production dynamic data. The algorithm should be combined with the actual reservoir parameters to reduce the multiplicity of inversion and improve the inversion efficiency.