1 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 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 low-gas-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 (Guo et al. 2017). 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.

Fig. 1
figure 1

Distribution of kerogen/organic matter (OM) and inorganic matter (iOM) in shale samples: a intact shale and b a 2D SEM image (figures are obtained from He et al. 2019 and Tahmasebi et al. 2015 and reprinted with the permission from He et al. Copyright Springer Nature 2019 and Tahmasebi et al. Copyright Springer Nature 2015)

Full size image

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, 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. 2015, 2018; Zeng 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, 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 (Sheng et al. 2015; 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.

Table 1 Features of matrix gas transfer models
Full size table

2 Conceptual models

2.1 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:

Fig. 2
figure 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)

Full size image
  1. (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. (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. (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. (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).

2.2 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 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, 2018). 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. (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. (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.

  3. (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.

3 Formulation of the conceptual models

3.1 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)

$$p_{\text{pD}} = \frac{{k_{\text{ref}} h_{\text{ref}} }}{{1422q_{\text{Ft}} T}}\left[ {p_{\text{p}} \left( {p_{\text{i}} } \right) - p_{\text{p}} \left( p \right)} \right],$$
(1)

where \(k_{\text{ref}}\) is the reference permeability in mD; \(h_{\text{ref}}\) is the reference height in ft; \(p\) is the pressure in psi; \(p_{\text{p}}\) is the pseudopressure in psi2/cP; \(T\) is the temperature in °R; \(q_{\text{Ft}}\) is the total rate of the MFHW under the standard condition in Mscf/d. The pseudopressure is defined as (Ozkan et al. 2010)

$$p_{\text{p}} \left( p \right) = 2\mathop \smallint \limits_{{p_{\text{ref}} }}^{p} \frac{k}{{k_{\text{i}} }}\frac{p}{{\mu_{\text{g}} Z}}{\text{d}}p,$$
(2)

where \(\mu_{\text{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

$$t_{\text{D}} = \frac{{\eta_{\text{ref}} t_{\text{a}} }}{{d_{\text{ref}}^{2} }},$$
(3)

where \(\eta_{\text{ref}}\) is the reference diffusivity in ft2/h; \(d_{\text{ref}}\) is the reference length in ft; and \(t_{\text{a}}\) is the pseudotime in h. The pseudotime is defined by (Anderson and Mattar 2007)

$$t_{\text{a}} = \mu_{\text{gi}} c_{\text{ti}} \mathop \smallint \limits_{0}^{t} \frac{{{\text{d}}t}}{{\mu_{\text{g}} c_{\text{t}} }}.$$
(4)

The diffusivity values of reference parameters, kerogen, inorganic matter, secondary fractures, and hydraulic fractures are given by

$$\eta_{\text{ref}} = \frac{{2.637 \times 10^{ - 4} k_{\text{ref}} }}{{\phi_{\text{ref}} \mu_{\text{ref}} c_{\text{tref}} }},$$
(5)
$$\eta_{{{\text{k}}n}} = \left( {\frac{{2.637 \times 10^{ - 4} k_{\text{appk}} }}{{\phi_{\text{k}} \mu_{\text{kgi}} c_{\text{kappi}} }}} \right)_{n} ,$$
(6)
$$\eta_{{{\text{m}}n}} = \left( {\frac{{2.637 \times 10^{ - 4} k_{\text{appm}} }}{{\phi_{\text{m}} \mu_{\text{mgi}} c_{\text{tmi}} }}} \right)_{n} ,$$
(7)
$$\eta_{{{\text{f}}n}} = \left( {\frac{{2.637 \times 10^{ - 4} k_{\text{f}} }}{{\phi_{\text{f}} \mu_{\text{fgi}} c_{\text{tfi}} }}} \right)_{n} ,$$
(8)
$$\eta_{\text{F}} = \frac{{2.637 \times 10^{ - 4} k_{\text{F}} }}{{\phi_{\text{F}} \mu_{\text{Fgi}} c_{\text{tFi}} }},$$
(9)

where \(\phi\) is the porosity, and \(c_{\text{t}}\) is the total compressibility in psi−1. Following Ozkan et al. (2010), \(\eta \approx \eta_{\text{i}}\). Then we have their dimensionless expressions:

$$\eta_{{{\text{kD}}n}} = \eta_{{{\text{k}}n}} /\eta_{\text{ref}} ,$$
(10)
$$\eta_{{{\text{mD}}n}} = \eta_{{{\text{m}}n}} /\eta_{\text{ref}} ,$$
(11)
$$\eta_{{{\text{fD}}n}} = \eta_{{{\text{f}}n}} /\eta_{\text{ref}} ,$$
(12)
$$\eta_{\text{FD}} = \eta_{\text{F}} /\eta_{\text{ref}} .$$
(13)

And the dimensionless distances are

$$x_{{1{\text{D}}}} = x_{\text{F}} /d_{\text{ref}} ,$$
(14)
$$x_{\text{eD}} = x_{\text{e}} /d_{\text{ref}} ,$$
(15)
$$y_{{1{\text{D}}}} = y_{1} /d_{\text{ref}} ,$$
(16)
$$y_{{2{\text{D}}}} = y_{2} /d_{\text{ref}} ,$$
(17)
$$z_{{1{\text{D}}}} = z_{1} /d_{\text{ref}} = h_{\text{F}} /\left( {2d_{\text{ref}} } \right),$$
(18)
$$z_{{2{\text{D}}}} = z_{2} /d_{\text{ref}} = h/\left( {2d_{\text{ref}} } \right),$$
(19)
$$w_{\text{D}} = w_{\text{F}} /d_{\text{ref}} ,$$
(20)
$$r_{\text{D}} = r/d_{\text{ref}} \;\left( {0 \le r \le r_{\text{m}} } \right),$$
(21)

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_{\text{F}}\) and \(h\) are fracture height and reservoir height in ft; and \(r\) is the distance to the spherical matrix center in ft.

3.2 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

$$- \frac{1}{{r^{2} }}\frac{\partial }{\partial r}\left( {r^{2} \rho_{\text{kg}} v_{\text{k}} } \right) = \frac{{\partial \left( {\rho_{\text{kg}} \phi_{\text{k}} } \right)}}{\partial t} + \frac{{\partial \left[ {\left( {1 - \phi_{\text{k}} } \right)\rho_{\text{scg}} V_{\text{Esc}} } \right]}}{\partial t},$$
(22)

where \(\rho_{\text{kg}}\) is the kerogen gas density in lbm/ft3; \(\rho_{\text{scg}}\) is the gas density under the standard condition in lbm/ft3; the modified real-gas equilibrium gas volume is (Zeng et al. 2020)

$$V_{\text{Esc}} = V_{\text{L}} \frac{p/\sqrt Z }{{p/\sqrt Z + p_{\text{L}} }},$$
(23a)

and the total gas velocity in kerogen can be written as

$$v_{\text{k}} = v_{\text{kfg}} + v_{\text{ks}} .$$
(23b)

Here \(V_{\text{L}}\) and \(p_{\text{L}}\) are the Langmuir volumetric concentration constant (scf/cf) and Langmuir pressure (psi) respectively; \(v_{\text{kfg}}\) is the kerogen free gas velocity in mD·psi/(cP·ft); and \(v_{\text{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)

$$v_{\text{kfg}} = - \frac{{k_{\text{kfg}} }}{{\mu_{\text{kg}} }}\frac{{\partial p_{\text{k}} }}{\partial r},$$
(24)

where

$$k_{\text{kfg}} = \xi_{1} \frac{{\phi_{\text{k}} }}{{\tau_{\text{k}} }}\frac{{r_{\text{pk}}^{2} }}{8}\left( {1 + \alpha_{\text{k}} Kn_{\text{k}} } \right)\left( {1 + \frac{{4Kn_{\text{k}} }}{{1 - bKn_{\text{k}} }}} \right) = \xi_{1} \frac{{\phi_{\text{k}} }}{{\tau_{\text{k}} }}\frac{{r_{\text{pk}}^{2} }}{8}\left( {1 + \alpha_{\text{k}} \frac{{\lambda_{\text{k}} }}{{r_{\text{pk}} }}} \right)\left( {1 + \frac{{4\frac{{\lambda_{\text{k}} }}{{r_{\text{pk}} }}}}{{1 - b\frac{{\lambda_{\text{k}} }}{{r_{\text{pk}} }}}}} \right),$$
(25)

and

$$v_{\text{ks}} = - \xi_{2} \frac{{MD_{\text{s}} }}{{\rho_{\text{kg}} }}\frac{{\partial C_{\text{k}} }}{\partial r} = - \xi_{2} \frac{{MD_{\text{s}} }}{{\rho_{\text{kg}} }}\frac{{\left( {1 - \phi_{\text{k}} } \right)\rho_{\text{scg}} }}{M}\frac{{V_{\text{L}} p_{\text{L}} }}{{2\sqrt {Z_{\text{k}} } }}\frac{{\left( {1 + p_{\text{k}} c_{\text{kg}} } \right)}}{{\left( {\frac{{p_{\text{k}} }}{{\sqrt {Z_{\text{k}} } }} + p_{\text{L}} } \right)^{2} }}\frac{{\partial p_{\text{k}} }}{\partial r} = - \xi_{2} \frac{{D_{\text{s}} \left( {1 - \phi_{\text{k}} } \right)p_{\text{sc}} Z_{\text{k}} T}}{{p_{\text{k}} T_{\text{sc}} }}\frac{{V_{\text{L}} p_{\text{L}} }}{{2\sqrt {Z_{\text{k}} } }}\frac{{\left( {1 + p_{\text{k}} c_{\text{kg}} } \right)}}{{\left( {\frac{{p_{\text{k}} }}{{\sqrt {Z_{\text{k}} } }} + p_{\text{L}} } \right)^{2} }}\frac{{\partial p_{\text{k}} }}{\partial r}.$$
(26)

Here Kn is the Knudsen number; \(k_{\text{kfg}}\) is the free gas apparent permeability in mD; \(\xi_{1}\) is a coefficient (9.4127 × 1013 mD/ft2) converting ft2 into mD; \(r_{\text{pk}}\) is the nano-pore radius of kerogen in ft; \(\tau_{\text{k}}\) is the kerogen tortuosity; \(\alpha_{\text{k}}\) is the coefficient for the permeability correction model (Beskok and Karniadakis 1999); \(\xi_{2}\) is a unit conversion factor, 158 mD·psi·d/(ft2·cP) (Ertekin et al. 1986; Zeng et al. 2017); \(D_{\text{s}}\) is the diffusivity for surface diffusion in ft2/D; \(c_{\text{kg}}\) is the kerogen gas compressibility in psi−1; and \(C_{\text{k}}\) is the kerogen gas molar concentration in mol/ft3 and is given by

$$C_{\text{k}} = \frac{{\left( {1 - \phi_{\text{k}} } \right)\rho_{\text{scg}} V_{\text{Esc}} }}{M}.$$
(27)

The unit of \(M\) here is lbm/mol, and \(Z_{\text{sc}} = 1\). Therefore, Eq. 23b can be further written as

$$v_{\text{k}} = - \frac{{k_{\text{kfg}} }}{{\mu_{\text{kg}} }}\frac{{\partial p_{\text{k}} }}{\partial r}\left[ {1 + \xi_{2} \frac{{\mu_{\text{kg}} D_{\text{s}} \left( {1 - \phi_{\text{k}} } \right)p_{\text{sc}} Z_{\text{k}} T}}{{k_{\text{kfg}} p_{\text{k}} T_{\text{sc}} }}\frac{{V_{\text{L}} p_{\text{L}} }}{{2\sqrt {Z_{\text{k}} } }}\frac{{\left( {1 + p_{\text{k}} c_{\text{kg}} } \right)}}{{\left( {\frac{{p_{\text{k}} }}{{\sqrt {Z_{\text{k}} } }} + p_{\text{L}} } \right)^{2} }}} \right].$$
(28)

Finally, combining Eqs. 22–28, we have the following mass balance equation

$$\begin{aligned} & \frac{1}{{r^{2} }}\frac{\partial }{\partial r}\left\{ {r^{2} \frac{{p_{\text{k}} }}{{Z_{\text{k}} }}\frac{{k_{\text{kfg}} }}{{\mu_{\text{kg}} }}\frac{{\partial p_{\text{k}} }}{\partial r}\left[ {1 + \xi_{2} \frac{{\mu_{\text{kg}} D_{\text{s}} \left( {1 - \phi_{\text{k}} } \right)p_{\text{sc}} Z_{\text{k}} T}}{{k_{\text{kfg}} p_{\text{k}} T_{\text{sc}} }}\frac{{V_{\text{L}} p_{\text{L}} }}{{2\sqrt {Z_{\text{k}} } }}\frac{{\left( {1 + p_{\text{k}} c_{\text{kg}} } \right)}}{{\left( {\frac{{p_{\text{k}} }}{{\sqrt {Z_{\text{k}} } }} + p_{\text{L}} } \right)^{2} }}} \right]} \right\} \\ & \quad = c_{\text{kg}} \phi_{\text{k}} \frac{{p_{\text{k}} }}{{Z_{\text{k}} }}\frac{{\partial p_{\text{k}} }}{\partial t} + \frac{{p_{\text{sc}} T}}{{T_{\text{sc}} }}\left( {1 - \phi_{\text{k}} } \right)\frac{{V_{\text{L}} p_{\text{L}} }}{{2\sqrt {Z_{\text{k}} } }}\frac{{\left( {p_{\text{k}} c_{\text{kg}} + 1} \right)}}{{\left( {\frac{{p_{\text{k}} }}{{\sqrt {Z_{\text{k}} } }} + p_{\text{L}} } \right)^{2} }}\frac{{\partial p_{\text{k}} }}{\partial t}, \\ \end{aligned}$$
(29)

where

$$c_{\text{kg}} = \frac{1}{{p_{\text{k}} }} - \frac{1}{{Z_{\text{k}} }}\frac{{{\text{d}}Z_{\text{k}} }}{{{\text{d}}p_{\text{k}} }}.$$
(30)

Converting Eq. 29 into the pseudopressure form, we obtain

$$\frac{1}{{r^{2} }}\frac{\partial }{\partial r}\left( {r^{2} k_{\text{appki}} \frac{{\partial p_{\text{pk}} }}{\partial r}} \right) = \mu_{\text{kg}} c_{\text{kapp}} \phi_{\text{k}} \left( {\frac{{k_{\text{appki}} }}{{k_{\text{appk}} }}} \right)\frac{{\partial p_{\text{pk}} }}{\partial t},$$
(31)

where the apparent permeability and compressibility are

$$k_{\text{appk}} = k_{\text{kfg}} \left[ {1 + \xi_{2} \frac{{\mu_{\text{kg}} D_{\text{s}} \left( {1 - \phi_{\text{k}} } \right)p_{\text{sc}} Z_{\text{k}} T}}{{k_{\text{kfg}} p_{\text{k}} T_{\text{sc}} }}\frac{{V_{\text{L}} p_{\text{L}} }}{{2\sqrt {Z_{\text{k}} } }}\frac{{\left( {1 + p_{\text{k}} c_{\text{kg}} } \right)}}{{\left( {\frac{{p_{\text{k}} }}{{\sqrt {Z_{\text{k}} } }} + p_{\text{L}} } \right)^{2} }}} \right],$$
(32)
$$c_{\text{kapp}} = c_{\text{kg}} + \frac{{p_{\text{sc}} TZ_{\text{k}} }}{{p_{\text{k}} T_{\text{sc}} }}\frac{{\left( {1 - \phi_{\text{k}} } \right)}}{{\phi_{\text{k}} }}\frac{{V_{\text{L}} p_{\text{L}} }}{{2\sqrt {Z_{\text{k}} } }}\frac{{\left( {p_{\text{k}} c_{\text{kg}} + 1} \right)}}{{\left( {\frac{{p_{\text{k}} }}{{\sqrt {Z_{\text{k}} } }} + p_{\text{L}} } \right)^{2} }}.$$
(33)

Converting Eq. 31 into the dimensionless form yields

$$\frac{1}{{r_{\text{D}}^{2} }}\frac{\partial }{{\partial r_{\text{D}} }}\left( {r_{\text{D}}^{2} \frac{{\partial p_{\text{pkD}} }}{{\partial r_{\text{D}} }}} \right) = \frac{1}{{\eta_{\text{kD}} }}\frac{{\partial p_{\text{pkD}} }}{{\partial t_{\text{D}} }}.$$
(34)

Assuming that the flow direction in secondary fractures is in the x-direction, we have the following initial and boundary conditions for Eq. 34

$$p_{\text{pkD}} \left( {r_{\text{D}} ,x_{\text{D}} ,t_{\text{D}} = 0} \right) = 0,$$
(35)
$$p_{\text{pkD}} \left( {r_{\text{D}} = 0,x_{\text{D}} ,t_{\text{D}} } \right) = p_{\text{pkD}} \left( {0,x_{\text{D}} ,t_{\text{D}} } \right)\;\left( {{\text{a}}\,{\text{finite}}\,{\text{value}}} \right),$$
(36)
$$p_{\text{pkD}} \left( {r_{\text{D}} = r_{\text{kD}} ,x_{\text{D}} ,t_{\text{D}} } \right) = p_{\text{pmD}} \left( {r_{\text{kD}} ,x_{\text{D}} ,t_{\text{D}} } \right).$$
(37)

In Eq. 36, the pressure in the kerogen core center is definitely a finite value. By defining \(F_{\text{kD}} \left( {r_{\text{D}} ,x_{\text{D}} ,t_{\text{D}} } \right) = r_{\text{D}} p_{\text{pkD}} \left( {r_{\text{D}} ,x_{\text{D}} ,t_{\text{D}} } \right)\) and conducting Laplace transform (Ozkan et al. 2010), Eqs. 34 and 36 to 37 become

$$\frac{{\partial^{2} \bar{F}_{\text{kD}} }}{{\partial r_{\text{D}}^{2} }} - \frac{s}{{\eta_{\text{kD}} }}\bar{F}_{\text{kD}} = 0\;(0 \le r_{\text{D}} \le r_{\text{kD}} ),$$
(38)
$$\bar{F}_{\text{kD}} \left( {r_{\text{D}} = 0,x_{\text{D}} ,s} \right) = 0,$$
(39)
$$F_{\text{kD}} \left( {r_{\text{kD}} ,x_{\text{D}} ,s} \right) = r_{\text{kD}} \bar{p}_{\text{pmD}} \left( {r_{\text{kD}} ,x_{\text{D}} ,s} \right).$$
(40)

Solving Eqs. 38–40 yields

$$\bar{F}_{\text{kD}} \left( {r_{\text{D}} ,x_{\text{D}} ,s} \right) = \frac{{r_{\text{kD}} { \sinh }\left( {\sqrt {s/\eta_{\text{kD}} } r_{\text{D}} } \right)}}{{{ \sinh }\left( {\sqrt {s/\eta_{\text{kD}} } r_{\text{kD}} } \right)}}\bar{p}_{\text{pmD}} \left( {r_{\text{kD}} ,x_{\text{D}} ,s} \right).$$
(41)

Therefore, the dimensionless pseudopressure is obtained

$$\bar{p}_{\text{pkD}} \left( {r_{\text{D}} ,x_{\text{D}} ,s} \right) = \frac{{r_{\text{kD}} { \sinh }\left( {\sqrt {s/\eta_{\text{kD}} } r_{\text{D}} } \right)}}{{r_{\text{D}} { \sinh }\left( {\sqrt {s/\eta_{\text{kD}} } r_{\text{kD}} } \right)}}\bar{p}_{\text{pmD}} \left( {r_{\text{kD}} ,x_{\text{D}} ,s} \right)$$
(42)

Equation 42 is the solution of the kerogen core.

3.3 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

$$\frac{1}{{r^{2} }}\frac{\partial }{\partial r}\left( {r^{2} k_{\text{appmi}} \frac{{\partial p_{\text{pm}} }}{\partial r}} \right) = \mu_{\text{mg}} c_{\text{mg}} \phi_{\text{m}} \left( {\frac{{k_{\text{appmi}} }}{{k_{\text{appm}} }}} \right)\frac{{\partial p_{\text{pm}} }}{\partial t},$$
(43)

where

$$k_{\text{appm}} = k_{\text{mfg}} = \xi_{1} \frac{{\phi_{\text{m}} }}{{\tau_{\text{m}} }}\frac{{r_{\text{pm}}^{2} }}{8}\left( {1 + \alpha_{\text{m}} Kn_{\text{m}} } \right)\left( {1 + \frac{{4Kn_{\text{m}} }}{{1 - bKn_{\text{m}} }}} \right) = \xi_{1} \frac{{\phi_{\text{m}} }}{{\tau_{\text{m}} }}\frac{{r_{\text{pm}}^{2} }}{8}\left( {1 + \alpha_{\text{m}} \frac{{\lambda_{\text{m}} }}{{r_{\text{pm}} }}} \right)\left( {1 + \frac{{4\frac{{\lambda_{\text{m}} }}{{r_{\text{pm}} }}}}{{1 - b\frac{{\lambda_{\text{m}} }}{{r_{\text{pm}} }}}}} \right),$$
(44)

Converting Eq. 43 into the dimensionless form, we have

$$\frac{1}{{r_{\text{D}}^{2} }}\frac{\partial }{{\partial r_{\text{D}} }}\left( {r_{\text{D}}^{2} \frac{{\partial p_{\text{pmD}} }}{{\partial r_{\text{D}} }}} \right) = \frac{1}{{\eta_{\text{mD}} }}\frac{{\partial p_{\text{pmD}} }}{{\partial t_{\text{D}} }}.$$
(45)

Assuming the flow direction in secondary fractures is in the x-direction, we have the following initial and boundary conditions for Eq. 45

$$p_{\text{pmD}} \left( {r_{\text{D}} ,x_{\text{D}} ,t_{\text{D}} = 0} \right) = 0,$$
(46)
$$\left. {k_{\text{appmi}} \frac{{\partial p_{\text{pmD}} }}{{\partial r_{\text{D}} }}} \right|_{{r_{\text{kD}} }} = \left. {k_{\text{appki}} \frac{{\partial p_{\text{pkD}} }}{{\partial r_{\text{D}} }}} \right|_{{r_{\text{kD}} }} ,$$
(47)
$$p_{\text{pmD}} \left( {r_{\text{D}} = r_{\text{mD}} ,x_{\text{D}} ,t_{\text{D}} } \right) = p_{\text{pfD}} \left( {r_{\text{mD}} ,x_{\text{D}} ,t_{\text{D}} } \right).$$
(48)

By defining \(F_{\text{mD}} \left( {r_{\text{D}} ,x_{\text{D}} ,t_{\text{D}} } \right) = r_{\text{D}} p_{\text{pmD}} \left( {r_{\text{D}} ,x_{\text{D}} ,t_{\text{D}} } \right)\) and conducting Laplace transform, Eqs. 45 and 47 to 48 become

$$\frac{{\partial^{2} \bar{F}_{\text{mD}} }}{{\partial r_{\text{D}}^{2} }} - \frac{s}{{\eta_{\text{mD}} }}\bar{F}_{\text{mD}} = 0,\;(r_{\text{kD}} \le r_{\text{D}} \le r_{\text{mD}} ),$$
(49)
$$\left. {\frac{{\partial \bar{F}_{\text{mD}} }}{{\partial r_{\text{D}} }}} \right|_{{r_{\text{kD}} }} = \bar{F}_{\text{mD}} \left( {r_{\text{kD}} ,x_{\text{D}} ,s} \right)\gamma ,$$
(50)
$$\bar{F}_{\text{mD}} \left( {r_{\text{mD}} ,x_{\text{D}} ,s} \right) = r_{\text{mD}} \bar{p}_{\text{pfD}} \left( {r_{\text{mD}} ,x_{\text{D}} ,s} \right),$$
(51)

where

$$\gamma = \frac{1}{{r_{\text{kD}} }}\left\{ {1 + \frac{{k_{\text{appki}} }}{{k_{\text{appmi}} }}\left[ {r_{\text{kD}} \sqrt {\frac{s}{{\eta_{\text{kD}} }}} { \coth }\left( {\sqrt {s/\eta_{\text{kD}} } r_{\text{kD}} } \right) - 1} \right]} \right\}.$$
(52)

Solving Eqs. 49–52 yields

$$\begin{aligned} \bar{F}_{\text{mD}} \left( {r_{\text{D}} ,x_{\text{D}} ,s} \right) \hfill \\ = \frac{{r_{\text{mD}} \left\{ {\left( {\frac{{\sqrt {\frac{s}{{\eta_{\text{mD}} }}} + \gamma }}{{\sqrt {\frac{s}{{\eta_{\text{mD}} }}} - \gamma }}} \right){ \exp }\left[ {\sqrt {\frac{s}{{\eta_{\text{mD}} }}} \left( {r_{\text{D}} - 2r_{\text{kD}} } \right)} \right] + { \exp }\left( { - \sqrt {\frac{s}{{\eta_{\text{mD}} }}} r_{\text{D}} } \right)} \right\}}}{{\left( {\frac{{\sqrt {\frac{s}{{\eta_{\text{mD}} }}} + \gamma }}{{\sqrt {\frac{s}{{\eta_{\text{mD}} }}} - \gamma }}} \right){ \exp }\left[ {\sqrt {\frac{s}{{\eta_{\text{mD}} }}} \left( {r_{\text{mD}} - 2r_{\text{kD}} } \right)} \right] + { \exp }\left( { - \sqrt {\frac{s}{{\eta_{\text{mD}} }}} r_{\text{mD}} } \right)}}\bar{p}_{\text{pfD}} \left( {r_{\text{mD}} ,x_{\text{D}} ,s} \right). \hfill \\ \end{aligned}$$
(53)

Therefore, the dimensionless pseudopressure is obtained

$$\begin{aligned} \bar{p}_{\text{pmD}} \left( {r_{\text{D}} ,x_{\text{D}} ,s} \right) \hfill \\ = \frac{{r_{\text{mD}} \left\{ {\left( {\frac{{\sqrt {\frac{s}{{\eta_{\text{mD}} }}} + \gamma }}{{\sqrt {\frac{s}{{\eta_{\text{mD}} }}} - \gamma }}} \right){ \exp }\left[ {\sqrt {\frac{s}{{\eta_{\text{mD}} }}} \left( {r_{\text{D}} - 2r_{\text{kD}} } \right)} \right] + { \exp }\left( { - \sqrt {\frac{s}{{\eta_{\text{mD}} }}} r_{\text{D}} } \right)} \right\}}}{{r_{\text{D}} \left\{ {\left( {\frac{{\sqrt {\frac{s}{{\eta_{\text{mD}} }}} + \gamma }}{{\sqrt {\frac{s}{{\eta_{\text{mD}} }}} - \gamma }}} \right){ \exp }\left[ {\sqrt {\frac{s}{{\eta_{\text{mD}} }}} \left( {r_{\text{mD}} - 2r_{\text{kD}} } \right)} \right] + { \exp }\left( { - \sqrt {\frac{s}{{\eta_{\text{mD}} }}} r_{\text{mD}} } \right)} \right\}}}\bar{p}_{\text{pfD}} \left( {r_{\text{mD}} ,x_{\text{D}} ,s} \right). \hfill \\ \end{aligned}$$
(54)

Equation 54 is the solution of the exterior inorganic sphere.

3.4 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)

$$f\left( {x,t} \right) = - \left. {\frac{2}{{h_{\text{f}} }}\left( {\rho_{\text{mg}} \frac{{k_{\text{appm}} }}{{\mu_{\text{mg}} }}\frac{{\partial p_{\text{m}} }}{\partial r}} \right)} \right|_{{\left( {r_{\text{m}} ,x,t} \right)}} .$$
(55)

Therefore, we have the mass balance equation for secondary fractures

$$- \frac{{\partial \left( {\rho_{\text{fg}} v_{\text{f}} } \right)}}{\partial x} + \left[ { - \left. {\frac{2}{{h_{\text{f}} }}\left( {\rho_{\text{mg}} \frac{{k_{\text{appm}} }}{{\mu_{\text{mg}} }}\frac{{\partial p_{\text{m}} }}{\partial r}} \right)} \right|_{{\left( {r_{\text{m}} ,x,t} \right)}} } \right] = \frac{{\partial \left( {\rho_{\text{fg}} \phi_{\text{f}} } \right)}}{\partial t}.$$
(56)

The dimensionless pseudopressure form of Eq. 56 is

$$\frac{{\partial^{2} p_{\text{pfD}} }}{{\partial x_{\text{D}}^{2} }} - \left. {\frac{{2k_{\text{appmi}} d_{\text{ref}} }}{{h_{\text{f}} k_{\text{f}} }}\frac{{\partial p_{\text{pmD}} }}{{\partial r_{\text{D}} }}} \right|_{{\left( {r_{\text{mD}} ,x_{\text{D}} ,t_{\text{D}} } \right)}} = \frac{1}{{\eta_{\text{fD}} }}\frac{{\partial p_{\text{pfD}} }}{{\partial t_{\text{D}} }}.$$
(57)

Converting Eq. 57 into the Laplace form yields

$$\frac{{\partial^{2} \bar{p}_{\text{pfD}} }}{{\partial x_{\text{D}}^{2} }} - c\left( s \right)\bar{p}_{\text{pfD}} = 0,$$
(58)

where

$$c\left( s \right) = \frac{s}{{\eta_{\text{fD}} }} + \frac{{2k_{\text{appmi}} d_{\text{ref}} }}{{h_{\text{f}} k_{\text{f}} }}a\left( s \right),$$
(59)
$$\begin{aligned} a\left( s \right) & = - \frac{1}{{r_{\text{mD}} }} \\ & \quad + \frac{{\sqrt {\frac{s}{{\eta_{\text{mD}} }}} \left\{ {\left( {\frac{{\sqrt {\frac{s}{{\eta_{\text{mD}} }}} + \gamma }}{{\sqrt {\frac{s}{{\eta_{\text{mD}} }}} - \gamma }}} \right){ \exp }\left[ {\sqrt {\frac{s}{{\eta_{\text{mD}} }}} \left( {r_{\text{mD}} - 2r_{\text{kD}} } \right)} \right] - { \exp }\left( { - \sqrt {\frac{s}{{\eta_{\text{mD}} }}} r_{\text{mD}} } \right)} \right\}}}{{\left( {\frac{{\sqrt {\frac{s}{{\eta_{\text{mD}} }}} + \gamma }}{{\sqrt {\frac{s}{{\eta_{\text{mD}} }}} - \gamma }}} \right){ \exp }\left[ {\sqrt {\frac{s}{{\eta_{\text{mD}} }}} \left( {r_{\text{mD}} - 2r_{\text{kD}} } \right)} \right] + { \exp }\left( { - \sqrt {\frac{s}{{\eta_{\text{mD}} }}} r_{\text{mD}} } \right)}}. \\ \end{aligned}$$
(60)

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.

3.4.1 Region 6

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

$$\frac{{\partial^{2} \bar{p}_{{{\text{pfD}}6}} }}{{\partial z_{\text{D}}^{2} }} - c\left( s \right)_{6} \bar{p}_{{{\text{pfD}}6}} = 0,$$
(61)
$$\left. {\frac{{\partial \bar{p}_{{{\text{pfD}}6}} }}{{\partial z_{\text{D}} }}} \right|_{{z_{\text{D}} = z_{{2{\text{D}}}} }} = 0,$$
(62)
$$\left. {\bar{p}_{{{\text{pfD}}6}} } \right|_{{z_{\text{D}} = z_{{1{\text{D}}}} }} = \left. {\bar{p}_{{{\text{pfD}}2}} } \right|_{{z_{\text{D}} = z_{{1{\text{D}}}} }} = \left. {\bar{p}_{{{\text{pfD}}4}} } \right|_{{z_{\text{D}} = z_{{1{\text{D}}}} }} .$$
(63)

Solving Eqs. 61–63 yields

$$\begin{aligned} \left. {\frac{{\partial \bar{p}_{{{\text{pfD}}6}} }}{{\partial z_{\text{D}} }}} \right|_{{z_{\text{D}} = z_{{1{\text{D}}}} }} & = - \left. {\bar{p}_{{{\text{pfD}}2}} } \right|_{{z_{\text{D}} = z_{{1{\text{D}}}} }} \sqrt {c\left( s \right)_{6} } { \tanh }\left[ {\sqrt {c\left( s \right)_{6} } \left( {z_{{2{\text{D}}}} - z_{{1{\text{D}}}} } \right)} \right] \\ & = - \left. {\bar{p}_{{{\text{pfD}}4}} } \right|_{{z_{\text{D}} = z_{{1{\text{D}}}} }} \sqrt {c\left( s \right)_{6} } { \tanh }\left[ {\sqrt {c\left( s \right)_{6} } \left( {z_{{2{\text{D}}}} - z_{{1{\text{D}}}} } \right)} \right]. \\ \end{aligned}$$
(64)

Equation 64 is the solution of region 6 and will be used in the diffusivity equations of regions 2 and 4.

3.4.2 Region 5

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

$$\frac{{\partial^{2} \bar{p}_{{{\text{pfD}}5}} }}{{\partial z_{\text{D}}^{2} }} - c\left( s \right)_{5} \bar{p}_{{{\text{pfD}}5}} = 0,$$
(65)
$$\left. {\frac{{\partial \bar{p}_{{{\text{pfD}}5}} }}{{\partial z_{\text{D}} }}} \right|_{{z_{\text{D}} = z_{{2{\text{D}}}} }} = 0,$$
(66)
$$\left. {\bar{p}_{{{\text{pfD}}5}} } \right|_{{z_{\text{D}} = z_{{1{\text{D}}}} }} = \left. {\bar{p}_{{{\text{pfD}}1}} } \right|_{{z_{\text{D}} = z_{{1{\text{D}}}} }} = \left. {\bar{p}_{{{\text{pfD}}3}} } \right|_{{z_{\text{D}} = z_{{1{\text{D}}}} }} .$$
(67)

Solving Eqs. 65–67 yields

$$\begin{aligned} \left. {\frac{{\partial \bar{p}_{{{\text{pfD}}5}} }}{{\partial z_{\text{D}} }}} \right|_{{z_{\text{D}} = z_{{1{\text{D}}}} }} & = - \left. {\bar{p}_{{{\text{pfD}}1}} } \right|_{{z_{\text{D}} = z_{{1{\text{D}}}} }} \sqrt {c\left( s \right)_{5} } { \tanh }\left[ {\sqrt {c\left( s \right)_{5} } \left( {z_{{2{\text{D}}}} - z_{{1{\text{D}}}} } \right)} \right] \\ & = - \left. {\bar{p}_{{{\text{pfD}}3}} } \right|_{{z_{\text{D}} = z_{{1{\text{D}}}} }} \sqrt {c\left( s \right)_{5} } { \tanh }\left[ {\sqrt {c\left( s \right)_{5} } \left( {z_{{2{\text{D}}}} - z_{{1{\text{D}}}} } \right)} \right]. \\ \end{aligned}$$
(68)

Equation 68 is the solution of region 5 and will be used in the diffusivity equations of regions 1 and 3.

3.4.3 Region 4

Region 4 involves 1-D flow in the x-direction. The diffusivity equation, no-flow boundary and pressure continuity conditions are

$$\frac{{\partial^{2} \bar{p}_{{{\text{pfD}}4}} }}{{\partial x_{\text{D}}^{2} }} + \frac{{k_{{{\text{f}}6}} }}{{k_{{{\text{f}}4}} z_{{1{\text{D}}}} }}\left. {\frac{{\partial \bar{p}_{{{\text{pfD}}6}} }}{{\partial z_{\text{D}} }}} \right|_{{z_{\text{D}} = z_{{1{\text{D}}}} }} - c\left( s \right)_{4} \bar{p}_{{{\text{pfD}}4}} = 0,$$
(69)
$$\left. {\frac{{\partial \bar{p}_{{{\text{pfD}}4}} }}{{\partial x_{\text{D}} }}} \right|_{{x_{\text{D}} = x_{\text{eD}} }} = 0,$$
(70)
$$\left. {\bar{p}_{{{\text{pfD}}4}} } \right|_{{x_{\text{D}} = x_{{1{\text{D}}}} }} = \left. {\bar{p}_{{{\text{pfD}}2}} } \right|_{{x_{\text{D}} = x_{{1{\text{D}}}} }} .$$
(71)

Solving Eqs. 69–71 yields

$$\left. {\frac{{\partial \bar{p}_{{{\text{pfD}}4}} }}{{\partial x_{\text{D}} }}} \right|_{{x_{\text{D}} = x_{{1{\text{D}}}} }} = - \left. {\bar{p}_{{{\text{pfD}}2}} } \right|_{{x_{\text{D}} = x_{{1{\text{D}}}} }} \sqrt {\alpha_{4} } { \tanh }\left[ {\sqrt {\alpha_{4} } \left( {x_{\text{eD}} - x_{{1{\text{D}}}} } \right)} \right],$$
(72)

where

$$\alpha_{4} = \frac{{k_{{{\text{f}}6}} }}{{k_{{{\text{f}}4}} z_{{1{\text{D}}}} }}\sqrt {c\left( s \right)_{6} } { \tanh }\left[ {\sqrt {c\left( s \right)_{6} } \left( {z_{{2{\text{D}}}} - z_{{1{\text{D}}}} } \right)} \right] + c\left( s \right)_{4} .$$
(73)

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

3.4.4 Region 3

Region 3 involves 1-D flow in the x-direction. The diffusivity equation, no-flow boundary and pressure continuity conditions are

$$\frac{{\partial^{2} \bar{p}_{{{\text{pfD}}3}} }}{{\partial x_{\text{D}}^{2} }} + \frac{{k_{{{\text{f}}5}} }}{{k_{{{\text{f}}3}} z_{{1{\text{D}}}} }}\left. {\frac{{\partial \bar{p}_{{{\text{pfD}}5}} }}{{\partial z_{\text{D}} }}} \right|_{{z_{\text{D}} = z_{{1{\text{D}}}} }} - c\left( s \right)_{3} \bar{p}_{{{\text{pfD}}3}} = 0,$$
(74)
$$\left. {\frac{{\partial \bar{p}_{{{\text{pfD}}3}} }}{{\partial x_{\text{D}} }}} \right|_{{x_{\text{D}} = x_{\text{eD}} }} = 0,$$
(75)
$$\left. {\bar{p}_{{{\text{pfD}}3}} } \right|_{{x_{\text{D}} = x_{{1{\text{D}}}} }} = \left. {\bar{p}_{{{\text{pfD}}1}} } \right|_{{x_{\text{D}} = x_{{1{\text{D}}}} }} .$$
(76)

Solving Eqs. 74–76 yields

$$\left. {\frac{{\partial \bar{p}_{{{\text{pfD}}3}} }}{{\partial x_{\text{D}} }}} \right|_{{x_{\text{D}} = x_{{1{\text{D}}}} }} = - \left. {\bar{p}_{{{\text{pfD}}1}} } \right|_{{x_{\text{D}} = x_{{1{\text{D}}}} }} \sqrt {\alpha_{3} } { \tanh }\left[ {\sqrt {\alpha_{3} } \left( {x_{\text{eD}} - x_{{1{\text{D}}}} } \right)} \right],$$
(77)

where

$$\alpha_{3} = \frac{{k_{{{\text{f}}5}} }}{{k_{{{\text{f}}3}} z_{{1{\text{D}}}} }}\sqrt {c\left( s \right)_{5} } { \tanh }\left[ {\sqrt {c\left( s \right)_{5} } \left( {z_{{2{\text{D}}}} - z_{{1{\text{D}}}} } \right)} \right] + c\left( s \right)_{3} .$$
(78)

Equation 77 is the solution of region 3 and will be used in the diffusivity equation of region 1.

3.4.5 Region 2

Region 2 involves 1-D flow in the y-direction. The diffusivity equation, no-flow boundary and pressure continuity conditions are

$$\frac{{\partial^{2} \bar{p}_{{{\text{pfD}}2}} }}{{\partial y_{\text{D}}^{2} }} + \frac{{k_{{{\text{f}}6}} }}{{k_{{{\text{f}}2}} z_{{1{\text{D}}}} }}\left. {\frac{{\partial \bar{p}_{{{\text{pfD}}6}} }}{{\partial z_{\text{D}} }}} \right|_{{z_{\text{D}} = z_{{1{\text{D}}}} }} + \frac{{k_{{{\text{f}}4}} }}{{k_{{{\text{f}}2}} x_{{1{\text{D}}}} }}\left. {\frac{{\partial \bar{p}_{{{\text{pfD}}4}} }}{{\partial x_{\text{D}} }}} \right|_{{x_{\text{D}} = x_{{1{\text{D}}}} }} - c\left( s \right)_{2} \bar{p}_{{{\text{pfD}}2}} = 0,$$
(79)
$$\left. {\frac{{\partial \bar{p}_{{{\text{pfD}}2}} }}{{\partial y_{\text{D}} }}} \right|_{{y_{\text{D}} = y_{{2{\text{D}}}} }} = 0,$$
(80)
$$\left. {\bar{p}_{{{\text{pfD}}2}} } \right|_{{y_{\text{D}} = y_{{1{\text{D}}}} }} = \left. {\bar{p}_{{{\text{pfD}}1}} } \right|_{{y_{\text{D}} = y_{{1{\text{D}}}} }} .$$
(81)

Solving Eqs. 79–81 yields

$$\left. {\frac{{\partial \bar{p}_{{{\text{pfD}}2}} }}{{\partial y_{\text{D}} }}} \right|_{{y_{\text{D}} = y_{{1{\text{D}}}} }} = - \left. {\bar{p}_{{{\text{pfD}}1}} } \right|_{{y_{\text{D}} = y_{{1{\text{D}}}} }} \sqrt {\alpha_{2} } { \tanh }\left[ {\sqrt {\alpha_{2} } \left( {y_{{2{\text{D}}}} - y_{{1{\text{D}}}} } \right)} \right],$$
(82)

where

$$\alpha_{2} = \frac{{k_{{{\text{f}}6}} }}{{k_{{{\text{f}}2}} z_{{1{\text{D}}}} }}\sqrt {c\left( s \right)_{6} } { \tanh }\left[ {\sqrt {c\left( s \right)_{6} } \left( {z_{{2{\text{D}}}} - z_{{1{\text{D}}}} } \right)} \right] + \frac{{k_{{{\text{f}}4}} }}{{k_{{{\text{f}}2}} x_{{1{\text{D}}}} }}\sqrt {\alpha_{4} } { \tanh }\left[ {\sqrt {\alpha_{4} } \left( {x_{\text{eD}} - x_{{1{\text{D}}}} } \right)} \right] + c\left( s \right)_{2} .$$
(83)

Equation 83 is the solution of region 2 and will be used in the diffusivity equation of region 1.

3.4.6 Region 1

Region 1 involves 1-D flow in the y-direction. The diffusivity equation, flux and pressure continuity conditions are

$$\frac{{\partial^{2} \bar{p}_{{{\text{pfD}}1}} }}{{\partial y_{\text{D}}^{2} }} + \frac{{k_{{{\text{f}}5}} }}{{k_{{{\text{f}}1}} z_{{1{\text{D}}}} }}\left. {\frac{{\partial \bar{p}_{{{\text{pfD}}5}} }}{{\partial z_{\text{D}} }}} \right|_{{z_{\text{D}} = z_{{1{\text{D}}}} }} + \frac{{k_{{{\text{f}}3}} }}{{k_{{{\text{f}}1}} x_{{1{\text{D}}}} }}\left. {\frac{{\partial \bar{p}_{{{\text{pfD}}3}} }}{{\partial x_{\text{D}} }}} \right|_{{x_{\text{D}} = x_{{1{\text{D}}}} }} - c\left( s \right)_{1} \bar{p}_{{{\text{pfD}}1}} = 0,$$
(84)
$$\left. {k_{{{\text{f}}1}} \frac{{\partial \bar{p}_{{{\text{pfD}}1}} }}{{\partial y_{\text{D}} }}} \right|_{{y_{\text{D}} = y_{{1{\text{D}}}} }} = \left. {k_{{{\text{f}}2}} \frac{{\partial \bar{p}_{{{\text{pfD}}2}} }}{{\partial y_{\text{D}} }}} \right|_{{y_{\text{D}} = y_{{1{\text{D}}}} }} ,$$
(85)
$$\left. {\bar{p}_{{{\text{pfD}}1}} } \right|_{{y_{\text{D}} = w_{\text{D}} /2}} = \left. {\bar{p}_{\text{pFD}} } \right|_{{y_{\text{D}} = w_{\text{D}} /2}} .$$
(86)

Solving Eqs. 84–86 yields

$$\left. {\frac{{\partial \bar{p}_{{{\text{pfD}}1}} }}{{\partial y_{\text{D}} }}} \right|_{{y_{\text{D}} = w_{\text{D}} /2}} = - \left. {\bar{p}_{\text{pFD}} } \right|_{{y_{\text{D}} = w_{\text{D}} /2}} \beta_{1} ,$$
(87)

where

$$\beta_{1} = \sqrt {\alpha_{1} } \frac{{{ \exp }\left( { - \sqrt {\alpha_{1} } \frac{{w_{\text{D}} }}{2}} \right) - { \exp }\left( { - 2\sqrt {\alpha_{1} } y_{{1{\text{D}}}} } \right){ \exp }\left( {\sqrt {\alpha_{1} } \frac{{w_{\text{D}} }}{2}} \right)\frac{{k_{{{\text{f}}1}} \sqrt {\alpha_{1} } - k_{{{\text{f}}2}} \sqrt {\alpha_{2} } { \tanh }\left[ {\sqrt {\alpha_{2} } \left( {y_{{2{\text{D}}}} - y_{{1{\text{D}}}} } \right)} \right]}}{{k_{{{\text{f}}1}} \sqrt {\alpha_{1} } + k_{{{\text{f}}2}} \sqrt {\alpha_{2} } { \tanh }\left[ {\sqrt {\alpha_{2} } \left( {y_{{2{\text{D}}}} - y_{{1{\text{D}}}} } \right)} \right]}}}}{{{ \exp }\left( { - \sqrt {\alpha_{1} } \frac{{w_{\text{D}} }}{2}} \right) + { \exp }\left( { - 2\sqrt {\alpha_{1} } y_{{1{\text{D}}}} } \right){ \exp }\left( {\sqrt {\alpha_{1} } \frac{{w_{\text{D}} }}{2}} \right)\frac{{k_{{{\text{f}}1}} \sqrt {\alpha_{1} } - k_{{{\text{f}}2}} \sqrt {\alpha_{2} } { \tanh }\left[ {\sqrt {\alpha_{2} } \left( {y_{{2{\text{D}}}} - y_{{1{\text{D}}}} } \right)} \right]}}{{k_{{{\text{f}}1}} \sqrt {\alpha_{1} } + k_{{{\text{f}}2}} \sqrt {\alpha_{2} } { \tanh }\left[ {\sqrt {\alpha_{2} } \left( {y_{{2{\text{D}}}} - y_{{1{\text{D}}}} } \right)} \right]}}}},$$
(88)
$$\alpha_{1} = \frac{{k_{{{\text{f}}3}} }}{{k_{{{\text{f}}1}} x_{{1{\text{D}}}} }}\sqrt {\alpha_{3} } { \tanh }\left[ {\sqrt {\alpha_{3} } \left( {x_{\text{eD}} - x_{{1{\text{D}}}} } \right)} \right] + \frac{{k_{{{\text{f}}5}} }}{{k_{{{\text{f}}1}} z_{{1{\text{D}}}} }}\sqrt {c\left( s \right)_{5} } { \tanh }\left[ {\sqrt {c\left( s \right)_{5} } \left( {z_{{2{\text{D}}}} - z_{{1{\text{D}}}} } \right)} \right] + c\left( s \right)_{1} .$$
(89)

Equation 87 is the solution of region 1 and will be used in the diffusivity equation of the hydraulic fracture region.

3.4.7 Hydraulic fracture region

The hydraulic fracture region involves 1-D flow in the x-direction. The diffusivity equation and boundary conditions are

$$\frac{{\partial^{2} \bar{p}_{\text{pFD}} }}{{\partial x_{\text{D}}^{2} }} + \frac{{k_{{{\text{f}}1{\text{b}}}} }}{{k_{\text{F}} \left( {w_{\text{D}} /2} \right)}}\left. {\frac{{\partial \bar{p}_{{{\text{pfD}}1}} }}{{\partial y_{\text{D}} }}} \right|_{{y_{\text{D}} = w_{\text{D}} /2}} - \frac{s}{{\eta_{\text{FD}} }}\bar{p}_{\text{pFD}} = 0,$$
(90)
$$\left. {\frac{{\partial \bar{p}_{\text{pFD}} }}{{\partial x_{\text{D}} }}} \right|_{{x_{\text{D}} = x_{{1{\text{D}}}} }} = 0,$$
(91)
$$\left. {\frac{{\partial \bar{p}_{\text{pFD}} }}{{\partial x_{\text{D}} }}} \right|_{{x_{\text{D}} = 0}} = - \frac{{\bar{q}_{{{\text{F}}j}} }}{{q_{\text{Ft}} }}\frac{{\pi k_{\text{ref}} h_{\text{ref}} d_{\text{ref}} }}{{k_{\text{F}} h_{\text{F}} w_{\text{F}} }}.$$
(92)

Here the subscript j represents the j-th hydraulic fracture. In Eq. 90, the bulk secondary fracture permeability of region 1 \((k_{{{\text{f}}1{\text{b}}}} )\) is used, while in other regions we equivalently use the intrinsic secondary fracture permeability as the fracture-matrix 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),

$$k_{\text{fb}} = k_{\text{f}} \frac{{V_{\text{f}} }}{{V_{\text{f}} + V_{\text{m}} }} = k_{\text{f}} \frac{{4\pi r_{\text{m}}^{2} h_{\text{f}} /2}}{{4\pi r_{\text{m}}^{2} h_{\text{f}} /2 + 4\pi r_{\text{m}}^{3} /3}} = \frac{{3k_{\text{f}} h_{\text{f}} }}{{3h_{\text{f}} + 2r_{\text{m}} }},$$
(93)

where \(k_{\text{f}}\) is the secondary fracture intrinsic permeability in mD; \(h_{\text{f}}\) is the secondary fracture layer thickness in ft; and \(r_{\text{m}}\) is the inorganic matrix radius in ft. Solving Eqs. 90–92 yields

$$\left. {\bar{p}_{{{\text{pFD}}j}} } \right|_{{x_{\text{D}} = 0}} = \frac{{\bar{q}_{{{\text{F}}j}} }}{{q_{\text{Ft}} }}\frac{{\pi k_{\text{ref}} h_{\text{ref}} d_{\text{ref}} }}{{k_{\text{F}} h_{\text{F}} w_{\text{F}} \sqrt {\alpha_{\text{F}} } { \tanh }\left( {\sqrt {\alpha_{\text{F}} } x_{{1{\text{D}}}} } \right)}},$$
(94)

where

$$\alpha_{\text{F}} = \frac{{k_{{{\text{f}}1{\text{b}}}} }}{{k_{\text{F}} \left( {w_{\text{D}} /2} \right)}}\beta_{1} + \frac{s}{{\eta_{\text{FD}} }}.$$
(95)

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

$$p_{{{\text{pD}}j}} = \frac{{\left( {p_{\text{pi}} - p_{\text{pwf}} } \right)}}{{\left( {p_{\text{pi}} - p_{\text{pwf}} } \right)_{\text{ref}} }},$$
(96)

where \(p_{\text{wf}}\) is the bottom-hole flowing pressure in psi. By assuming \(p_{\text{i}} - p_{\text{wf}} = \left( {p_{\text{i}} - p_{\text{wf}} } \right)_{\text{ref}}\) and ignoring the pressure drop in the horizontal wellbore for gas reservoirs, we have

$$\left. {\bar{p}_{{{\text{pFD}}j}} } \right|_{{x_{\text{D}} = 0}} = \frac{1}{s}.$$
(97)

Then the final solution of constant pressure cases is

$$\bar{q}_{{{\text{D}}j}} = \frac{{\bar{q}_{{{\text{F}}j}} }}{{q_{\text{Ft}} }}\frac{1}{{\left. {\bar{p}_{{{\text{pFD}}j}} } \right|_{{x_{\text{D}} = 0}} }}\frac{1}{s},$$
(98)

Note that, in Eq. 98, \(\bar{p}_{{{\text{pFD}}j}}\) is obtained from Eq. 94. And the total rate of an MFHW is expressed by

$$\bar{q}_{\text{D}} = \mathop \sum \limits_{j = 1}^{n} \bar{q}_{{{\text{D}}j}} .$$
(99)

Consequently, the cumulative production is given by

$$\bar{Q}_{\text{D}} = \frac{{\bar{q}_{\text{D}} }}{s}.$$
(100)

According to Callard and Schenewerk (1995) and Chaudhry (2003), under constant pressure production conditions, the dimensionless rate and cumulative production are defined as follows:

$$q_{\text{D}} = \frac{1422qT}{{k_{\text{ref}} h_{\text{ref}} \left( {p_{\text{pi}} - p_{\text{pwf}} } \right)_{\text{ref}} }},$$
(101)
$$Q_{\text{D}} = \frac{9.009QT}{{\phi_{\text{ref}} c_{\text{tref}} h_{\text{ref}}\mu_{\text{ref}} d_{\text{ref}}^{2} \left( {p_{\text{pi}} - p_{\text{pwf}} } \right)_{\text{ref}} }}.$$
(102)

Here \(Q\) is the cumulative production in Mscf.

3.5 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_{\text{m}} = r_{\text{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. (2017, 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. 2017, 2018, 2020): (1) the width of region 1 is no smaller than 10% of the half fracture spacing \((y_{1} \ge 0.1y_{2} )\); (2) the fracture length is at least 10% of the reservoir width \((x_{\text{F}} \ge 0.1x_{\text{e}} )\); (3) at least 60% of the formation is vertically penetrated by hydraulic fractures \((h_{\text{F}} \ge 0.6h)\). Note that if the fracture height is too small \((h_{\text{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.

Table 2 Input parameters used for model validation
Full size table
Fig. 3
figure 3

Comparison of this model and the published model

Full size image

3.6 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. The well is an MFHW producing under a constant bottom-hole pressure condition with 28 hydraulic fractures (Al-Ahmadi and Wattenbarger 2011; Yu and Sepehrnoori 2014). The basic reservoir and well properties, obtained from the original paper and related studies (Al-Ahmadi and Wattenbarger 2011; Brown et al. 2011; Apaydin et al. 2012; Yu and Sepehrnoori 2014), are listed in 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 \times 10^{ - 4}\) ft to \(5 \times 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): \(\tau = \phi^{ - 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.

Fig. 4
figure 4

Comparison of simulation results and Barnett Shale production data

Full size image
Table 3 Input parameters used in history matching
Full size table

4 Parametric investigations

4.1 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).

Fig. 5
figure 5

Typical flow regimes of an MFHW in a shale gas reservoir with closed outer boundaries

Full size image
Table 4 Input parameters used for the flow regime diagnosis
Full size table

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).

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 matrix-fracture 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 pseudo-steady state has been achieved.

4.2 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 (xD = 1/2x1D and yD = 1/2y1D) 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 (\(r > 3.0\) ft) where the dimensionless pseudopressure values almost overlap each other with nearly the same pressure drop, there are still pseudopressure differences at different locations within the kerogen core (\(r \le 3.0\) ft). Even the whole matrix reaches the dynamic equilibrium state, the dimensionless rate profile (\(t = 13.6\) years) shows that the kerogen core (\(r \le 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/2y1) 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.

Fig. 6
figure 6

Dimensionless pseudopressure distribution of a matrix block located at xD = 1/2x1D and yD = 1/2y1D

Full size image
Fig. 7
figure 7

Dimensionless rate profile of a matrix block located at xD = 1/2x1D and yD = 1/2y1D under a dynamic equilibrium state

Full size image

4.3 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 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. 8
figure 8

Effects of the kerogen pore radius on pressure behavior

Full size image
Fig. 9
figure 9

Effects of the inorganic matrix pore radius on pressure behavior

Full size image
Fig. 10
figure 10

Effects of the inorganic matrix pore radius on rate behavior

Full size image
Fig. 11
figure 11

Effects of the inorganic matrix pore radius on cumulative production

Full size image

4.4 Effects of matrix porosity

In this section, the effects of both kerogen and inorganic matrix porosity are analyzed. The kerogen porosity values used here are greater than or equal to those of inorganic matter because pores are better developed in organic matter (Kuchinskiy 2013). The kerogen porosity can be five times higher than that of inorganic matrices (Wang and Reed 2009). Figures 12 and 13 illustrate that the influence of kerogen porosity starts from the late matrix-fracture inter-porosity flow regime. With the kerogen porosity increases, the dimensionless pseudopressure and derivative values become smaller, and the dimensionless rates are slightly larger in these affected regimes. The matrix-fracture inter-porosity regime also ends later. However, it only slightly affects the bilinear flow regime of the pressure-derivative curve because this regime is controlled by natural fracture properties. Higher kerogen porosity also puts off the arrival of the unit-slope boundary-dominated regime because more free gas is stored in kerogen and can contribute to production even though the adsorbed gas amount is a little smaller. Consequently, this results in the increases in cumulative production as shown in Fig. 14. Although the inorganic porosity is smaller than that of kerogen, inorganic matter occupies a larger matrix volume (87.5% of the total matrix volume in the given case). Compared with kerogen porosity, the inorganic porosity overall has a heavier impact on well performance. It controls the whole production process except the early-mid time fracture linear flow regime, as shown in Figs. 15 and 16. However, its impact becomes marginal during the later matrix-fracture inter-porosity flow regime and the transient bilinear flow regime because these regimes are dominated by kerogen and natural fracture properties respectively. As the inorganic matrix occupies a larger portion of the total matrix volume, the delay in the arrival of boundary effects and the increase of cumulative production are more noticeable with the increase of inorganic porosity, as shown in Figs. 15, 16 and 17.

Fig. 12
figure 12

Effects of kerogen porosity on pressure behavior

Full size image
Fig. 13
figure 13

Effects of kerogen porosity on rate behavior

Full size image
Fig. 14
figure 14

Effects of kerogen porosity on cumulative production

Full size image
Fig. 15
figure 15

Effects of inorganic matrix porosity on pressure behavior

Full size image
Fig. 16
figure 16

Effects of inorganic matrix porosity on rate behavior

Full size image
Fig. 17
figure 17

Effects of inorganic matrix porosity on cumulative production

Full size image

4.5 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_{\text{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 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.

Fig. 18
figure 18

Effects of the Langmuir volume on pressure behavior

Full size image
Fig. 19
figure 19

Effects of the Langmuir volume on rate behavior

Full size image
Fig. 20
figure 20

Effects of the Langmuir volume on cumulative production

Full size image

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 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.

Fig. 21
figure 21

Effects of the Langmuir pressure on pressure behavior

Full size image
Fig. 22
figure 22

Effects of the Langmuir pressure on rate behavior

Full size image
Fig. 23
figure 23

Effects of the Langmuir pressure on cumulative production

Full size image

4.6 Contributions of different linear flow regions on gas production

In essence, the reservoir-fracture flow model is a seven-region 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_{\text{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.

Fig. 24
figure 24

Contributions of different linear flow regions on gas production

Full size image

5 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. (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. (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. (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. (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. (5)

    The occurrence of contributions from different linear flow regions depends on their relative location, while the cumulative contribution is dominated by their sizes.