© M.F. El-Amin et al., published by IFP Energies nouvelles, 2020

Nomenclature

b_d : Factor of slippage

k_0,d: Intrinsic permeability

${\varphi }_m$ : Porosity of the matrix blocks

${\varphi }_f$ : Porosity of the fractures

${\rho }_m$ : Mass density of the gas in the matrix blocks

${\rho }_f$ : Mass density of the gas in the fractures

${\rho }_s$ : Mass density of the core

$M_w$ : Gas molar weight

$P_L$ : Langmuir pressure

$p_m$ : Gas pressure in matrix blocks

$p_f$ : Gas pressure in fractures

$p_c$ : Critical pressure

$p_{\mathrm{fe}}$ : Average fractures pressure around the well

$V_L$ : Langmuir volume

$V_{\mathrm{std}}$ : Mole volume under standard conditions

$R$ : General constant of gases

$T$ : Temperature

$T_c$ : Critical temperature

$Z$ : Gas compressibility factor

$r_w$ : Well-radius

$r_e$ : Drainage-radius

$n$ : Set of normal fractures

${\sigma }_m$ : Mean stress in matrix blocks

${\sigma }_f$ : Mean stress in fractures

${\sigma }_m^{\mathrm{\prime }}$ : Effective stress in matrix blocks

${\sigma }_f^{\mathrm{\prime }}$ : Effective stress in fractures

${\alpha }_d$ : Biot’s effective parameter

$L_x$ : Fracture spacing in $x$

$L_y$ : Fracture spacing in $y$

$L_z$ : Fracture spacing in $z$

Superscripts and subscripts

$f$ : Fractures

$m$ : Matrix

$d$ : $m$ or $f$

1 Introduction

Finite Element Methods (FEMs) are effective numerical techniques for solving the complex engineering problems. FEMs appeared in the middle of the last century and were widely used in solid mechanical applications [1]. The FEM has been extensively investigated and developed in the last five decades [2], and now is used for highly nonlinear problems [3]. However, the standard FEM faced certain instabilities for solving elliptic problem that describe the flow in porous media [4, 5]. On the other hand, the Mixed Finite Element Methods (MFEMs) have succeeded in eliminating such instabilities [6, 7] as it may be extended to higher-order approximations as well as it is a locally conservative method [8]. Different physical variables can be treated differently in the MFEM from other variables. For example, the enhanced linear Brezzi-Douglas-Marini (BDM1) space is used for speed approximation, which requires more precision, while the piece-wise constant space is used for pressure approximation. The models of shale-gas transport were mainly developed by adapting the traditional models of flow in fractured porous media. Warren and Root [9] have developed an idealized geometric model to investigate the behavioral flow in fractured porous media. The dual porosity model including stress of fracture was used to investigate the flow of gas in the shale [10, 11], while, the dual-continuum model was used to describe gas flow Kerogen shale [12, 13]. Several versions of the dual-continuum model have been developed to include, for instance, adsorption, heat, deformation and Knudsen diffusion [14–17]. Several authors presented research work in shale gas reservoirs with geomechanical effects, such as Wang et al. [18], Lin et al. [19], and Yang et al. [20]. Wang et al. [18] developed a general model including embedded discrete fractures, multiple interacting continua and geomechanics in shale gas reservoirs with multi-scale fractures. In [19], authors provide insights on fracture propagation using reservoir flow simulation integrated with geomechanics in unconventional tight reservoirs. Effects of multicomponent adsorption and geomechanics are investigated in a shale condensate reservoir [20]. Girault et al. [21] have introduced a priori error estimates for a discretized poro-elastic-elastic system, in which the flow pressure equation is discretized by either a continuous Galerkin scheme or a mixed scheme, while, elastic displacement equations are discretized by a continuous Galerkin scheme. El-Amin et al. [22] have used the MFEM with stability analysis to simulate the problem of natural gas transport in a low-permeability reservoir without considering fractures. The authors [23] extended their work to cover fractured porous media and the rock stress-sensitivity with considered stability analysis of the MFEM. In this work, we present a theoretical basis with proofs of the stability analysis of the MFEM (in Ref. [23]) including the necessary lemmas and theorem.

2 Modeling and formulation

