A Review of Macroscopic Modeling for Shale Gas Production: Gas Flow Mechanisms, Multiscale Transport, and Solution Techniques

: The boost of shale gas production in the last decade has reformed worldwide energy structure. The macroscale modeling of shale gas production becomes particularly important as the economic development of such resources relies on the deployment of expensive hydraulic fracturing and the reasonable planning of well schedules. A ﬂood of literature was therefore published focused on accurately and efﬁciently simulating the production performance of shale gas and better accounting for the various geological features or ﬂow mechanisms that control shale gas transport. In this regard, this paper presents a holistic review of the macroscopic modeling of gas transport in shale. The review is carried out from three important points of view, which are the modeling of the gas ﬂow mechanisms, the representation of multiscale transport, and solution techniques for the mathematical models. Firstly, the importance of gas storage and ﬂow mechanisms in shale is discussed, and the various theoretical models used to characterize these effects in the continuum scale are introduced. Then, based on the intricate pore structure and various pore types of shale gas reservoirs, this review summarizes the multiple-porosity models in the literature to represent multiscale gas transport, and discusses the applicability of each model. Finally, the numerical and analytical/semi-analytical approaches used to solve the macroscopic mathematical model governing shale gas production are reviewed, with a focus on the treatment of the complex fracture network formed after multistage hydraulic fracturing.


Introduction
The rapid growth of shale gas exploitation in the last two decades has revolutionized the global energy landscape.It is commonly believed that the independence of US energy is attributed mainly to the vast development of shale gas reservoirs [1].Originally, shale was only known as a fine-grained sedimentary rock with no natural productivity due to its ultra-low intrinsic permeability, typically between 0.00001 to 0.001 mD [2,3].The advancements of horizontal drilling and hydraulic fracturing significantly enhanced well productivity by creating pathways for fluid flow, which removed the barriers hindering the economic development of such reserves and ignited the shale revolution [4,5].
The accurate modeling of shale gas production is crucial for the efficient and sustainable extraction of such resources.However, the traditional theories for fluid transport in macroporous materials are no longer applicable, considering the unique characteristics of shale.The main features of shale compared to conventional sandstone reservoirs are threefold.Firstly, the gas storage and transport mechanisms in shale are different from those of conventional gas reservoirs, e.g., shale contains a significant amount of adsorbed gas [6], and gas transport in the shale matrix is shown to deviate strongly from Darcy's law [7].Secondly, shale reservoirs exhibit significant heterogeneity in terms of mineralogy, organic content, and fracture distribution [8].This heterogeneity leads to variations in pore types (round pores, slit pores, fissures) and pore sizes (from nanometer kerogen pores to micrometer fractures) [9].Thirdly, gas flow and well performance are highly dependent on the properties of hydraulic fractures in shale, and hydraulic fractures may communicate with local natural fractures to form a complicated fracture network, further amplifying reservoir heterogeneity [10][11][12].The accurate modeling of a fracture network is another essential component for the forecast of shale gas production.
As conventional methods fall short in the simulation of shale gas production, numerous new models have been proposed for the macroscopic simulation of gas production in shale.These models can be roughly divided into two classes, which are numerical simulation and analytical/semi-analytical models.Both numerical simulation and analytical/semi-analytical models have their strengths and limitations.Numerical simulation offers detailed and accurate representations of reservoir behavior, but requires more computational resources, while analytical models may have to adopt some simplifying assumptions such as regular boundary and linear flow, but most of these assumptions can be justified in unconventional shale gas reservoirs [13].
The unique characteristics of shale, including but not limited to pore space confinement, multiscale pore structure, gas sorption, and stress sensitivity, cause fluid flow behaviors to substantially deviate from Darcy's law [14,15].Regardless of numerical simulation or simplified analytical solutions, a profound understanding and correct descriptions of the aforementioned gas storage and transport mechanisms are essential for the modeling of shale gas reservoirs.Macroscopic models solve a set of field equations derived based on a single represented element volume (REV) [16], thus REV-scale modeling serves as a bridge that ensures the flow principles are correctly upscaled from the microscale to the continuum scale (as shown in Figure 1) [17].Popular REV-scale modeling includes direct numerical simulation (DNS) [18,19], which directly performs computational fluid dynamics or Lattice Boltzmann simulations [20][21][22] on a real porous structure obtained through digital rock imaging, pore network modeling (PNM) [23,24] that idealizes the real rock skeleton into a pore network for improved simulation efficiency, as well as discrete element modeling [25] that accounts for particle flow or rock grain deformation.For a shale matrix with nanoscale pores, interactions between fluid molecules and rock atoms play a significant role, thus molecular dynamics simulation (MDS) [26][27][28] and density function theory (DFT) [29,30] have emerged as powerful tools to investigate the molecular-level processes occurring in shale nanopores.Despite their rapid advancements, it is difficult to directly couple microscale simulation within macroscale models due to the explosion of computational costs and the incompatibility of numerical solution procedures.Therefore, microscopic simulations are usually combined with experiments to calibrate the traditional theories used to describe the various flow mechanisms in macroscopic governing equations [17,26].These theories usually adopt closed-form equations to express the macroscopic parameters (e.g., apparent permeability, porosity, adsorbed gas content) in the continuum-scale simulations.
The past decade has witnessed a drastic increase in the literature related to the theoretical and experimental discussion, as well as the simulation, of shale gas transport at various scales.It is necessary to summarize the key advancements and findings on the related topics.In this work, we present a critical review of the macroscopic modeling techniques of shale gas flow, as well as how complex transport mechanisms and multiscale porous systems in shale can be represented in such models.The past decade has witnessed a drastic increase in the literature related to th theoretical and experimental discussion, as well as the simulation, of shale gas transpo at various scales.It is necessary to summarize the key advancements and findings on th related topics.In this work, we present a critical review of the macroscopic modelin techniques of shale gas flow, as well as how complex transport mechanisms an multiscale porous systems in shale can be represented in such models.

