Shale gas reservoir modeling and production evaluation considering complex gas transport mechanisms and dispersed distribution of kerogen

Stimulated shale reservoirs consist of kerogen, inorganic matter, secondary and hydraulic fractures. The dispersed distribution of kerogen within matrices and complex gas flow mechanisms make production evaluation challenging. Here we establish an analytical method that addresses kerogen-inorganic matter gas transfer, dispersed kerogen distribution, and complex gas flow mechanisms to facilitate evaluating gas production. The matrix element is defined as a kerogen core with an exterior inorganic sphere. Unlike most previous models, we merely use boundary conditions to describe kerogen-inorganic matter gas transfer without the instantaneous kerogen gas source term. It is closer to real inter-porosity flow conditions between kerogen and inorganic matter. Knudsen diffusion, surface diffusion, adsorption/desorption, and slip corrected flow are involved in matrix gas flow. Matrix-fracture coupling is realized by using a seven-region linear flow model. The model is verified against a published model and field data. Results reveal that inorganic matrices serve as a major gas source especially at early times. Kerogen provides limited contributions to production even under a pseudo-steady state. Kerogen properties’ influence starts from the late matrix-fracture inter-porosity flow regime, while inorganic matter properties control almost all flow regimes except the early-mid time fracture linear flow regime. The contribution of different linear flow regions is also documented.