In this section, the mathematical model of the problem under consideration is developed. The Dual Porosity Dual Permeability (DPDP) model is employed to describe the gas transport in fractured porous media which consists of matrix blocks and fractures. The absorbed gas and the free gas coexist in the matrix blocks, however, the free gas exists in fractures only. The mass of free gas accumulation is ${\varphi }_m{\rho }_m$, and the adsorbed gas accumulation (Langmuir isotherm model) is [24, 25],$$\frac{(1-{\varphi }_m){\rho }_s{M}_w{V}_L{p}_m}{V_{\mathrm{std}}(P_L+p_m)}.$$(1)

The mass density of the gas can be written as,$${\rho }_{g,d}=\frac{p_d{M}_w}{\mathrm{ZRT}}, \hspace{1em}d=m,f.$$(2)

The gas compressibility factor, $Z$ which is given by Peng-Robinson equation of state [26],$$Z^3-(1-B)Z^2+(A-3B^3-2B)Z-(\mathrm{AB}-B^2-B^3)=0,$$(3)$$A=\frac{a_{T}p}{R^2{T}^2},\hspace{1em} B=\frac{b_{T}p}{\mathrm{RT}},$$(4)$$a_T=0.45724\frac{R^2{T}_c^2}{p_c}, \hspace{1em}{b}_T=0.0778\frac{R{T}_c}{p_c}.$$(5)

In low-permeability formations, Klinkenberg effect, which is described by the apparent permeability, takes place and written as,$$k_{\mathrm{app},d}=k_{0,d}\left(1+\frac{b_d}{p_d}\right),\hspace{1em} d=m,f.$$(6)