Gas Adsorption Isotherms
Gas molecules in shale may attach to the pore wall due to the high adsorptio affinity of the organic pore surface.According to experimental and field evidence, th amount of adsorbed gas accounts for 20% to 85% of the total gas in place in shale [9].Th accurate characterization of gas adsorption is paramount not only for the reliab estimation of gas in place, but also for the accurate prediction of gas productivity.The are two main types of gas adsorption mechanisms, namely, physical adsorption an chemical adsorption [31].Chemical adsorption is the type of adsorption that involves chemical reaction between the surface and the adsorbate.Physical adsorption is the typ of adsorption during which the electronic structure of the atom or molecule is bare perturbed [32].
The amount of gas adsorption is known to depend on pressure and temperatur The variation of adsorption concentration with pressure change at a consta temperature is termed adsorption isotherm [33].According to the shape of the adsorptio vs. pressure curve, there are six classes of adsorption isotherm for physical adsorptio which are shown in Figure 2 [34].The differences in each curve depend on the pore siz the number of adsorption layers, and the strength of gas-solid interactions durin adsorption.For example, type I characterizes gas adsorption in microporous media, an types II and IV correspond to multilayer adsorption in macroporous and mesoporou

Gas Adsorption Isotherms
Gas molecules in shale may attach to the pore wall due to the high adsorption affinity of the organic pore surface.According to experimental and field evidence, the amount of adsorbed gas accounts for 20% to 85% of the total gas in place in shale [9].The accurate characterization of gas adsorption is paramount not only for the reliable estimation of gas in place, but also for the accurate prediction of gas productivity.There are two main types of gas adsorption mechanisms, namely, physical adsorption and chemical adsorption [31].Chemical adsorption is the type of adsorption that involves a chemical reaction between the surface and the adsorbate.Physical adsorption is the type of adsorption during which the electronic structure of the atom or molecule is barely perturbed [32].
The amount of gas adsorption is known to depend on pressure and temperature.The variation of adsorption concentration with pressure change at a constant temperature is termed adsorption isotherm [33].According to the shape of the adsorption vs. pressure curve, there are six classes of adsorption isotherm for physical adsorption, which are shown in Figure 2 [34].The differences in each curve depend on the pore size, the number of adsorption layers, and the strength of gas-solid interactions during adsorption.For example, type I characterizes gas adsorption in microporous media, and types II and IV correspond to multilayer adsorption in macroporous and mesoporous materials, respectively.Types III and V correspond to gas adsorption in microporous material and mesoporous material, both with weak gas-solid interactions.The hysteresis loop in type IV and V is a result of pore condensation, and the amount of adsorbed molecules differs between condensation and evaporation at the same pressure.Type VI is rarely seen in porous materials; it usually occurs during multilayer adsorption on nonporous surfaces.Out of the six types of adsorption curves in Figure 2, the type I curve is the most commonly observed adsorption isotherm for gas sorption in shale [31], where the adsorbed gas concentration increases fast at low pressure and reaches a plateau at high pressure.Some authors have also reported the possibility of multilayer sorption in shale [35], which can be represented by type II or type IV isotherms.
molecules differs between condensation and evaporation at the same pressure.Type VI is rarely seen in porous materials; it usually occurs during multilayer adsorption on nonporous surfaces.Out of the six types of adsorption curves in Figure 2, the type I curve is the most commonly observed adsorption isotherm for gas sorption in shale [31], where the adsorbed gas concentration increases fast at low pressure and reaches a plateau at high pressure.Some authors have also reported the possibility of multilayer sorption in shale [35], which can be represented by type II or type IV isotherms.Various theoretical models have been proposed to quantitatively describe the relationship between adsorption concentration and pressure; Table 1 shows a brief summary of some commonly used models.The Langmuir model is a theoretical adsorption isotherm that can be derived through adsorption kinetics with the assumption of monolayer adsorption and a finite number of sorption sites, thus the Langmuir model describes the type I adsorption curve.The Brunauer, Emmett and Teller (BET) model is another widely applied theory used to describe gas-solid adsorption, particularly with multilayer sorption; the theory is usually adopted to measure the specific pore surface area of porous materials.The Dubinin model, on the other hand, adopts the assumption that the completion of the adsorption is signified by the filling of the pore, which is more applicable for microporous solids.Apart from the theoretical models, there are also empirical models obtained through fitting the experimental data, which include the two-parameter Freundlich model and three-parameter Toth model.The Langmuir theory, out of all the isotherms, is by far the most popular model for gas adsorption in the literature, considering its simplicity.However, as multilayer adsorption of methane in shale nanopores was observed [35], the BET isotherm is viewed as another plausible alternative for the modeling of shale gas adsorption.
Table 1.Commonly used gas adsorption models in the literature (the explanation of all variables in the equations can be found in the Nomenclature section).

Model Equation Feature
Langmuir [36]  (1) Monolayer adsorption Various theoretical models have been proposed to quantitatively describe the relationship between adsorption concentration and pressure; Table 1 shows a brief summary of some commonly used models.The Langmuir model is a theoretical adsorption isotherm that can be derived through adsorption kinetics with the assumption of monolayer adsorption and a finite number of sorption sites, thus the Langmuir model describes the type I adsorption curve.The Brunauer, Emmett and Teller (BET) model is another widely applied theory used to describe gas-solid adsorption, particularly with multilayer sorption; the theory is usually adopted to measure the specific pore surface area of porous materials.The Dubinin model, on the other hand, adopts the assumption that the completion of the adsorption is signified by the filling of the pore, which is more applicable for microporous solids.Apart from the theoretical models, there are also empirical models obtained through fitting the experimental data, which include the two-parameter Freundlich model and threeparameter Toth model.The Langmuir theory, out of all the isotherms, is by far the most popular model for gas adsorption in the literature, considering its simplicity.However, as multilayer adsorption of methane in shale nanopores was observed [35], the BET isotherm is viewed as another plausible alternative for the modeling of shale gas adsorption.
Table 1.Commonly used gas adsorption models in the literature (the explanation of all variables in the equations can be found in the Nomenclature section).