Introduction
Recent years have seen a dramatic increase in shale gas production as shale reservoir development using multistage fractured horizontal wells (MFHWs) has achieved great success and is one of the main focuses in the petroleum industry (Li et al. 2017;Wang et al. 2017). A stimulated shale reservoir can be classified into four sub-systems, including kerogen, inorganic matrix, secondary and hydraulic fracture systems. Due to different spatial structures of these sub-systems, their properties control well performance at different timescales. Generally, gas production first comes from fractures. At that time, the matrix gas inflow is negligible (Wasaki and Akkutlu 2015). In a later stage that occupies a long production period, the gas stored in the matrix contributes to production. As a result, the overall production curve will first experience a sharp decline which is followed by a long tail (Olorode et al. 2017). Although we know well about this production behavior, there are still many challenges in accurately modeling and evaluating the well performance in shale plays due to their unique characteristics. Shale rocks consist of pore networks with a wide range of sizes and different types of minerals. These pore networks lead to complex gas transport mechanisms, involving the slip flow, Knudsen diffusion, surface diffusion, and desorption at multiple scales (Olorode et al. 2017;Wang and Reed 2009;Akkutlu and Fathi 2012). Apart from that, the SEM images of shale samples also reveal the discontinuous distribution of kerogen within the inorganic matrix (Wasaki and Akkutlu Edited by Yan-Hua Sun 2015; Ambrose et al. 2012;Ruppel and Loucks 2008;Yang et al. 2019;Liu et al. 2019). The kerogen presents as finely dispersed phases within the inorganic matrix (Olorode et al. 2017). Figure 1 shows the typical distribution of kerogen in SEM images of shale samples. According to Olorode et al. (2017), the element sizes in reservoir simulation can easily be five orders of magnitude larger than the SEM image size. Gas movement from the kerogen surface to the lowgas-concentration areas of the inorganic matrix is not an instantaneous process but a gradual procedure depending on the inorganic matrix transport properties. Moreover, after fracturing treatments, hydraulic fractures are induced and connected with existing natural fractures ). The macro-hydraulic fractures and self-supported natural (secondary) fractures serve as the main flow channels of gas (Tian 2014;Li et al. 2018). These characteristics provide unique challenges for shale gas extraction simulation.
In the literature, a number of analytically or semi-analytically based models using fast and gridless approaches have been proposed for tight formations with MFHWs because the refined numerical models are usually time-consuming, inconvenient in the iterative calculation, and require extensive input data. Most importantly, analytical models can be employed to validate numerical models. Brown et al. (2011) modified the trilinear flow solution with an inner dual-porosity region to simulate MFHWs. Based on Brown's trilinear flow model, Stalgorova and Mattar (2013) established a five-region linear flow model where the region between adjacent hydraulic fractures includes a stimulated reservoir volume and an unstimulated reservoir volume. Recently, seven-region and ten-region linear flow models were published for shale and tight sand reservoirs (Zeng et al. 2017(Zeng et al. , 2019a. With additional flow regions, composite heterogeneous reservoirs that are partially penetrated by fractures in vertical and horizontal directions and fracture damage can be addressed. By characterizing unconventional reservoirs with the fractal theory, alternative linear flow idealization formulas were proposed (Ozcan et al. 2014;Wang et al. 2015Wang et al. , 2018Zeng et al. 2020). Although they are applicable for interpreting pressure and rate responses, their reliability is heavily dependent on the fractal properties they choose. To model the fluid transfer from the shale matrix to fracture networks, Ozkan et al. (2010) followed the approaches of Ertekin et al. (1986) and developed a spherical matrix element model that includes the direct superposition of Darcy's velocity and slip corrected velocity in the matrix. By utilizing similar methods, Apaydin et al. (2012) examined the effects of microfractures on matrix effective permeability. Although they used spherical matrix blocks, the fluid flow in the outer-layer microfractures is simplified as linear flow. This assumption is acceptable when the microfractures only exist in a very thin layer outside the matrix core and is not suitable for the interspersed kerogen core with a relatively thick exterior inorganic matrix layer. Javadpour and Ghanbarnezhad Moghanloo (2013) used a concentric-sphere model to simulate the gas diffusion process from kerogen surfaces to inorganic matter and calculate inorganic matter effective diffusion coefficients. However, they ignored the gas diffusion within kerogen and assumed that the gas concentration values on kerogen core and inorganic matter surfaces are constants, which is not applicable for the actual gas production process. Semi-analytical models that involve Green's function and the source/sink method were also used to test and evaluate MFHWs' performance with desorption effects or fracture stress sensitivity (Yao et al. 2013(Yao et al. , 2016. Based on the superposition theory, Zhang et al. (2018) developed a semi-analytical method considering the irregular fracture geometry and complex fracture patterns of fractured vertical wells. Zhao et al. (2013) derived a triple-porosity semi-analytical model for both free gas and adsorbed gas transport. However, gas transport mechanisms of the above spherical matrix and semi-analytical models are too simple compared with the actual gas flow process. In terms of multiple gas transport mechanisms, Huang et al. (2015) proposed a triple-porosity model for unstimulated vertical wells considering the inter-porosity flow between kerogen and inorganic clay. However, they merely used a kerogen source term in the inorganic matrix flow equation indicating that the kerogen gas is instantaneously and uniformly distributed in the surrounding inorganic matrix. And the unstimulated vertical well is not suitable for shale reservoir development. Some other analytical and semi-analytical quadruple-porosity models Guo et al. 2015;He et al. 2016;Fan and Ettehadtavakkol 2017;Zeng et al. 2019b) seem to be versatile enough to deal with the predominant flow behavior of fractured shale formations. However, they also utilized the source term method or slab matrix blocks. Using slab blocks means the cross-sectional areas for kerogen and inorganic matrix flow are the same, which will overestimate kerogen's contribution. Besides, Kucuk and Sawyer (1980) compared different element shapes when analyzing the Devonian gas shale data. They found that spherical and cylindrical elements perform better than Warren-Root and Kazemi cubic models. And some gas transport mechanisms and real-gas effects are neglected or are not properly described. After a thorough literature review, we find that typical reservoir models based on single-porosity, dual-porosity, triple-porosity, and even quadruple-porosity assumptions only model limited characteristics of shale reservoirs. Table 1 further summarizes the features of typical models that handle kerogen, inorganic matrix, and fracture inter-porosity flow within shale formations. The biggest advantage of this model is that it considers the dispersed distribution of kerogen within inorganic matrices and uses the flux and pressure continuity boundary conditions to characterize the gradual gas transfer from kerogen to inorganic matrices instead of using a source term. Besides, complex gas transport mechanisms are included as well.

Matrix model
The fractured shale reservoir is divided into kerogen, inorganic matrices, natural (secondary) fractures, and hydraulic fractures, as shown in Fig. 2. The matrix model is extended from de Swaan's sphere-element model (de Swaan 1976), as shown in Fig. 2a. We utilize two concentric spheres to represent the kerogen core and the surrounding inorganic matrix respectively. To derive the matrix model, the following assumptions are made: (1) The finely dispersed kerogen is represented by the interior sphere, while the surrounding inorganic matrix is depicted by the exterior sphere with different properties and gas transport mechanisms. The size of the inorganic sphere, which can be interpreted from field data, is selected from the literature (Apaydin et al. 2012). In reservoir simulation, the spherical matrix element size can easily be up to five orders of magnitude larger than the SEM images (Olorode et al. 2017). After determining the exterior sphere size, the kerogen core size is calculated by the total organic carbon (TOC) volume fraction.
(2) The single gas-phase transport in kerogen involves slip corrected flow, Knudsen diffusion, surface diffusion, adsorption/desorption, while adsorption/desorption and surface diffusion are not included in the inorganic matrix due to its hydrophilic nature. (3) The flow in kerogen and inorganic matter is both strict spherical flow. The gas within kerogen first moves to the kerogen surface from its center and then enters the inorganic sphere. The two porous systems are coupled at their boundary, and the mass exchange between them is realized through flux and pressure continuity conditions at their interface without using a source term to describe kerogen's contribution. Real-gas effects for both free gas and adsorbed gas movement are considered. (4) The matrix element is evenly covered by a thin secondary fracture layer. The gas from the matrix element is instantaneously and uniformly distributed in one-half the fracture volume as the fracture permeability is far larger than that of matrix blocks (Ozkan et al. 2010).

Fracture and reservoir models
We use a linear flow model to couple fracture-reservoir gas flow because linear flow is the predominant flow regime in shale reservoirs with hydraulic fractures. Transient linear flow usually occurs in unconventional formations and may last for many years (Arevalo-Villagran et al. 2006;Tabatabaie and Pooladi-Darvish 2017). The major reason is that linear flow is associated with hydraulic fractures, and the well-fracture geometry results in linear flow behavior (Arevalo-Villagran et al. 2006). MFHWs drilled and fractured in shale gas plays exhibit long-term linear flow (Nobakht and Clarkson 2012). Linear flow tends to be the dominant flow regime during the production life of tight gas wells (Wattenbarger et al. 1998). Decline curves of some tight gas wells show that linear flow can last for 20 years, and the boundary flow is directly reached without experiencing pseudo-radial flow (Wattenbarger et al. 1998). Stright and Gordon (1983) also reported that the long-term production of fractured low-permeability gas wells can be  (2017) and Kang et al. (2011) 2 This is a kerogen and inorganic matrix parallel flow model without kerogeninorganic matrix inter porosity flow He et al. (2016) 3 The parallel flow model ignores the kerogen and inorganic matter mass exchange. The contribution from kerogen in pure inorganic matrix blocks is ignored for the combination of series and parallel flow. It is also difficult to determine the weighting factors of the three gas transport patterns for the combination of series and parallel flow Kang et al. (2011) and Zeng et al. (2019b) 4 Kerogen is surrounded by inorganic matrices to describe the dispersed distribution of kerogen. However, the kerogeninorganic matrix gas transfer is simply achieved by a source term, indicating that kerogen gas is instantaneously and evenly distributed in the inorganic matter domain. For tight shale matrices, it is not possible to achieve the instantaneous and homogeneous distribution of kerogen gas within inorganic matter. Besides, spherical matrices are more suitable for shale and coal compared with cubic matrices (Kucuk and Sawyer 1980) Olorode et al. (2017) and Huang et al. (2020) approximated by using linear flow equations. Bone et al. (2017) stated that linear post linear flow, which refers to one linear flow followed by another, is widely observed in the Permian Basin, the Eagle Ford Shale, Bakken Shale, and some formations in Oklahoma. Previously, seven-region linear flow models have been used to analyze several field cases (Zeng et al. 2017. The seven-region linear flow model can be used to estimate reservoir and fracture properties, such as formation permeability, hydraulic fracture length, permeability, secondary fracture conductivity, and reservoir sizes, from production histories (history matching) as long as we have some basic input data. By matching with a short-time production history, we can also predict the long-term well behavior. A validated linear flow model can be used to conduct extensive sensitivity analyses to examine the reservoir performance (Tabatabaie and Pooladi-Darvish 2017). By further considering the dispersed distribution of kerogen and avoiding the use of the instantaneous kerogen gas source term, the simulation would be closer to real reservoir conditions. Besides, due to the analytical nature of this model, the calculation is not time-consuming compared with semi-analytical and numerical models. Therefore, it can serve as a fast and effective tool to interpret field data, evaluate reservoir properties, and predict gas production. Due to the symmetry, the seven-region linear flow model can be described through one-eighth of a fracture drainage volume, as shown in Fig. 2b. Specifically, the fracture-reservoir model is derived under the following assumptions: (1) The MFHW is located at the center of a shale gas reservoir with no-flow outer boundaries. When the gas is released from the matrix, it flows within the natural fractures and finally enters the wellbore through hydraulic fractures. Regions 6 and 5 represent the upper/lower reservoir volumes beyond the fracture height involving 1-D vertical flow. Regions 3 and 4 are the reservoir volumes beyond the fracture tip with 1-D flow in the x-direction. Regions 2 and 1 are two inner reservoir regions with 1-D flow in the y-direction. Finally, the hydraulic fracture region is connected with region 1 through the fracture face with 1-D flow in the x-direction. The flow directions assumed in these flow regions of the seven-region linear flow model have been validated in the literature (Zeng et al. 2017). The composite linear flow model is practical for common unconventional tight formations. In this way, the four porous systems are finally connected in series.
(2) The reservoir can be homogeneous or heterogeneous with different properties in different linear flow regions. Between two hydraulic fractures, there is a no-flow boundary that describes the fracture interference. In both hydraulic and secondary fractures, Darcy's law is applied.  Fig. 2 a Schematic of the multi-continuum quadruple-porosity model and b schematic of the seven-region linear flow model (after Zeng et al. 2020) (3) The reservoir contains a single gas phase. The gravity effects can be ignored for the single-phase flow case. Besides, the low gas density and extremely tight texture of shale rocks make the gravity effects become marginal. The whole shale formation is under an isothermal condition during production.
Solutions are all analytical equations in the Laplace domain. The Stehfest algorithm (Stehfest 1970) and possible iteration processes (for heterogeneous reservoirs only) are finally used to obtain the real-time domain solutions. For convenience, the fundamental formulas are expressed in the dimensionless form. We will introduce the definition of dimensionless variables and gas transport mechanisms in each system, and give the derivation of solutions.

Dimensionless variables
In this paper, the subscripts i, sc, n , k, m, f, F, g, and app indicate the properties of the initial condition, standard condition, the n-th reservoir region, kerogen, inorganic matrix, secondary fractures, hydraulic fractures, gas, and apparent parameters. We let subscripts pk, pm, D, and ref represent kerogen and inorganic nano-pore properties, dimensionless parameters, and reference parameters for dimensionless variable definition. For the real-gas flow, the dimensionless pseudopressure is (Stalgorova and Mattar 2013) where k ref is the reference permeability in mD; h ref is the reference height in ft; p is the pressure in psi; p p is the pseudopressure in psi 2 /cP; T is the temperature in °R; q Ft is the total rate of the MFHW under the standard condition in Mscf/d. The pseudopressure is defined as (Ozkan et al. 2010) where g is the gas viscocity in cP, and Z is the Z-factor. Equation 2 deals with pressure-dependent properties and is used to linearize the diffusivity equation. The dimensionless time is given by where ref is the reference diffusivity in ft 2 /h; d ref is the reference length in ft; and t a is the pseudotime in h. The pseudotime is defined by (Anderson and Mattar 2007) (1) The diffusivity values of reference parameters, kerogen, inorganic matter, secondary fractures, and hydraulic fractures are given by where is the porosity, and c t is the total compressibility in psi −1 . Following Ozkan et al. (2010), ≈ i . Then we have their dimensionless expressions: And the dimensionless distances are where x , y , and z are reservoir and fracture dimensions in ft, as shown in Fig. 2b; w is the fracture width in ft; h F and h are fracture height and reservoir height in ft; and r is the distance to the spherical matrix center in ft.

Gas transport in the kerogen core
Kerogen is a nanoporous organic material that is interspersed within the inorganic matrix (Akkutlu and Fathi 2012), which leads to multiscale gas transport phenomena. Besides, as nano-pores in kerogen have large surface areas with the affinity to methane, a certain portion of methane is adsorbed in kerogen, inducing gas desorption and surface diffusion during gas extraction (Wasaki and Akkutlu 2015). In our kerogen gas transport model, the slip corrected flow, Knudsen diffusion, gas adsorption/desorption, and surface diffusion are considered. The mass balance equation of the gas flow in kerogen is given by where kg is the kerogen gas density in lbm/ft 3 ; scg is the gas density under the standard condition in lbm/ft 3 ; the modified real-gas equilibrium gas volume is (Zeng et al. 2020) and the total gas velocity in kerogen can be written as Here V L and p L are the Langmuir volumetric concentration constant (scf/cf) and Langmuir pressure (psi) respectively; v kfg is the kerogen free gas velocity in mD·psi/(cP·ft); and v ks is the kerogen surface diffusion velocity in mD·psi/(cP·ft). They can be expressed as follows (Beskok and Karniadakis 1999;Zeng et al. 2020) where Here Kn is the Knudsen number; k kfg is the free gas apparent permeability in mD; 1 is a coefficient (9.4127 × 10 13 mD/ ft 2 ) converting ft 2 into mD; r pk is the nano-pore radius of kerogen in ft; k is the kerogen tortuosity; k is the coefficient for the permeability correction model (Beskok and Karniadakis 1999); 2 is a unit conversion factor, 158 mD·psi·d/ (ft 2 ·cP) (Ertekin et al. 1986;Zeng et al. 2017); D s is the diffusivity for surface diffusion in ft 2 /D; c kg is the kerogen gas compressibility in psi −1 ; and C k is the kerogen gas molar concentration in mol/ft 3 and is given by The unit of M here is lbm/mol, and Z sc = 1 . Therefore, Eq. 23b can be further written as Finally, combining Eqs. 22-28, we have the following mass balance equation where Converting Eq. 29 into the pseudopressure form, we obtain where the apparent permeability and compressibility are Converting Eq. 31 into the dimensionless form yields Assuming that the flow direction in secondary fractures is in the x-direction, we have the following initial and boundary conditions for Eq. 34 In Eq. 36, the pressure in the kerogen core center is definitely a finite value. By defining Therefore, the dimensionless pseudopressure is obtained Equation 42 is the solution of the kerogen core.

Gas transport in the inorganic matrix sphere
In inorganic matrices, the pore size is larger than that of kerogen although the inorganic porosity is normally lower than kerogen porosity (Wang and Reed 2009;Kang et al. 2011;Wang et al. 2014). The permeability correction model (Beskok and Karniadakis 1999) is still applicable here. Different from kerogen, the hydrophilic nature of clay minerals in the inorganic matrix reduces the gas adsorption capacity (Zhang et al. 2012). The connate water molecules are sorbing to the hydrophilic pore surfaces of the inorganic matrix, resulting in few sorption sites for methane molecules. Therefore, the surface diffusion and gas desorption effects are ignored in the inorganic matrix, which is in accordance with the assumption made by other related studies (Akkutlu and Fathi 2012;Cao et al. 2016). Similarly, the mass balance equation of the gas flow in inorganic matter in the pseudopressure form is where Converting Eq. 43 into the dimensionless form, we have Assuming the flow direction in secondary fractures is in the x-direction, we have the following initial and boundary conditions for Eq. 45

Gas transport in secondary fractures and hydraulic fractures
In secondary and hydraulic fractures, gas transport simply obeys Darcy's law. The matrix inflow mass rate per unit secondary fracture volume per unit time is given by (Zeng et al. 2017) Therefore, we have the mass balance equation for secondary fractures The dimensionless pseudopressure form of Eq. 56 is Converting Eq. 57 into the Laplace form yields 2p pfD Then, the composite linear flow model of this paper will be introduced. The schematic of the flow regions and flow directions is given in Fig. 2b. The equations for linear flow regions are similar to those of the classical seven-region linear flow model (Zeng et al. 2017). To ensure the completeness of the paper and avoid the confusion generated by the similar symbols of previous seven-region models, we follow Zeng et al. (2017) to demonstrate the linear flow equations in each region.

Region 6
Region 6 involves 1-D flow in the vertical direction. The diffusivity equation, no-flow boundary and pressure continuity conditions are Solving Eqs. 61-63 yields Equation 64 is the solution of region 6 and will be used in the diffusivity equations of regions 2 and 4.

Region 5
Region 5 involves 1-D flow in the vertical direction. The diffusivity equation, no-flow boundary and pressure continuity conditions are (60)

Solving Eqs. 65-67 yields
Equation 68 is the solution of region 5 and will be used in the diffusivity equations of regions 1 and 3.

Region 4
Region 4 2p pfD4 p pfD4 p pfD4 Equation 72 is the solution of region 4 and will be used in the diffusivity equation of region 2.

Region 2
Region 2 involves 1-D flow in the y-direction. The diffusivity equation, no-flow boundary and pressure continuity conditions are 2p pfD3 p pfD3 p pfD3 2p pfD2 2p pfD1 Equation 87 is the solution of region 1 and will be used in the diffusivity equation of the hydraulic fracture region.

Hydraulic fracture region
The hydraulic fracture region involves 1-D flow in the x-direction. The diffusivity equation and boundary conditions are Here the subscript j represents the j-th hydraulic fracture. In Eq. 90, the bulk secondary fracture permeability of region 1 (k f1b ) is used, while in other regions we equivalently use the intrinsic secondary fracture permeability as the fracturematrix volume ratios are the same in all reservoir regions. The bulk secondary fracture permeability is defined below for spherical matrix elements (Apaydin et al. 2012), where k f is the secondary fracture intrinsic permeability in mD; h f is the secondary fracture layer thickness in ft; and r m is the inorganic matrix radius in ft. Solving Eqs. 90-92 yields where Equation 94 is the final solution for constant flow rate cases. For constant bottom-hole pressure conditions, the definition of dimensionless pressure is different as given below 2p pFD p pFD where p wf is the bottom-hole flowing pressure in psi. By assuming p i − p wf = p i − p wf ref and ignoring the pressure drop in the horizontal wellbore for gas reservoirs, we have Then the final solution of constant pressure cases is Note that, in Eq. 98, p pFDj is obtained from Eq. 94. And the total rate of an MFHW is expressed by Consequently, the cumulative production is given by According to Callard and Schenewerk (1995) and Chaudhry (2003), under constant pressure production conditions, the dimensionless rate and cumulative production are defined as follows: Here Q is the cumulative production in Mscf.

Model validation
In this section, this model is verified against a published analytical model with single-continuum spherical matrix blocks (Zeng et al. 2017). The published model has been validated by comparing with commercial well-testing software KAPPA (Zeng et al. 2017). To perform a fair comparison, the proposed model must be degraded into the published one by assuming that the kerogen core size is equal to the matrix block size (r m = r k ) . The input parameters are listed in Table 2. Figure 3 shows that a good agreement between the calculation results of the two models has been reached. Because the results are from the degraded form of our model, flow regime diagnostic is not conducted here. Flow regimes will be analyzed in the discussion section using the original form of the proposed model. In fact, there could be alternative flow directions in these linear flow regions. Stalgorova and Mattar (2013) and Zeng et al. (2017Zeng et al. ( , 2018 found that the current flow direction assumptions of the seven-region linear flow model are accurate enough for reservoir simulation if the following conditions are met (Stalgorova and Mattar 2013; Zeng et al. 2017Zeng et al. , 2018Zeng et al. , 2020: (1) the width of region 1 is no smaller than 10% of the half fracture spacing (y 1 ≥ 0.1y 2 ) ; (2) the fracture length is at least 10% of the reservoir width (x F ≥ 0.1x e ) ; (3) at least 60% of the formation is vertically penetrated by hydraulic fractures (h F ≥ 0.6h) . Note that if the fracture height is too small (h F < 0.6h) , the hydraulic fracturing process is regarded as unsuccessful, and there is no need to conduct further analyses. These conditions are generally met by reservoirs with MFHWs (Stalgorova and Mattar 2013). Next, we will provide a field application example to further ensure the reliability of this model.

History matching of a Barnett Shale MFHW
To further ensure its reliability and applicability, an application example of history matching of Barnett Shale field data is provided, as shown in Fig. 4 Table 3. For unconventional tight formations, the outer spherical block size normally ranges from 3 ft to 6 ft, and the fracture thickness is between 5 × 10 −4 ft to 5 × 10 −3 ft (Apaydin et al. 2012). As mentioned before, the spherical matrix element size can easily be up to five orders of magnitude larger than the SEM images in reservoir simulation (Olorode et al. 2017). Then the kerogen core size is calculated based on the kerogen volume fraction. As demonstrated before, kerogen normally has a smaller pore size and a larger porosity value compared with those of the inorganic matrix. The nano-pore size is selected to satisfy matrix permeability. In this field case, only the total matrix porosity and matrix permeability are available, we, therefore, assume that kerogen and inorganic matrices have the same porosity and pore radius. In our later sensitivity analyses, the input porosity of kerogen is larger than inorganic matter porosity, and the pore radius of kerogen is smaller than the inorganic pore radius. It is also worth noting that the water saturation in the original case is 0.3. Therefore, the effective porosity for gas flow is 0.042. However, when we calculate adsorption/desorption effects and surface diffusion, we need to use the 0.06 porosity as the volume ratio of the rock grains that store the adsorbed-phase gas will not change with water saturation. The tortuosity is calculated from the porosity according to the method of Chen et al. (2015): = −0.5 . And the surface diffusion coefficient is selected from Akkutlu and Fathi (2012). The conductivity values of natural and hydraulic fractures are evaluated through the matching process. The fracture closure stress of this case is between 929 psi and 3379 psi (Yu and Sepehrnoori 2014). The estimated hydraulic fracture conductivity (0.17 mD·ft) and unpropped natural fracture conductivity (0.003 mD·ft) are reasonable compared with the laboratory measurements (0.037-0.23 mD·ft under 4000 psi closure stress conditions and 0.0017-0.01 mD·ft at 3000 psi confining pressure for hydraulic and natural fractures respectively) (Zhang 2016;Zhang et al. 2016). Therefore, this model is practical and can serve as an effective tool for MFHW performance prediction and evaluation in shale gas reservoirs.

Flow regimes
The flow regimes observed in our simulation are introduced in this section, as shown in Fig. 5. Key input parameters are listed in Table 4. Five typical flow regimes can be seen from the log-log dimensionless pseudopressure and derivative curves, including the fracture linear flow regime, matrix-fracture inter-porosity flow regime, transient bilinear flow regime, transient fracture interference regime, and the boundary dominant flow regime. This is in accordance with the single-continuum spherical matrix solution (Zeng et al. 2020).

Regime I:
This regime is an early linear flow regime in the hydraulic fracture. The first two parallel straight lines in the dimensionless pseudopressure and derivative curves represent the fracture linear flow period which is evidenced by the 1/2 slope. Note that the slopes are utilized to identify flow regimes. From point A to point B, the fluids within hydraulic fracture systems contribute to production. The section from point B to point C is a transient regime between fracture linear flow and matrix-fracture inter-porosity flow.

Regime II:
This regime is an inter-porosity flow regime involving the mass exchange between matrices and natural fractures. During this period, the derivative curve forms a concave-up curve. The pressure difference between shale matrices and natural fractures will first increase and then decrease due to their permeability difference. The section from point C to point D indicates that matrix gas supplies the fracture fluid flow, and the derivative increment slows down. This period ends when the matrix and natural fracture systems reach a dynamic equilibrium (pseudo-steady) state (Kim and Lee 2015). Calculation results Field gas rates Fig. 4 Comparison of simulation results and Barnett Shale production data Regime III: This regime is a transient bilinear flow regime as evidenced by a slope of 1/4 in the derivative curve. The section from point D to point E indicates that simultaneous linear flow in the hydraulic fracture and to the hydraulic fracture plane (bilinear flow) exists. With hydraulic fracture conductivity increases, the slope of this transient regime will become larger. If the hydraulic fracture conductivity approaches an equivalent level of infinite-conductivity fractures, the slope of this period becomes 1/2, indicating an inner formation linear flow response, which is in accordance with Kim and Lee (2015). If the matrix permeability is extremely low, this regime can be masked by the matrixfracture inter-porosity flow regime that occupies a long production period before the boundary flow.

Regime IV:
This regime is a transition regime representing the fracture interference (the section from point E to point F). In this regime, the pressure wave has reached the no-flow boundary between two hydraulic fractures but has not reached the outer reservoir boundary. The slope of this regime in the derivative curve of the given case is between 1/2 and 1, which agrees with Kim and Lee (2015). The range of slopes indicates that it is between the formation linear flow and boundary flow. The duration of this period depends on the relative magnitude of hydraulic fracture length, fracture spacing, and the distances from the horizontal wellbore to external reservoir boundaries.

Regime V:
This regime is the boundary-dominated flow regime. In the section from point F to point G, the pressure wave induced by production has encountered the external reservoir boundary. The dimensionless pseudopressure and derivative curves gradually converge. Finally, they overlap each other with a unit-slope, which indicates that a pseudosteady state has been achieved.

Dimensionless pseudopressure and rate profiles of the matrix
In this section, the radial dimensionless pseudopressure distribution and the rate profile of a matrix block located in region 1 (x D = 1/2x 1D and y D = 1/2y 1D ) are presented, as shown in Figs. 6 and 7. The radii of the kerogen core and inorganic sphere are 3 ft and 6 ft respectively according to Apaydin et al. (2012) and the TOC volume fraction. The dimensionless pseudopressure reveals the magnitude of the pressure drop. When the dimensionless pseudopressure value is too small, there are fluctuations in the Stehfest numerical inversion process. Therefore, these fluctuant points are ignored in our analysis and are not presented in Fig. 6. The pressure profile shows that even the dynamic equilibrium (pseudo-steady) state has been achieved in the inorganic matrix (   . Even the whole matrix reaches the dynamic equilibrium state, the dimensionless rate profile ( t = 13.6 years) shows that the kerogen core ( r ≤ 3.0 ft) only contributes limited gas compared with the inorganic matrix, as shown in Fig. 7. This indicates that when the inorganic matrix has been nearly depleted, a certain portion of the kerogen core is still not effectively drained. And most of the gas is produced from the inorganic matrix layer and secondary fractures for a relatively long period. By comparing Figs. 5 and 6, one can also find that there is a delay in arrival of the dynamic equilibrium state for the whole matrix block that is 62.5 ft (y = 1/2y 1 ) away from the hydraulic fracture because the dynamic equilibrium state (the end of the concave regime) occurs earlier on the fracture face (Fig. 5). Thus, even the well has produced for a certain period, only limited matrix volumes in matrix blocks far away from the hydraulic fracture contribute to well production.

Effects of the matrix pore size
The pore size determines organic and inorganic matrix apparent permeability. In this study, the range of the pore sizes in inorganic matrices and kerogen is selected from the literature (Wu et al. 2016). As mentioned before, the average pore radius in kerogen is normally much smaller than that of the inorganic matrices (Kang et al. 2011). The average inorganic pore size is larger because inorganic matter contains more microfractures. The limiting case in our analysis assumes that the kerogen and inorganic matrix pore sizes are equal. Figure 8 shows that the variation of the kerogen pore radius mainly influences the late matrix-fracture inter-porosity flow regime and slightly affects the transient bilinear flow regime. The larger the kerogen pore radius is, the lower the dimensionless pseudopressure and derivative values in the late concave regime will be. And no obvious influence can be seen in the dimensionless rate and cumulative production curves which are not shown here. The effects of the inorganic pore radius are more noticeable. Figure 9 demonstrates that the inorganic pore radius controls the late fracture linear flow regime, concave inter-porosity flow regime, and the transient bilinear flow regime. The bigger the inorganic pore radius is, the earlier the concave and bilinear flow regimes will occur, and the lower the dimensionless pseudopressure will be. When the inorganic matrix and kerogen pore radii are both small (5 × 10 −9 ft), the mass exchange between matrices and natural fractures takes a long time before the dynamic equilibrium state, and the bilinear flow regime disappears. Actually, if the permeability is further lower, the fracture interference regime may also be masked. Figures 9 and 10 also report that with the inorganic pore radius increases, its influences on dimensionless pressure and rate curves become smaller. In terms of cumulative production, only the limiting case with the smallest inorganic Pressure (rpk = 5e-9 ft, rpm = 2e-7 ft) Derivative (rpk = 5e-9 ft, rpm = 2e-7 ft) Pressure (rpk = 2e-8 ft, rpm = 2e-7 ft) Derivative (rpk = 2e-8 ft, rpm = 2e-7 ft) Pressure (rpk = 5e-8 ft, rpm = 2e-7 ft) Derivative (rpk = 5e-8 ft, rpm = 2e-7 ft) Pressure (rpk = 2e-7 ft, rpm = 2e-7 ft) Derivative (rpk = 2e-7 ft, rpm = 2e-7 ft) p pwD , dp pwD /dlnt D t D Fig. 8 Effects of the kerogen pore radius on pressure behavior pore radius has slightly smaller cumulative production, as shown in Fig. 11. It is worth noting that the ultimate cumulative production for all cases is the same. This is because the pore radius only depicts apparent permeability, the porosity of each porous system and the adsorption capability of kerogen are the same in all cases. Note that the leading zero in a decimal less than 1 is omitted in the figures.

Fig. 21
Effects of the Langmuir pressure on pressure behavior the dimensionless rates increase. However, these effects are not as significant as those of kerogen porosity. When the Langmuir volume changes from 0 scf/cf to 100 scf/cf, the late-time cumulative production increases almost linearly, as shown in Fig. 20. Another important parameter is the Langmuir pressure. Figures 21, 22 and 23 demonstrate the influence of Langmuir pressure on well responses. Here we keep the Langmuir volume as a constant and change the Langmuir pressure. The range of Langmuir pressure is selected from Chen et al. (2018). It is shown that with the increase in Langmuir pressure, the dimensionless pseudopressure and derivative at late times become slightly lower, indicating a delay in the occurrence of boundary flow. For dimensionless rate responses, the gas rate decreases slower at late times for cases with larger Langmuir pressure. The

Effects of the Langmuir volume and Langmuir pressure
The adsorbed gas release process in coal and shale reservoirs is normally described by the Langmuir isotherm theory. The Langmuir volume V L in the Langmuir isotherm equation is a critical parameter for reservoir performance evaluation. In this section, the range of the Langmuir volume is selected from the literature (Ali 2012) and converted into the unit of scf/cf. Similar to kerogen porosity, the Langmuir volume effects become noticeable from the late matrix-fracture inter-porosity flow regime to the end of well life, as shown in Figs. 18 and 19. The adsorbed-phase gas serves as an additional gas source for well production. With the Langmuir volume increases, the dimensionless pressure and derivative values of the affected regimes decrease, while ultimate cumulative production increases when Langmuir pressure turns larger. Our results are in accordance with the numerical results of Cao et al. (2016). Theoretically, a higher Langmuir pressure value represents a lower adsorbability at identical pressure, while the adsorbability difference of cases with variable Langmuir pressure decreases with the increase of reservoir pressure. The initial reservoir pressure is high, and the adsorbability difference among cases with different Langmuir pressure is relatively small. Under the depleted reservoir condition, the reservoir pressure is very low. Cases with larger Langmuir pressure can release more gas for production although their initial amount of adsorbed gas is smaller than cases with lower Langmuir pressure.

Contributions of different linear flow regions on gas production
In essence, the reservoir-fracture flow model is a sevenregion linear flow model. The contributions of individual flow regions are finally analyzed. The sizes of different regions (length × width × height) are: 400 ft × 100 ft × 100 ft for regions 1 and 2, 150 ft × 100 ft × 100 ft for regions 3 and 4, and 550 ft × 100 ft × 25 ft for regions 5 and 6. Other input parameters are the same as those in Table 4. The log-log curve (Fig. 24) shows the cumulative contributions of different regions. It can be seen that region 1, which is the closest region to the hydraulic fracture, contributes to production first. This is followed by regions 5, 3, and 2 which are directly connected to region 1. Among the three regions, the contribution of region 5 occurs first because region 3's contribution starts only when the production-induced pressure wave propagates to the fracture tip (x = x F ) , and region 2's contribution occurs latest when the pressure wave reaches the interface between regions 1 and 2. Region 5 is directly connected to region 1 along the hydraulic fracture and hence its contribution starts earlier than regions 2 and 3. In a similar manner, the contribution of region 6 starts slightly earlier than that of region 4. The ultimate cumulative contribution depends on the reservoir volume of each region. Therefore, the ultimate cumulative production of different regions has the following relationship: region 1 = region 2 > region 3 = region 4 > region 5 = region 6, which is in accordance with their volumes.

Conclusions and final remarks
A new analytical model for shale gas extraction considering dispersed distribution of kerogen and complex gas flow mechanisms is built. In essence, the gas inflow from the dispersed kerogen to inorganic matter is realized through flux continuity conditions at the kerogen core surface instead of using the kerogen source term in the inorganic matrix flow equation, which is closer to the actual gas transport situation. Based on the simulation results, we can draw the following conclusions: (1) The inorganic matrix plays as the major gas source during production. When the whole matrix reaches the dynamic equilibrium drainage state, kerogen still contributes limited gas to production compared with inorganic matter. (2) The pore radius of kerogen mainly affects the late matrix-fracture inter-porosity flow regime and slightly influences the following transient bilinear flow regime. In contrast, the inorganic pore radius dominates the shape and duration of the late fracture linear flow regime, concave inter-porosity flow regime, and the transient bilinear flow regime. The cumulative production curve is not sensitive to the pore radius change unless the pore radius is extremely low when the matrix-fracture inter-porosity regime masks the transient bilinear flow regime and even the transient fracture interference regime. (3) In terms of porosity, kerogen porosity's influences start from the late matrix-fracture inter-porosity flow regime, while inorganic matter porosity manipulates almost all flow regimes except the early-mid fracture linear flow regime. When the porosity increases, both pressure and derivative values decrease, and the dimensionless rates increase with a delay in the arrival of boundary flow regimes. The porosity, especially the inorganic porosity, also significantly influences the cumulative production. Overall, the influence of kerogen porosity is less important than that of inorganic porosity. (4) Gas desorption from kerogen severs as an additional gas source during production. The Langmuir volume influences the well responses in a similar manner as the Langmuir pressure. (5) The occurrence of contributions from different linear flow regions depends on their relative location, while the cumulative contribution is dominated by their sizes.