The DPDP model of gas transport in fracture shale strata is represented as,$$f_1\left(p_m\right)\frac{\mathrm{\partial }{p}_m}{\mathrm{\partial }t}-\nabla \cdot \left[\frac{{\rho }_{g,m}\left(p_m\right)}{\mu }{k}_{0,m}\left(1+\frac{b_m}{p_m}\right)\nabla p_m\right]=-S(p_m,p_f),$$(7)and$$f_2\left(p_f\right)\frac{\mathrm{\partial }{p}_f}{\mathrm{\partial }t}-\nabla \cdot \left[\frac{{\rho }_{g,f}\left(p_f\right)}{\mu }{k}_{0,f}\left(1+\frac{b_f}{p_f}\right)\nabla p_f\right]=S(p_m,p_f)-Q\left(p_f\right),$$(8)where$$f_1\left(p_m\right)=\frac{M_w}{\mathrm{ZRT}}\left[{\varphi }_m\left(p_m\right)+{\varphi }_m^{\text{'}}\left(p_m\right)p_m\right]+\frac{M_w{V}_L{\rho }_s}{V_{\mathrm{std}}}\left(\frac{P_L\left(1-{\varphi }_m\left(p_m\right)\right)}{{\left(P_L-p_m\right)}^2}-\frac{{\varphi }_m^{\text{'}}\left(p_m\right)}{P_L+p_m}\right),$$(9)and$$f_2\left(p_f\right)=\frac{M_w}{\mathrm{ZRT}}\left[{\varphi }_f\left(p_f\right)+{\varphi \text{'}}_f\left(p_f\right)p_f\right],$$(10)such that ${\varphi }_m$ and ${\varphi }_f$ is the matrix and fractures porosities, respectively,$${\varphi \text{'}}_d\left(p_d\right)=\frac{\mathrm{d}{\varphi }_d}{\mathrm{d}{p}_d}, \hspace{1em}d=m,f.$$(11)

Moreover, based on the following definitions, $b_m$ and $b_f$ are treated as constants [17],$$b_m=\sqrt{\frac{8\pi \mathrm{RT}}{M_w}}\frac{1}{r}\left(\frac{2}{\alpha }-0.995\right){\mu }_g,$$(12)and$$b_f=\sqrt{\frac{\pi \mathrm{RT}{\varphi }_{0,f}}{M_w{k}_{0,f}}} {\mu }_g.$$(13)

The Knudsen diffusion coefficient is defined as,$$D_{\mathrm{kf}}=\sqrt{\frac{\pi \mathrm{RT}{k}_{0,f}{\varphi }_{0,f}}{M_w}}.$$(14)Here, $S(p_m,p_f)$ is the transfer parameter which is connecting of the matrix-fracture domains, and defined as,$$S\left(p_m,p_f\right)=\frac{\sigma {\rho }_{g,m}{k}_m}{\mu }\left(p_m-p_f\right).$$(15)

One may provide the definition of the source of the production well by,$$Q\left(p_f\right)=\frac{\theta {\rho }_{g,f}{k}_f}{\mu \ln \frac{r_e}{r_w}}\left(p_{\mathrm{fe}}-p_w\right),$$(16)$$\theta =\{\begin{array}{l}\theta =2\pi \hspace{1em}\mathrm{if} \mathrm{the} \mathrm{production} \mathrm{well} \mathrm{at} \mathrm{the} \mathrm{center} \mathrm{of} \mathrm{the} \mathrm{reservoir}\\ \theta =\frac{\pi }{2}\hspace{1em}\mathrm{if} \mathrm{the} \mathrm{production} \mathrm{well} \mathrm{at} \mathrm{the} \mathrm{corner} \mathrm{of} \mathrm{the} \mathrm{reservoir}\end{array}.$$

Given the constant, $r_c$, the drainage-radius is represented as,$$r_e=0.14\sqrt{(\Delta x{)}^2+(\Delta y{)}^2}.$$(17)

The coefficient of crossflow between the matrix and fracture domains is defined as [9],$$\sigma =\frac{2n(n+1)}{l^2},$$(18)such that $l$ is given by,$$l=\frac{3L_x{L}_y{L}_z}{L_x{L}_y+L_y{L}_z+L_x{L}_z}.$$(19)

2.1 Geomechanical effects

Invoking the rock stress-sensitivity [27], the porosity can be expressed as a function of the mean effective-stress,$${{{\varphi }_d(\sigma \text{'}}_d\left(p_d\right))={\varphi }_{r,d}+({\varphi }_{0,d}-{\varphi }_{r,d}\left)\exp \right(-\mathrm{a\sigma }\text{'}}_d), \hspace{1em}d=m,f,$$(20)$${\sigma \text{'}}_d={\sigma }_d+{\alpha }_d{p}_d, \hspace{1em}d=m,f.$$(21)From (20) and (21), we can obtain,$${\varphi }_d\left(p_d\right)={\varphi }_{r,d}+C_{1,d}\exp (-a{p}_d), C_{1,d}=({\varphi }_{0,d}-{\varphi }_{r,d})\exp (-a{\sigma }_d), \hspace{1em}d=m,f.$$(22)

As the porosity increases the difference $({\varphi }_{0,d}-{\varphi }_{r,d})$ becomes positive and vice versa. Also, the coefficient $C_{1,d}$ becomes small based on the change in the porosity, and it becomes positive if porosity increases and negative if porosity decreases.

2.2 Initial and boundary conditions

The initial conditions can be represented by,$$p_m(\cdot ,0)=p_f(\cdot ,0)=p_0\hspace{1em}\mathrm{i}\mathrm{n}\hspace{1em} {\mathrm{\Omega }}_m\cup {\mathrm{\Omega }}_f.$$(23)

The pressure on the production well (boundary condition) is given as,$$p_f\left(\cdot ,t\right)=p_w\hspace{1em}\mathrm{o}\mathrm{n}\hspace{1em}{\mathrm{\Gamma }}_D^f\times \left(0,T\right).$$(24)

The no-flow boundary conditions on matrix are given as,$$u_m\cdot n=0\hspace{1em}\mathrm{o}\mathrm{n}\hspace{1em}{\mathrm{\Gamma }}_N^m\cup {\mathrm{\Gamma }}_N^f\times (0,T),$$(25)while the no-flow boundary conditions for fracture domain are,$$u_f\cdot n=0\hspace{1em}\mathrm{o}\mathrm{n}\hspace{1em}{\mathrm{\Gamma }}_N^f\times \left(0,T\right).$$(26)

3 MFEM spaces

The MFEM contains two spaces for a scalar variable and it’s flux. The MFEM was developed to approximate both variables simultaneously and to give a higher order approximation for the flux. A compatibility condition must hold for the two spaces to insure stability, consistency and convergence of the mixed method and by considering more constraints to the numerical discretization.

Now, let us define the inner product in $\mathrm{\Omega }$ as,$$(f,g{)}_{\mathrm{\Omega }}={\int }_{\mathrm{\Omega }}\mathrm{}f(\mathbf{x}\left)g\right(\mathbf{x})\mathrm{d}\mathbf{V} \hspace{1em}\forall f,g:\mathrm{\Omega }\to \mathbb{R},$$(27)and the inner product on $\mathrm{\partial }\mathrm{\Omega }$ as,$$\langle f,g{\rangle }_{\mathrm{\partial }\mathrm{\Omega }}={\int }_{\mathrm{\partial }}\mathrm{}\mathrm{fg}\mathrm{d}S \hspace{1em}\forall f,g:\mathrm{\partial }\mathrm{\Omega }\to \mathbb{R}.$$(28)Given,$${\mathbb{L}}^2\left(\mathrm{\Omega }\right)=\{f:\mathrm{\Omega }\to \mathbb{R}:{\int }_{\mathrm{\Omega }}\mathrm{}{f}^2\mathrm{d}x<+\infty \},$$(29)is the largest Hilbert space, such that $(D^{-1}\mathbf{u},w)$ and $(p,\nabla \cdot w)$ are well defined, therefore, $p,\nabla \cdot w\in {\mathbb{L}}^2\left(\mathrm{\Omega }\right)$. It is needed that $p,\phi \in {\mathbb{L}}^2\left(\mathrm{\Omega }\right)$ and $\mathbf{u},w\in \mathbb{H}(\mathrm{div},\mathrm{\Omega })$. The Hilbert space with norm given by,$$\mathbb{H}(\mathrm{div},\mathrm{\Omega })=\{\mathbf{u}:\nabla \cdot \mathbf{u}\in {\mathbb{L}}^2(\mathrm{\Omega }\left)\right\},$$(30)$${\Vert \mathbf{u}\Vert }_{\mathbf{H}(\mathrm{d}\mathrm{iv},\mathrm{\Omega })}^2={\Vert \mathbf{u}\Vert }^{\mathbf{2}} + {\Vert \nabla \cdot \mathbf{u}\Vert }^2.$$(31)

The seeking solution is,$$(p,\mathbf{u})\in {\mathbb{L}}^2\left(\mathrm{\Omega }\right)\times \mathbb{H}(\mathrm{div},\mathrm{\Omega }),$$(32)such that $p$ and $\mathbf{u}$ are smooth. For the approximate numerical solution, the two spaces become,$$W_h\subset {\mathbb{L}}^2\left(\mathrm{\Omega }\right), \hspace{1em}\dim \left(W_h\right)<+\infty$$(33)and,$$V_h\subset \mathbb{H}\left(\mathrm{div},\mathrm{\Omega }\right),\hspace{1em}\dim \left(V_h\right)<+\infty .$$(34)

Thus, the normal component $\mathbf{u}\cdot \mathbf{n}$ is continuous across the inter-element boundaries.

Moreover, it is useful to present the $\mathbf{R}{\mathbf{T}}_r$ elements that are designed to approximate $\mathbb{H}(\mathrm{div},\mathrm{\Omega })$ [7], which satisfies,$$W_h=\left\{w\in {\mathbb{L}}^2\left(\mathrm{\Omega }\right):w{|}_E\in {\mathbb{P}}_0\left(E\right),E\in {\mathcal{E}}_h\right\},$$(35)and,$$V_h=\left\{\mathbf{u}\in \mathbb{H}\left(\mathrm{div},\mathrm{\Omega }\right):\mathbf{u}{|}_E\in {\mathbb{P}}_1\left(E\right),E\in {\mathcal{E}}_h\right\},$$(36)where $w$ is discontinuous piecewise constant and $\mathbf{u}$ is piecewise linear.

4 Mixed finite element approximation

Let ${\mathrm{\Omega }}_m$ and ${\mathrm{\Omega }}_f$ are, respectively, the matrix and fracture in a polygonal/polyhedral Lipschitz domain $\mathrm{\Omega }\subset {\mathbb{R}}^d,\hspace{1em}d\in \{\mathrm{1,2},3\}$ (on which we define the standard $L^2\left(\mathrm{\Omega }\right)$ space such that ${\mathbb{L}}^2\left(\mathrm{\Omega }\right)\equiv {\left(L^2\left(\mathrm{\Omega }\right)\right)}^d$), with the boundaries, $\mathrm{\partial }{\mathrm{\Omega }}_m={\mathrm{\Gamma }}_D^m\cup {\mathrm{\Gamma }}_N^m$ and $\mathrm{\partial }{\mathrm{\Omega }}_f={\mathrm{\Gamma }}_D^f\cup {\mathrm{\Gamma }}_N^f$. One may write the above dual-porosity model in the following general form,$$f_1\left(p_m\right)\frac{\mathrm{\partial }{p}_m}{\mathrm{\partial }t}+\nabla \cdot u_m=-S\left(p_m,p_f\right)\hspace{1em}\mathrm{i}\mathrm{n}\hspace{1em}{\mathrm{\Omega }}_m\times (0,T),$$(37)$${\mathbf{D}}_m(p_m{)}^{-1}{u}_m=-\nabla p_m\hspace{1em}\mathrm{i}\mathrm{n}\hspace{1em}{\mathrm{\Omega }}_m\times (0,T),$$(38)$$f_2\left(p_f\right)\frac{\mathrm{\partial }{p}_f}{\mathrm{\partial }t}+\nabla \cdot u_f=S(p_m,p_f)-Q\left(p_f\right)\hspace{1em}\mathrm{i}\mathrm{n}\hspace{1em}{\mathrm{\Omega }}_f\times (0,T),$$(39)$${\mathbf{D}}_f(p_f{)}^{-1}{u}_f=-\nabla p_f\hspace{1em}\mathrm{i}\mathrm{n}\hspace{1em}{\mathrm{\Omega }}_f\times (0,T),$$(40)where$${\mathbf{D}}_m\left(p_m\right)=\frac{\rho \left(p_m\right)}{\mu }{k}_{0,m}\left(1+\frac{b_m}{p_m}\right),$$(41)and$${\mathbf{D}}_f\left(p_f\right)=\frac{\rho \left(p_f\right)}{\mu }{k}_{0,f}\left(1+\frac{b_f}{p_f}\right).$$(42)

The functions ${\mathbf{D}}_m(p_m{)}^{-1}$ and ${\mathbf{D}}_f(p_f{)}^{-1}$ are moved to the left hand side to avoid discontinuity when we integrate $\nabla p_m$ and $\nabla p_f$ by parts. Selecting any $\phi \in W_h$ and $\omega \in V_h$, the mixed finite element weak formulation can be written in the following form,$$\left(f_1\left(p_m\right)\frac{\mathrm{\partial }{p}_m}{\mathrm{\partial }t},\phi \right)+(\nabla \cdot {\mathbf{u}}_m,\phi )+\left(\mathcal{S}\right(p_m,p_f),\phi )=0,$$(43)$$\left({\mathbf{D}}_m\right(p_m{)}^{-1}{\mathbf{u}}_m,\omega )=(p_m,\nabla \cdot \omega ),$$(44)$$\left(f_2\left(p_f\right)\frac{\mathrm{\partial }{p}_f}{\mathrm{\partial }t},\phi \right)+(\nabla \cdot {\mathbf{u}}_f,\phi )-\left(\mathcal{S}\right(p_m,p_f),\phi )=-\left(Q\right(p_f),\phi ),$$(45)$$\left({\mathbf{D}}_f\right(p_f{)}^{-1}{\mathbf{u}}_f,\omega )=(p_f,\nabla \cdot \omega )-{\langle p_w,\omega \rangle }_{{\mathrm{\Gamma }}_f^D}.$$(46)

Now, let the approximating subspace duality $V_h\subset H(\mathrm{\Omega };\mathrm{div})$ and $W_h\subset L^2\left(\mathrm{\Omega }\right)$ be the $r$-th order ($r\ge 0$) Raviart–Thomas space (RT_r) on the partition ${\mathcal{T}}_h$. The mixed finite element formulations are stated as below: find $p_m^h,p_f^h\in W_h$ and $u_m^h,u_f^h\in V_h$ such that,$$\left(f_1\left(p_m^h\right)\frac{\mathrm{\partial }{p}_m^h}{\mathrm{\partial }t},\phi \right)+(\nabla \cdot {\mathbf{u}}_m^h,\phi )+\left(\mathcal{S}\right(p_m^h,p_f^h),\phi )=0,$$(47)$$\left({\mathbf{D}}_m\right(p_m^h{)}^{-1}{\mathbf{u}}_m^h,\omega )=(p_m^h,\nabla \cdot \omega ),$$(48)$$\left(f_2\left(p_f^h\right)\frac{\mathrm{\partial }{p}_f^h}{\mathrm{\partial }t},\phi \right)+(\nabla \cdot {\mathbf{u}}_f^h,\phi )-\left(\mathcal{S}\right(p_m^h,p_f^h),\phi )=-\left(Q\right(p_f^h),\phi ),$$(49)$$\left({\mathbf{D}}_f\right(p_f^h{)}^{-1}{\mathbf{u}}_f^h,\omega )=(p_f^h,\nabla \cdot \omega )-{\langle p_w,\omega \rangle }_{{\mathrm{\Gamma }}_f^D},$$(50)for any $\phi \in W_h$ and $\omega \in V_h$.

In order to obtain an explicit formulation for the flux, we employ a quadrature rule along with the MFE method to [8]. The total time interval $[0,T]$ is divided into $N_T$ time steps with length $\Delta t^n=t^{n+1}-t^n$. The superscript $n+1$ denotes for the current time step, while $n$ denotes for the previous one. We use a backward Euler semi-implicit discretization for the time derivative terms. The following scheme has been developed,$$\left(f_1\left(p_m^{h,n}\right)\frac{p_m^{h,n+1}-p_m^{h,n}}{\Delta t},\phi \right)+(\nabla \cdot {\mathbf{u}}_m^{h,n+1},\phi )+\left(\mathcal{S}\right(p_m^{h,n+1},p_f^{h,n}),\phi )=0,$$(51)$$\left({\mathbf{D}}_m\right(p_m^{h,n}{)}^{-1}{\mathbf{u}}_m^{h,n+1},\omega )=(p_m^{h,n+1},\nabla \cdot \omega ),$$(52)$$\left(f_2\left(p_f^{h,n}\right)\frac{p_f^{h,n+1}-p_f^{h,n}}{\Delta t},\phi \right)+(\nabla \cdot {\mathbf{u}}_f^{h,n+1},\phi )-\left(\mathcal{S}\right(p_m^{h,n+1},p_f^{h,n+1}),\phi )=-\left(Q\right(p_f^{h,n+1}),\phi ),$$(53)$$\left({\mathbf{D}}_f\right(p_f^{h,n}{)}^{-1}{\mathbf{u}}_f^{h,n+1},\omega )=(p_f^{h,n+1},\nabla \cdot \omega )-{\langle p_w,\omega \rangle }_{{\mathrm{\Gamma }}_f^D}.$$(54)Given $p_m^{h,n}$ and $p_f^{h,n}$, the numerical procedure to calculate pressure and velocity is presented here,

1. Update the thermodynamical variables explicitly.

2. Solve the equations (51) and (52) to obtain $p_m^{h,n+1}$ and ${\mathbf{u}}_m^{h,n+1}$.

3. Solve the equations (53) and (54) to obtain $p_f^{h,n+1}$ and ${\mathbf{u}}_f^{h,n+1}$.

5 Stability analysis

In this section, we carry out the stability analysis of the proposed MFE method, which ensures that the discrete solutions are bounded in a physically reasonable range. The key issue encountered in the stability analysis is the nonlinearity of the matrix and fracture pressures. In order to resolve this issue, we need to define some auxiliary functions and analyze their boundedness. Now let us define the following functions,$$F_1\left(p_d^h\right)={\int }_0^{p_d^h}\mathrm{}{f}_1\left(p_d^h\right)\mathrm{d}{p}_d^h, G_1\left(p_d^h\right)={\int }_0^{p_d^h}\mathrm{}{F}_1\left(p_d^h\right)\mathrm{d}{p}_d^h,$$(55)$$F_2\left(p_d^h\right)={\int }_0^{p_d^h}\mathrm{}{f}_2\left(p_d^h\right)\mathrm{d}{p}_d^h, G_2\left(p_d^h\right)={\int }_0^{p_d^h}\mathrm{}{F}_2\left(p_d^h\right)\mathrm{d}{p}_d^h.$$(56)

Therefore, one may$$\begin{array}{l}{F}_1\left(p_m^h\right)=\left({\varphi }_{r,m}+C_{1,m}\right)\frac{M_w}{\mathrm{ZRT}}{p}_m^h{e}^{-a{p}_m^h}+\frac{M_w{V}_L{\rho }_s}{V_{\mathrm{std}}}[C_{1,m}\left(a{e}^{a{P}_L}\left(1+P_L\right)\mathrm{Ei}(-x){|}_{a{P}_L}^{a(P_L+p_m^h)}+\frac{P_L{e}^{-a{p}_m^h}}{P_L+p_m^h}-1\right)+(1-{\varphi }_{r,m})\frac{p_m^h}{P_L+p_m^h}],\\ \end{array}$$(57)$$\begin{array}{l}{G}_1\left(p_m^h\right)=\left({\varphi }_{r,m}+C_{1,m}\right)\frac{M_w}{a^2\mathrm{ZRT}}\left[1-(a{p}_m^h+1)e^{-a{p}_m^h}\right]+\frac{M_w{V}_L{\rho }_s}{V_{\mathrm{std}}}\left[C_{1,m}\left(a{e}^{a{P}_L}\left(1+P_L\right)\left((P_L+p_m^h)\mathrm{Ei}(-x){|}_{a{P}_L}^{a(P_L+p_m^h)}-\frac{1}{a}(1-e^{-a{p}_m^h})\right)P_L{e}^{a{P}_L}+\mathrm{Ei}(-x){|}_{a{P}_L}^{a(P_L+p_m^h)}\right)+(1-{\varphi }_{r,m}-C_{1,m})p_m^h+(1-{\varphi }_{r,m})P_L\ln \left(\frac{P_L}{P_L+p_m^h}\right)\right],\\ \\ \end{array}$$(58)$$F_2\left(p_f^h\right)=\left({\mathrm{\varphi }}_{r,f}+C_{1,f}\right)\frac{M_w}{\mathrm{ZRT}}{p}_f^h{e}^{-a{p}_f^h},$$(59)and$$G_2\left(p_f^h\right)=\left({\varphi }_{r,f}+C_{1,f}\right)\frac{M_w}{a^2\mathrm{ZRT}}\left[1-(a{p}_f^h+1)e^{-a{p}_f^h}\right].$$(60)

$\mathrm{Ei}\left(x\right)$ is a special function called the exponent integral function. As stated above in the model formulation that the coefficient $C_{1,d},\hspace{1em}d=m,f$ is positive in the case of increasing porosity and negative in the case of decreasing porosity.

Lemma 5.1

For the case of increasing matrix porosity, $C_{1,m}>0$, and sufficiently large $p_m^h$, there exists a positive constant,$${\gamma }_1=\left({\varphi }_{r,m}+C_{1,m}\right)C_4^l\frac{M_w}{\mathrm{ZRT}}+\frac{M_w{V}_L{\rho }_s}{V_{\mathrm{std}}}(1-{\varphi }_{r,m})C_3^l,$$(61)such that,$$\left(F_1\right(p_m^h),p_m^h)\ge {\gamma }_1{\Vert p_m^h\Vert }^2.$$(62)

Proof. Note that the value of the exponent integral function $\mathrm{Ei}(-x)$ is negative and (small for big values of $x$), therefore the quantity, $\mathrm{Ei}(-x){|}_{a{P}_L}^{a(P_L+p_m^h)}$ has always a positive value which has lower and upper bounds, namely,$$C_{2,m}^u\ge \mathrm{Ei}\left(-x\right){|}_{a{P}_L}^{a\left(P_L+p_m^h\right)}\ge C_{2,m}^l,$$(63)where $C_{2,m}^l,C_{2,m}^u>0$. On the other hand, in shale reservoir, it is well known that the initial pressure has the maximum value, i.e., $\max \left(p_m^h\right)=P_0$, while, the pressure of the well has the minimum value, i.e., $\min \left(p_m^h\right)=P_w$. Therefore, we have,$$C_3^u=\frac{1}{P_L+p_w}\ge \frac{1}{P_L+p_m^h}\ge \frac{1}{P_L+p_0}=C_3^l,$$(64)$$C_4^u=e^{-a{p}_w}\ge e^{-a{p}_m^h}\ge e^{-a{p}_0}=C_4^l,$$(65)and,$$C_5^u=\frac{e^{-a{p}_w}}{P_L+p_w}\ge \frac{e^{-a{p}_m^h}}{P_L+p_m^h}\ge \frac{e^{-a{p}_0}}{P_L+p_0}=C_5^l,$$(66)where $C_3^l,C_3^u,C_4^l,C_4^u,C_5^l,C_5^u\ge 0$. Therefore, for the case of increasing matrix porosity, $C_{1,m}>0$, and sufficiently large $p_m^h$, and holding (63)–(66), the following coefficient is positive, i.e.,$$a{e}^{a{P}_L}\left(1+P_L\right)\mathrm{Ei}\left(-x\right){|}_{a{P}_L}^{a\left(P_L+p_m^h\right)}+\frac{P_L{e}^{-a{p}_m^h}}{P_L+p_m^h}-1>,$$(67)$$a{e}^{a{P}_L}\left(1+P_L\right)C_{2,m}^l+P_L{C}_5^l-1=C_{3,m}>0.$$(68)Therefore,$$\begin{array}{cc}\left(F_1\right(p_m^h),p_m^h)& \ge ({\varphi }_{r,m}+C_{1,m})C_4^l\frac{M_w}{\mathrm{ZRT}}(p_m^h,p_m^h)+\frac{M_w{V}_L{\rho }_s}{V_{\mathrm{std}}}\left[C_{1,m}{C}_{3,m}+(1-{\varphi }_{r,m})C_3^l(p_m^h,p_m^h)\right]\\ \begin{array}{c}\\ \end{array}& \begin{array}{c}\ge \left({\varphi }_{r,m}+C_{1,m}\right)C_4^l\frac{M_w}{\mathrm{ZRT}}(p_m^h,p_m^h)+\frac{M_w{V}_L{\rho }_s}{V_{\mathrm{std}}}(1-{\varphi }_{r,m})C_3^l(p_m^h,p_m^h)\\ =\left[\left({\varphi }_{r,m}+C_{1,m}\right)C_4^l\frac{M_w}{\mathrm{ZRT}}+\frac{M_w{V}_L{\rho }_s}{V_{\mathrm{std}}}(1-{\mathrm{\varphi }}_{r,m})C_3^l\right]{\Vert p_m^h\Vert }^2={\gamma }_1{\Vert p_m^h\Vert }^2,\end{array}\end{array}$$(69)and this completes Lemma 5.1 proof.

Lemma 5.2

For the case of decreasing matrix porosity, $C_{1,m}<0$, assuming,$$a{e}^{a{P}_L}\left(1+P_L\right)\mathrm{Ei}(-x){|}_{a{P}_L}^{a\left(P_L+p_m^h\right)}+\frac{P_L{e}^{-a{p}_m^h}}{P_L+p_m^h}-1<0,$$(70)there exists a positive constant ${\gamma }_1$, such that,$$\left(F_1\right(p_m^h),p_m^h)\ge {\gamma }_1{\Vert p_m^h\Vert }^2.$$(71)

Proof. It is clear that, $C_{1,m}<{\varphi }_{r,m}$, so, we always have, ${\varphi }_{r,m}-C_{1,m}>0$. Holding the assumption (70), one can easily prove this lemma in a similar way to Lemma 5.1.

Lemma 5.3

For the case of increasing matrix porosity, $C_{1,m}>0$, and for sufficiently large $p_m^h$, there exist two positive constants,$$\begin{array}{l}{\gamma }_2=\frac{M_w}{a^2\mathrm{ZRT}}\left({\varphi }_{r,m}+C_{1,m}\right)\left(1-C_{4,m}^u\right)+\frac{M_w{V}_L{\rho }_s}{V_{\mathrm{std}}}[C_{1,m}(a{e}^{a{P}_L}\left(1+P_L\right)P_L{C}_{2,m}^u+P_L{e}^{a{P}_L}{C}_{2,m}^u)+\left(1-{\varphi }_{r,m}\right)P_L\ln \left(P_L{C}_3^u\right)],\\ \end{array}$$(72)and$${\gamma \text{'}}_2^{}=\frac{M_w{V}_L{\rho }_s}{V_{\mathrm{std}}}\left[C_{1,m}a{e}^{a{P}_L}(1+P_L)C_{2,m}^u+(1-{\varphi }_{r,m}-C_{1,m})\right],$$(73)such that,$$\left(G_1\left(p_m^h\right),1\right)\le {\gamma }_2\left|\mathrm{\Omega }\right|+{{\gamma }^{\text{'}}}_2^{}\Vert p_m^h\Vert .$$(74)

Proof. Holding (63)–(66) with $C_{4,m}^u=\frac{a{p}_m^h+1}{e^{a{p}_m^h}}<1$ for sufficiently large $p_m^h$, and integrate (58) over ${\mathrm{\Omega }}_m$ for the case of increasing porosity, $C_{1,m}>0$, it is easy to prove (74), which completes the lemma proof.

Lemma 5.4

For the case of decreasing matrix porosity, $C_{1,m}<0$, with assuming that,$$C_{1,m}\left(a{e}^{a{P}_L}(1+P_L)P_L{C}_{2,m}^u+P_L{e}^{a{P}_L}{C}_{2,m}^u\right)+(1-{\varphi }_{r,m})P_L\ln \left(P_L{C}_3^u\right)\ge 0,$$(75)and$$C_{1,m}a{e}^{a{P}_L}(1+P_L)C_{2,m}^u+(1-{\varphi }_{r,m}-C_{1,m})\ge 0,$$(76)one can prove,$$\left(G_1\right(p_m^h),1)\le {\gamma }_2\left|\mathrm{\Omega }\right|+{\gamma \text{'}}_2^{}\Vert p_m^h\Vert .$$(77)