Model Equation Feature
Langmuir [ Semi-empirical based on Langmuir model

Gas Transport in Shale Nanopores
Advection and diffusion are both important transport processes for fluid in any medium.For conventional gas reservoirs with relatively high rock permeability, advection driven by pressure gradient dominates the overall gas transport.However, diffusive flux caused by the chemical potential gradient becomes significant in the shale matrix due to the extremely low permeability [41,42].This diffusive flux may be a joint effect of bulk diffusion due to the gradient of bulk phase concentration, and Knudsen diffusion due to the collision between the gas molecule and pore wall, as well as the surface diffusion of the adsorbed molecules (as shown in Figure 3) [43,44].Fick's law can be applied to calculate mass flux due to the above diffusive mechanisms.Meanwhile, the slippage effect (or Klingenberg effect) of gas flow along the pore wall also becomes notable when the pore size becomes comparable to the collision diameter of the gas molecules.As a result, the measured mass flux in the shale matrix is substantially higher than the value predicted by Darcy's law for the same pressure gradient.The magnitude of this flow enhancement depends on the Knudsen number kn, which is defined in Equation ( 6), and a large Knudsen number indicates a larger extent of flow enhancement compared to Darcy's law [45].The Knudsen number can also be used to identify flow regimes, which sequentially go through continuum flow, slippage flow, a transition regime, and free molecular flow with the increase in Knudsen number.Within different flow regimes, the dominating transport mechanism of the overall mass flux differs, and it is reported that gas transport in shale mainly falls under a slippage flow or transitional regime, wherein all the above-mentioned transport mechanisms may be influential [46].
where λ is the mean free path of the gas molecules, r p is the average radius of the matrix pore, k B is the Boltzmann constant, T is temperature in Kelvin, δ is the collision diameter of gas molecules, and p is gas pressure.
Processes 2023, 10, x FOR PEER REVIEW 6 of 22 In order to characterize this flow enhancement and simultaneously model the above processes, a convenient method is to define a lumped permeability coefficient to combine the advective flux with diffusive flux and slippage enhancements, which would result in apparent permeability as a function of gas pressure.Civan et al. [47] proposed an empirical formula for the apparent permeability during gas transport in nanotubes, which achieved a unified description of all four flow regimes that may take place during shale gas transport.Different from the empirical model by Civan, Javadpour et al. [48] modeled the mass fluxes due to slippage, Knudsen diffusion, and viscous flow separately, and the overall gas flux equals the superposition of all the considered mechanisms.Wu et al. [46,49] improved the direct summation of different transport fluxes in Javadpour's model as a weighted summation and proposed the kn-dependent weight calculation method.Based on a similar workflow, numerous authors further incorporated complex pore geometry [50,51], surface diffusion [52,53], natural fracture contribution [54], and water distribution [55,56] into the apparent permeability model by In order to characterize this flow enhancement and simultaneously model the above processes, a convenient method is to define a lumped permeability coefficient to combine the advective flux with diffusive flux and slippage enhancements, which would result in apparent permeability as a function of gas pressure.Civan et al. [47] proposed an empirical formula for the apparent permeability during gas transport in nanotubes, which achieved a unified description of all four flow regimes that may take place during shale gas transport.Different from the empirical model by Civan, Javadpour et al. weighted summation and proposed the kn-dependent weight calculation method.Based on a similar workflow, numerous authors further incorporated complex pore geometry [50,51], surface diffusion [52,53], natural fracture contribution [54], and water distribution [55,56] into the apparent permeability model by improving the flux calculation of each transport mechanism.Moreover, factors such as adsorption layer thickness [53,57], effective stress variation [58,59], and desorption-induced shrinkage [44,59] are also honored by relating the average pore radius in Equation ( 6) with pressure or stress.

Stress Sensitivity and Rock Deformation Effects
The pressure drawdown induced by shale gas production would cause significant variations in effective stress, which may result in the closure of fractures and impairments in reservoir permeability.The variation of effective stress may cause a 30% decline in well productivity for reservoirs with a low bulk modulus [60,61].There are two main approaches used to account for the effect of rock deformation induced by effective stress.The first and more rigorous approach is to solve the coupled governing equations of fluid flow and solid deformation.The simultaneous solution of coupled partial differential equations usually involves finite element discretization and involves high computational costs.Thus, a cheaper alternative to treat stress sensitivity and rock deformation is to relate permeability with pressure, and only solve the fluid transport equation.Several empirical models have been proposed to correlate permeability with pore pressure, which include linear function [62], exponential function [63], power-law function [64], as well as more complex ones [65][66][67].Out of these, the most popular one is the exponential permeability model shown in Equation (8).
where k is the current rock permeability, k i is the initial or reference permeability, p is pore pressure, p i is the initial or reference pressure, and γ is the permeability modulus.At a constant level of confining stress, the γ value can be obtained by fitting the curve of different permeability values measured at different pore pressure levels.
Apart from rock deformation due to the change of effective stress, another mechanism that may deform the shale rock is adsorption-induced dilation or desorption-induced shrinkage.For sorptive rocks such as shale or coal, gas sorption may cause substantial rock strain [68,69], which may, in turn, alter the sorption capacity of the rock.Numerous constitutive relations have been proposed for the described coupling between deformation and gas adsorption.Empirical models assume a linear relationship between matrix strain and the amount of gas adsorption [70][71][72], thus the majority of empirical models follow a similar form to the Langmuir model.Another type of constitutive correlation is given by the thermodynamics-based models; such models combine gas adsorption thermodynamics with the theory of poromechanics to incorporate the influence of gas sorption into the coupling between pore fluid and solid phase deformation [73,74].The thermodynamics-based constitutive models were built for a variety of different rock types, including microporous rock [75], mesoporous rock [74], as well as dual-porosity rocks [76].
The two mechanisms mentioned above contribute to the overall permeability variation in contradicting manners: the reduction in effective stress during pressure buildup will dilate the pore space and open the fractures, while the adsorption-induced deformation will compact the pores and fractures.As shown in Figure 4, the permeability variation due to sorption strain changes abruptly at low pressures, but less significantly at high pressures, which is similar to the variation in adsorption amount.Therefore, the overall permeability variation due to the joint effect of gas adsorption and effective stress would be a net reduction followed by a net increase.Further, these rock deformation mechanisms can be coupled with the various gas slippage and diffusion mechanisms in the tight matrix [59], resulting in an intricate evolution of the shale's apparent permeability during pressure disturbances.
variation due to sorption strain changes abruptly at low pressures, but less significantly at high pressures, which is similar to the variation in adsorption amount.Therefore, the overall permeability variation due to the joint effect of gas adsorption and effective stress would be a net reduction followed by a net increase.Further, these rock deformation mechanisms can be coupled with the various gas slippage and diffusion mechanisms in the tight matrix [59], resulting in an intricate evolution of the shale's apparent permeability during pressure disturbances.

Representation of Multiscale Gas Transport
The complexity of the inherent porous structure and the induced hydraulic fractures contributes to the multiscale and multi-porosity feature of shale gas reservoirs [77].

Representation of Multiscale Gas Transport
The complexity of the inherent porous structure and the induced hydraulic fractures contributes to the multiscale and multi-porosity feature of shale gas reservoirs [77].Massive efforts have been made to characterize multiscale gas transport in shale at the continuum scale.The traditional theory uses the single-porosity model to homogenize the reservoir rock as an equivalent medium.However, this simplistic treatment fails to accurately capture the complexity of shale reservoirs, particularly the dynamic flow between matrix pores and fractures.To address these limitations, multi-porosity models have been developed, such as the dual-porosity model, the triple-porosity model, the quadruple-porosity model, etc.These models leverage the distinct characteristics of fractures or the matrix to provide a more comprehensive understanding of the multiscale phenomena associated with shale gas transport.
For the shale gas reservoir, the single-porosity model simplifies diverse pore structure types, including organic matrix, inorganic micropores, and microfractures, as a single continuous medium.The gas adsorption amount throughout the entire space can be determined using the equivalent pore pressure.In order to simplify the complexity of forward simulation and improve computational efficiency, the single-porosity model is commonly used to estimate reservoir parameters and account for multiphase flow.This model has been proven useful in various applications, such as analyzing the flowback data [78,79], extracting fracture information [10,80], and simulating the early production period [81][82][83] in shale gas reservoirs.However, the single-porosity model is too simple to adequately characterize the multi-porosity and multiscale transport phenomena observed in shale reservoirs.
The pioneering work of Warren and Root [84] introduced the dual-porosity model to simulate fractured oil/gas reservoirs.Subsequently, efforts were made to extend this model to investigate shale gas reservoirs [85,86].In this context, the porous media, including inorganic pores and organic pores, are incorporated into a matrix system, while the fractured media, comprising natural fractures and secondary fractures, are integrated into a fracture system.Two mechanisms, namely, pseudo-steady-state crossflow and unsteady-state crossflow, have been widely used to describe the mass transfer between the matrix and fracture.The pseudo-steady-state inter-porosity transfer model is a relatively simplified method of characterizing the fluid transport between the matrix and fractures.It assumes that the fluid supply rate from the matrix to the fractures is directly proportional to the pressure difference between them.This model is particularly suitable for reservoirs involving small matrix particles or high matrix permeability.The unsteady-state inter-porosity transfer model further considers the transient seepage process of the fluid within the matrix particles.It can be categorized into three types based on the shape of the matrix particles, including cubic, cylindrical, and spherical.The unsteady model is particularly suitable for the simulation of shale gas transport, considering the extremely low permeability of shale reservoirs.
In addition, the unsteady dual-porosity model has further incorporated various complex seepage mechanisms within the matrix, such as adsorption, desorption, slippage effects, and Knudsen diffusion [87][88][89].However, the introduction of these mechanisms would significantly increase the nonlinearity of the mathematical model.Parameters such as adsorption amount and formation permeability are all functions of pore pressure, posing challenges for the solution of the mathematical model.Existing studies commonly calculate the values of the adsorption compression factor and apparent permeability at the initial pressure, and regard them as constants when solving the model.Consequently, the effects of adsorption and slippage diffusion on the compression factor and apparent permeability are not rigorously accounted for.
As shown in Figure 5a,b, in order to further improve the accuracy of shale gas seepage simulation, the matrix and fracture system of the dual-porosity mode have been further divided into triple-and quadruple-porosity models.These advanced models offer a more detailed characterization of the complex pore structure and flow behavior within shale reservoirs.Ci-qun [90,91] developed the triple-porosity model to investigate the radial flow of micro-compressible fluids under pseudo-steady-state flow conditions.In addition, a triple-porosity model composed of two fracture systems and one matrix system has been proposed [92,93], which less closely considered the effects of desorption, diffusion, slip flow, and the geometric characteristics of fractures.Al-Ahmadi and Wattenbarger [94] also established a triple-porosity model consisting of two fracture systems and a matrix system; the model honors the linear crossflow between the two fracture media and the gas desorption mechanism within matrix pores.In order to account for the flow mechanism ranging from the nanoscale (Kundsen diffusion and slip flow) to the macroscale (adsorption, desorption, and Darcy flow), another type of triple-porosity model was developed for shale gas reservoirs [95,96].These models incorporate three distinct porous media: kerogen, matrix, and fractures.Furthermore, a quadruple-porosity media model was developed for shale gas reservoirs [97,98].The multiple interacting continua (MINC) concept [99] is usually adopted to treat the interaction among different continua when multiple porosity types are included.This model incorporates four systems: organic matter, kerogen pore, inorganic matrix, and natural fractures.In their work, the porous kerogen is assumed to comprise kerogen pores and organic matters; kerogen pores may contain free gas or adsorbed gas, while organic matters may contain dissolved gas.Gas first dissolves from the solid organic matter into the kerogen pores, and then diffuses into the inorganic matrix, and finally into the natural fractures.
While introducing additional types of porous media can potentially enhance the accuracy of shale gas seepage simulation, it is important to note that indiscriminately increasing the number of porosity types can reduce the practicality of the model.The inclusion of more pore types substantially increases the computational cost and complexity of solutions.Furthermore, the multi-porosity model involves a large number of input parameters, which can be challenging to accurately obtain.The efficiency and accuracy of input parameters greatly influence the practicality and reliability of the model.Considering both the complexity and the practicability of the model, the unsteady-state dual-porosity model is by far the most adopted approach for shale gas transport in the literature and major simulation codes.This model can accurately characterize the multiscale seepage characteristics of shale reservoirs while ensuring the efficient solution of the model.interaction among different continua when multiple porosity types are included.This model incorporates four systems: organic matter, kerogen pore, inorganic matrix, and natural fractures.In their work, the porous kerogen is assumed to comprise kerogen pores and organic matters; kerogen pores may contain free gas or adsorbed gas, while organic matters may contain dissolved gas.Gas first dissolves from the solid organic matter into the kerogen pores, and then diffuses into the inorganic matrix, and finally into the natural fractures.While introducing additional types of porous media can potentially enhance the accuracy of shale gas seepage simulation, it is important to note that indiscriminately increasing the number of porosity types can reduce the practicality of the model.The inclusion of more pore types substantially increases the computational cost and complexity of solutions.Furthermore, the multi-porosity model involves a large number of input parameters, which can be challenging to accurately obtain.The efficiency and accuracy of input parameters greatly influence the practicality and reliability of the model.Considering both the complexity and the practicability of the model, the unsteady-state dual-porosity model is by far the most adopted approach for shale gas transport in the literature and major simulation codes.This model can accurately characterize the multiscale seepage characteristics of shale reservoirs while ensuring the efficient solution of the model.

Macroscopic Modeling and Solution Techniques of Shale Gas Flow
The theoretical models concluded in Section 2 characterize the various gas transport and storage mechanisms by relating the continuum-scale flow parameters (e.g., porosity, permeability, adsorption rate) with gas concentration or pressure, and inserting these theoretical models into the governing equations of each porosity system constitutes the final mathematical model describing multiscale gas transport in fractured shale gas reservoirs.Therefore, the continuum-scale simulation of shale gas production involves the solution of nonlinear partial differential equations due to the dependence of model parameters on gas pressure.Two distinct approaches, analytical and numerical methods, may be adopted to solve the obtained mathematical model and determine the evolution of gas pressure or production rate.This section introduces the overall procedure of each approach, and summarizes the advances of each approach over recent years.

Numerical Simulation
Numerical simulation methods involve discretizing the reservoir into a grid system so that the original governing differential equations for gas flow can be transformed into a system of algebraic equations, which will then be solved numerically at each grid block and time step to simulate the gas flow behavior.
The major numerical simulation techniques for shale gas flow include the finite difference method (FDM), finite volume method (FVM), and finite element method (FEM).The finite difference method (FDM) involves discretizing the continuous domain into a set of discrete points and approximating derivatives using finite difference approximations on Cartesian grids [16].Unlike the FDM, the FVM considers the conservation laws integrated over each control volume.The method accounts for the fluxes across the boundaries of these volumes and conserves the quantities within each controlled volume [100].The FEM solves PDEs by dividing the domain into a collection of smaller subdomains called finite elements.It approximates the solution within each element using a set of basis functions [101].For the simulation of shale gas production, the FDM is the most widely adopted approach considering its high computational efficiency; the FVM may be used when the accurate representation of a complex fracture or boundary is required.For stress-sensitive shale gas reservoirs, the simultaneous solution of pore pressure and rock stress usually requires the FEM so as to achieve the accurate solution of mechanics equations [101][102][103].
As mentioned in the previous sections, the production performance of shale gas is highly dependent on the properties and the distribution of fractures, thus the accurate representation of the fracture network is a key issue for the numerical simulation of shale gas production.Early studies [104][105][106] usually assumed that, except for the main fractures created by hydraulic fracturing, secondary fractures can activate and connect with the natural fractures to form a continuous fracture network.Therefore, the dual-porosity model [104][105][106] has been the favored method to represent the fracture network so that all fractures are collectively treated by an equivalent continuum.The dual-porosity model allows for simulating fractured reservoirs using a simple orthogonal grid system, and is computationally efficient.However, the dual-porosity model presumes the existence of an REV containing both matrix and fractures, thus the method may be inaccurate when the fracture density is small, as the geometry of each fracture cannot be treated explicitly.With the advancements in micro-seismic inversion technology, the number and distribution of fractures may be determined.To accurately represent the characteristics of each fracture, the traditional discrete fracture network (DFN) model [107][108][109][110] is adopted to represent a fractured rock mass as a network of discrete fractures or fracture sets.These fractures can have different orientations, lengths, and hydraulic or mechanical properties.The model takes into account the geometrical characteristics and connectivity of fractures within the rock mass, providing a detailed representation of the complex fracture network.However, such detailed modeling of fracture characteristics can only be accomplished on unstructured grid systems, and FVM or FEM must be adopted to solve the governing equations.Meanwhile, the permeability of the fracture nodes/grids is usually several magnitudes higher than the neighboring matrix nodes/grids, which may create severe convergence issues during the solution.Therefore, the high accuracy of the DFN method comes with the expense of the high computational cost.
To balance simulation accuracy and efficiency, Moinfar et al. [111] extended the idea proposed by Lee et al. [112], and developed a new method named embedded discrete fracture model (EDFM) to model the hydraulic fractures.An illustration of EDFM and the comparison with other methods is depicted in Figure 6.The EDFM treats the fractures as discrete fractures embedded in the matrix grids, and the mass transfer between matrix and fractures is calculated in a way similar to the dual-porosity approach, while the mass flux from the matrix grid to the local fracture depends on how the fracture penetrates the local grid.This allows for the accurate representation of the fracture orientation, length, and aperture on the Cartesian system without local grid refinement, which significantly reduces the computational cost of the numerical model.However, the EDFM approach is reported to be inaccurate when the fracture permeability is close to the matrix permeability, due to the approximation of the matrix-fracture mass flux.Tene et al. [113] proposed a projectionbased embedded discrete fracture model (pEDFM) by revising the transmissibility between matrix and fractures; such an improvement ensures that pEDFM applies to any permeability contrast between the two media, and it was also shown that pEDFM can effectively reduce the error of EDFM during multiphase flow conditions [114].The EDFM method has gained great popularity over recent years due to its high accuracy and efficiency in fracture modeling; plenty of efforts have been made to improve the EDFM method, which include the incorporation of geomechanics and fracture aperture changes [115], the combination of EDFM with dual-porosity models to account for fractures with different scales [116], as well as the application of EDFM in compositional simulations [117,118].
Wang and Fidelibus [119] gave an example of the application of EDFM to simulate the production of a multi-fractured horizontal well in Barnett shale.The input reservoir and well parameters are listed in Table 2.The numerical simulation codes were adopted from an open-source EDFM simulator developed by Wang and Fidelibus [119]; the simulator incorporates gas desorption, slippage, diffusion, as well as stress-sensitive fracture permeability.Besides the above mechanisms, the numerical simulator can also account for the complex geometry of the hydraulic/natural fractures, demonstrating a strong predictive capability even with various reservoir complexities.Figure 7 shows the pressure distribution after 1000 days of well production, as well as the history-matching results.
pEDFM applies to any permeability contrast between the two media, and it was also shown that pEDFM can effectively reduce the error of EDFM during multiphase flow conditions [114].The EDFM method has gained great popularity over recent years due to its high accuracy and efficiency in fracture modeling; plenty of efforts have been made to improve the EDFM method, which include the incorporation of geomechanics and fracture aperture changes [115], the combination of EDFM with dual-porosity models to account for fractures with different scales [116], as well as the application of EDFM in compositional simulations [117,118].Wang and Fidelibus [119] gave an example of the application of EDFM to simulate the production of a multi-fractured horizontal well in Barnett shale.The input reservoir and well parameters are listed in Table 2.The numerical simulation codes were adopted from an open-source EDFM simulator developed by Wang and Fidelibus [119]; the simulator incorporates gas desorption, slippage, diffusion, as well as stress-sensitive fracture permeability.Besides the above mechanisms, the numerical simulator can also account for the complex geometry of the hydraulic/natural fractures, demonstrating a  strong predictive capability even with various reservoir complexities.Figure 7 shows the pressure distribution after 1000 days of well production, as well as the history-matching results.

Analytical and Semi-Analytical Methods
Analytical or semi-analytical models aim to provide simplified, closed-form solutions to describe the flow behavior in shale gas reservoirs.These models are based on mathematical equations and assumptions that simplify the complex physics involved.While they may not capture all the intricacies of the reservoir, analytical and

Analytical and Semi-Analytical Methods
Analytical or semi-analytical models aim to provide simplified, closed-form solutions to describe the flow behavior in shale gas reservoirs.These models are based on mathematical equations and assumptions that simplify the complex physics involved.While they may not capture all the intricacies of the reservoir, analytical and semi-analytical models are often used as screening tools, providing quick estimates and insights into key parameters and production behavior.They are computationally efficient and can help in initial assessments and quick sensitivity analyses.Thus, analytical/semi-analytical processes have gained great popularity in production forecasting, reservoir characterization, and production optimization applications due to their high computational efficiency.
The analytical solution to the partial differential equation (PDE) governing gas flow typically involves the following steps: The first step is to convert the governing PDE into ordinary differential equations (ODE) as the direct solution of PDE is usually difficult.For this purpose, the Laplace transform is frequently adopted to convert the PDE in the time domain into an ODE in Laplace space, and the Fourier transform is usually employed to eliminate the space variables and is especially useful when vertical flow occurs.Besides this, several works have applied the Boltzmann transformation [120][121][122] to model fluid flow in tight formations during infinite-acting transient flow, but this technique imposes special requirements on the boundary conditions of the investigated problem.The second step is to solve the ODE and obtain the closed-form expression of pressure distribution and flow rate evolution, and perform the inverse transform to attain the solution in the time domain.However, not all problems are analytically solvable despite the simplifying assumptions, especially for the modeling of shale gas transport due to the nonlinear nature of gas transport.To treat these nonlinearities, the pseudo-variable approach [123,124] is predominantly adopted to remove the nonlinearity caused by the dependence of gas properties on pressure, while the perturbation technique proves applicable in treating the nonlinearities caused by stress-sensitive permeability [60,61].Other methods for solving nonlinear differential equations also include Green's function approach [83,125] and the piecewise constant diffusivity approach [81].
Based on the adopted assumptions, analytical/semi-analytical models can be divided into two groups, namely, linear-flow-based models and line-source-based models.As shown in Figure 8, linear-flow-based models assume that gas transport into the horizontal well is involves separate linear flows in different regions, and such an assumption is justified by the large permeability contrast among these regions.The line-source-based model is built by superposing the pressure variation caused by a line source along the trajectory of hydraulic fractures, and the overall pressure drawdown after superposition is the pressure drop caused by the multi-fractured horizontal well.
Processes 2023, 10, x FOR PEER REVIEW 13 of 22 forecasting, reservoir characterization, and production optimization applications due to their high computational efficiency.
The analytical solution to the partial differential equation (PDE) governing gas flow typically involves the following steps: The first step is to convert the governing PDE into ordinary differential equations (ODE) as the direct solution of PDE is usually difficult.For this purpose, the Laplace transform is frequently adopted to convert the PDE in the time domain into an ODE in Laplace space, and the Fourier transform is usually employed to eliminate the space variables and is especially useful when vertical flow occurs.Besides this, several works have applied the Boltzmann transformation [120][121][122] to model fluid flow in tight formations during infinite-acting transient flow, but this technique imposes special requirements on the boundary conditions of the investigated problem.The second step is to solve the ODE and obtain the closed-form expression of pressure distribution and flow rate evolution, and perform the inverse transform to attain the solution in the time domain.However, not all problems are analytically solvable despite the simplifying assumptions, especially for the modeling of shale gas transport due to the nonlinear nature of gas transport.To treat these nonlinearities, the pseudo-variable approach [123,124] is predominantly adopted to remove the nonlinearity caused by the dependence of gas properties on pressure, while the perturbation technique proves applicable in treating the nonlinearities caused by stress-sensitive permeability [60,61].Other methods for solving nonlinear differential equations also include Green's function approach [83,125] and the piecewise constant diffusivity approach [81].
Based on the adopted assumptions, analytical/semi-analytical models can be divided into two groups, namely, linear-flow-based models and line-source-based models.As shown in Figure 8, linear-flow-based models assume that gas transport into the horizontal well is involves separate linear flows in different regions, and such an assumption is justified by the large permeability contrast among these regions.The line-source-based model is built by superposing the pressure variation caused by a line source along the trajectory of hydraulic fractures, and the overall pressure drawdown after superposition is the pressure drop caused by the multi-fractured horizontal well.The simplest linear-flow-based model is the 1D linear flow of gas from the reservoir to the infinitely conductive hydraulic fracture.Most of the rate-transient analysis models adopt this configuration to extract formation properties [78][79][80], and this simplified geometry allows the mathematical model to incorporate more complex physics, such as multiphase flow [80], fracture extension [126,127], adsorption-induced deformation [43], and gas slippage and diffusion [128].Brown et al. [129] proposed a tri-linear flow model by dividing the hydrocarbon transport into a sequence of linear flows: linear flow in the primary fracture; linear flow in the stimulated reservoir volume (SRV) region between two adjacent primary fractures; linear flow in the non-stimulated region.To better account for the inherent heterogeneity of the fractured shale gas formation, several works further divided the SRV region into more regions with different permeability, such as the penta-linear model developed by Stalgorova and Mattar [130], and the septa-linear models proposed by Yuan et al. [131] and Zeng [132].However, the abovementioned models assume that the reservoir properties remain constant in the same region, but change abruptly across the boundaries, which is not consistent with the continuous permeability variation in the actual shale formation.To address this issue, Fan and Ettehadtavakkol [133,134], as well as Wang et al. [135,136], assumed that the fracture network in the SRV region shows a fractal geometry, and adopted power-law functions to describe the continuous variation of fracture properties.Shahamat et al. [137] applied the perturbation theory to deal with the exponential deterioration of reservoir permeability with distance in the SRV region.Zhang and Xu [138] proposed a composite linear flow model to account for the permeability heterogeneity in the main fractures and SRV region, and any arbitrary permeability functions can be treated as inputs to account for the described heterogeneities.Figure 9 shows the evolution of linear-flow-based models used in treating the heterogeneity of the fractured reservoir.
Processes 2023, 10, x FOR PEER REVIEW 14 of 22 The simplest linear-flow-based model is the 1D linear flow of gas from the reservoir to the infinitely conductive hydraulic fracture.Most of the rate-transient analysis models adopt this configuration to extract formation properties [78][79][80], and this simplified geometry allows the mathematical model to incorporate more complex physics, such as multiphase flow [80], fracture extension [126,127], adsorption-induced deformation [43], and gas slippage and diffusion [128].Brown et al. [129] proposed a tri-linear flow model by dividing the hydrocarbon transport into a sequence of linear flows: linear flow in the primary fracture; linear flow in the stimulated reservoir volume (SRV) region between two adjacent primary fractures; linear flow in the non-stimulated region.To better account for the inherent heterogeneity of the fractured shale gas formation, several works further divided the SRV region into more regions with different permeability, such as the penta-linear model developed by Stalgorova and Mattar [130], and the septa-linear models proposed by Yuan et al. [131] and Zeng [132].However, the abovementioned models assume that the reservoir properties remain constant in the same region, but change abruptly across the boundaries, which is not consistent with the continuous permeability variation in the actual shale formation.To address this issue, Fan and Ettehadtavakkol [133,134], as well as Wang et al. [135,136], assumed that the fracture network in the SRV region shows a fractal geometry, and adopted power-law functions to describe the continuous variation of fracture properties.Shahamat et al. [137] applied the perturbation theory to deal with the exponential deterioration of reservoir permeability with distance in the SRV region.Zhang and Xu [138] proposed a composite linear flow model to account for the permeability heterogeneity in the main fractures and SRV region, and any arbitrary permeability functions can be treated as inputs to account for the described heterogeneities.Figure 9 shows the evolution of linear-flow-based models used in treating the heterogeneity of the fractured reservoir.Different from the linear-flow-based models, the line-source-based models do not depend on the assumption of linear flow and can reflect more flow regimes.This is especially necessary when the hydraulic fractures are not perpendicular to the horizontal well.Early linear-flow-based models assumed that the fracture is infinitely conductive, such that the pressure variation within the reservoir can be readily obtained by superposing the line source solution [139,140].Later works managed to account for gas Different from the linear-flow-based models, the line-source-based models do not depend on the assumption of linear flow and can reflect more flow regimes.This is especially necessary when the hydraulic fractures are not perpendicular to the horizontal well.Early linear-flow-based models assumed that the fracture is infinitely conductive, such that the pressure variation within the reservoir can be readily obtained by superposing the line source solution [139,140].Later works managed to account for gas flow in main hydraulic fractures by solving a system of equations constraining the flow rate distribution along the hydraulic fractures, and the same method was later extended to model finite-conductivity fractures with irregular shapes [141,142].To account for the reservoir heterogeneity caused by SRV, the line source solutions for a composite system comprising a high-permeability inner zone and a low-permeability outer zone are proposed [143,144], but such treatment can only account for heterogeneity in the radial direction.Ren et al. [145] further proposed the line-source-based model for shale gas reservoirs with heterogeneity along the wellbore direction.Besides this, massive efforts have been made to couple the line-source-based models with multiple-porosity models [99,139,146], such that complex transport mechanisms can be considered in the multiscale porous system, e.g., gas diffusion, adsorption, and stress-sensitivity.
Here, we also demonstrate the application of an analytical model to extract formation information from a shale gas reservoir.The target reservoir is located in the southwest region of China.The five-region analytical model developed by Stalgorova and Mattar [130] is employed to perform the forward prediction of the gas production rate.Taking advantage of the speed of the analytical model, the simulated annealing algorithm is employed to automatically interpret the unknown fracture parameters that lead to the best matching between the analytical model and the actual production data.Figure 10 displays the matching between the analytical solution and the actual production data, and Table 3 summarizes all the input data of the analytical model and the interpreted information of hydraulic fractures.
Here, we also demonstrate the application of an analytical model to extract formation information from a shale gas reservoir.The target reservoir is located in the southwest region of China.The five-region analytical model developed by Stalgorova and Mattar [130] is employed to perform the forward prediction of the gas production rate.Taking advantage of the speed of the analytical model, the simulated annealing algorithm is employed to automatically interpret the unknown fracture parameters that lead to the best matching between the analytical model and the actual production data.Figure 10 displays the matching between the analytical solution and the actual production data, and Table 3 summarizes all the input data of the analytical model and the interpreted information of hydraulic fractures.

Conclusions
This review covers the fundamental flow and storage mechanisms, characterization methods for multiscale gas transport, as well as macroscopic modeling approaches for shale gas transport.The careful treatment of these three aspects is required for the reliable prediction of shale gas production.
Several storage and flow mechanisms are found to be influential for shale gas production, which mainly include gas adsorption, gas slippage and diffusion, as well as rock deformation.Although various theoretical models have been proposed to describe the gas adsorption process, the classical Langmuir model remains the most popular theory used to model shale gas adsorption.Gas slippage, gas diffusion, and rock deformation are completely different physical processes from gas advection (Darcy flow); the accurate modeling of these processes requires the simultaneous solution of the gas diffusion equation and solid deformation equation, which requires considerably high computational resources.Therefore, the current approach is to lump the effects of all these mechanisms into an apparent permeability term.
The major efforts made to account for the multiscale transport of shale gas have been undertaken through the multi-porosity model.However, separating the formation rock into systems with greater porosity may result in impractical modeling and redundant input parameters.By far the most popular theory for the macroscopic modeling of shale gas transport is the unsteady-state dual-porosity model.The macroscopic modeling approach can also be categorized into two types based on the solution strategy, which are numerical simulations and analytical/semi-analytical models.FDM, FVM, and FEM are common solution methods used to discretize the reservoir and numerically solve the governing PDEs, and the EDFM method has become the most popular approach to modeling the fractures via numerical simulation due to its balanced simulation accuracy and efficiency.Analytical/semi-analytical models, on the other hand, directly solve the governing equation without discretizing the domain.The linear-flow-based models are usually adopted to obtain analytical solutions for strongly heterogeneous systems, while the line-source-based models are more suitable for cases with irregular fracture geometries.
With the rapid development of computational capabilities and parallel computation algorithms, the overall trend of macroscopic modeling for shale gas reservoirs will be towards fine-mesh-based multi-physics numerical simulation with precise fracture representation.All the coupled processes, such as gas diffusion and rock deformation, may be rigorously solved together with the gas flow equation.However, this does not necessarily mean that analytical/semi-analytical models would be completely eliminated, as the high computational efficiency of such models would still make them useful in inverse problems, i.e., the rate-transient analysis of shale gas reservoirs.

Figure 1 .
Figure 1.Various modeling approaches for shale gas transport at different scales.

Figure 1 .
Figure 1.Various modeling approaches for shale gas transport at different scales.

Figure 2 .
Figure 2. Six types of adsorption isotherm curves (with each curve explained in the text).

Figure 2 .
Figure 2. Six types of adsorption isotherm curves (with each curve explained in the text).

Figure 3 .
Figure 3. Illustration of different gas transport mechanisms in shale matrix.

Figure 3 .
Figure 3. Illustration of different gas transport mechanisms in shale matrix.
[48]  modeled the mass fluxes due to slippage, Knudsen diffusion, and viscous flow separately, and the overall gas flux equals the superposition of all the considered mechanisms.Wu et al.[46,49] improved the direct summation of different transport fluxes in Javadpour's model as a Processes 2023, 11, 2766 6 of 21

Figure 4 .
Figure 4. Possible evolution pattern of shale permeability under rock deformation and slippage effects.

Figure 4 .
Figure 4. Possible evolution pattern of shale permeability under rock deformation and slippage effects.

Figure 5 .
Figure 5. Illustration of different multi-porosity models in the literature.

Figure 5 .
Figure 5. Illustration of different multi-porosity models in the literature.

Figure 7 .
Figure 7. History-matching results and pressure distribution at 1000 days.

Figure 7 .
Figure 7. History-matching results and pressure distribution at 1000 days.

Figure 8 .
Figure 8.The two types of analytical/semi-analytical models for the modeling of a multi-fractured horizontal well.

Figure 8 .
Figure 8.The two types of analytical/semi-analytical models for the modeling of a multi-fractured horizontal well.

Figure 9 .
Figure 9. Evolution of composite linear flow models.

Figure 9 .
Figure 9. Evolution of composite linear flow models.

Figure 10 .
Figure 10.Matching between the analytical solution and actual production data.Figure 10.Matching between the analytical solution and actual production data.

Figure 10 .
Figure 10.Matching between the analytical solution and actual production data.Figure 10.Matching between the analytical solution and actual production data.

Table 2 .
The input reservoir and well parameters.

Table 2 .
The input reservoir and well parameters.

Table 3 .
The input parameters and interpretation results.

Table 3 .
The input parameters and interpretation results.
maximum adsorption concentration in the Langmuir model, kg/kg.n m equivalent mono-layer adsorption concentration in BET model, kg/kg.n D maximum adsorption concentration in Dubinin adsorption model, kg/kg.n F empirical distribution coefficient in Freundlich adsorption model, dimensionless.p pore pressure, Pa.p 0 vapor pressure of the adsorptive bulk liquid phase in BET model, Pa.p L pressure at which the adsorption concentration reaches half of the maximum concentration in Langmuir model, Pa.p S saturation pressure of the adsorbent in Dubinin adsorption model, Pa.r p average radius of the matrix pore, m.T temperature, K. x empirical power index in Freundlich adsorption model, dimensionless.
λ mean free path of the gas molecules, m. δ collision diameter of gas molecules, m. γ permeability modulus, Pa.