Theoretical stability analysis of mixed ﬁ nite element model of shale-gas ﬂ ow with geomechanical effect

. In this work, we introduce a theoretical foundation of the stability analysis of the mixed ﬁ nite element solution to the problem of shale-gas transport in fractured porous media with geomechanical effects. The differential system was solved numerically by the Mixed Finite Element Method (MFEM). The results include seven lemmas and a theorem with rigorous mathematical proofs. The stability analysis presents the boundedness condition of the MFE solution.


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 dualcontinuum model have been developed to include, for instance, adsorption, heat, deformation and Knudsen diffusion [14][15][16][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.

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 / m q m , and the adsorbed gas accumulation (Langmuir isotherm model) is [24,25], The mass density of the gas can be written as, The gas compressibility factor, Z which is given by Peng-Robinson equation of state [26], In low-permeability formations, Klinkenberg effect, which is described by the apparent permeability, takes place and written as, The DPDP model of gas transport in fracture shale strata is represented as, and where and such that / m and / f is the matrix and fractures porosities, respectively, Moreover, based on the following definitions, b m and b f are treated as constants [17], The Knudsen diffusion coefficient is defined as, Here, Sðp m ; p f Þ is the transfer parameter which is connecting of the matrix-fracture domains, and defined as, One may provide the definition of the source of the production well by, Given the constant, r c , the drainage-radius is represented as, The coefficient of crossflow between the matrix and fracture domains is defined as [9], such that l is given by,

Geomechanical effects
Invoking the rock stress-sensitivity [27], the porosity can be expressed as a function of the mean effective-stress, From (20) and (21), we can obtain, As the porosity increases the difference ð/ 0;d À / 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.

Initial and boundary conditions
The initial conditions can be represented by, The pressure on the production well (boundary condition) is given as, The no-flow boundary conditions on matrix are given as, while the no-flow boundary conditions for fracture domain are,

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 X as, and the inner product on oX as, Given, is the largest Hilbert space, such that ðD À1 u; wÞ and ðp; r Á wÞ are well defined, therefore, p; r Á w 2 L 2 ðXÞ. It is needed that p; u 2 L 2 ðXÞ and u; w 2 Hðdiv; XÞ. The Hilbert space with norm given by, The seeking solution is, ðp; uÞ 2 L 2 ðXÞ Â Hðdiv; XÞ; such that p and u are smooth. For the approximate numerical solution, the two spaces become, and, Thus, the normal component u Á n is continuous across the inter-element boundaries. Moreover, it is useful to present the RT r elements that are designed to approximate Hðdiv; XÞ [7], which satisfies, and, where w is discontinuous piecewise constant and u is piecewise linear.

Mixed finite element approximation
Let X m and X f are, respectively, the matrix and fracture in a polygonal/polyhedral Lipschitz domain X & R d ; d 2 f1; 2; 3g (on which we define the standard L 2 ðXÞ space such that L 2 ðXÞ L 2 ðXÞ À Á d ), with the boundaries, oX m ¼ One may write the above dual-porosity model in the following general form, where and The functions D m ðp m Þ À1 and D f ðp f Þ À1 are moved to the left hand side to avoid discontinuity when we integrate rp m and rp f by parts. Selecting any u 2 W h and x 2 V h , the mixed finite element weak formulation can be written in the following form, ð45Þ for any u 2 W h and x 2 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 Á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, Given p h;n m and p h;n f , the numerical procedure to calculate pressure and velocity is presented here,

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, Therefore, one may and EiðxÞ is a special function called the exponent integral function. As stated above in the model formulation that the coefficient C 1;d ; 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 h m , there exists a positive constant, such that, Proof. Note that the value of the exponent integral function EiðÀxÞ is negative and (small for big values of x), therefore the quantity, EiðÀxÞj aðP L þp h m Þ aP L has always a positive value which has lower and upper bounds, namely, where C l 2;m ; C u 2;m > 0. On the other hand, in shale reservoir, it is well known that the initial pressure has the maximum value, i.e., maxðp h m Þ ¼ P 0 , while, the pressure of the well has the minimum value, i.e., minðp h m Þ ¼ P w . Therefore, we have, and, where C l 3 ; C u 3 ; C l 4 ; C u 4 ; C l 5 ; C u 5 ! 0. Therefore, for the case of increasing matrix porosity, C 1;m > 0, and sufficiently large p h m , and holding (63)-(66), the following coefficient is positive, i.e., Therefore, and this completes Lemma 5.1 proof.
Lemma 5.2 For the case of decreasing matrix porosity, C 1;m < 0, assuming, there exists a positive constant c 1 , such that, Proof. It is clear that, C 1;m < / r;m , so, we always have, / 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 h m , there exist two positive constants, and such that, Proof. Holding (63)-(66) with C u 4;m ¼ ap h m þ1 e ap h m < 1 for sufficiently large p h m , and integrate (58) over X 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 ae aP L ð1 þ P L ÞP L C u 2;m þ P L e aP L C u and C 1;m ae aP L ð1 þ P L ÞC u 2;m þ ð1 À / r;m À C 1;m Þ ! 0; ð76Þ one can prove, Proof. It is clear that, C 1;m < / r;m , so,we always have, / r;m À C 1;m > 0. Holding the assumptions (75) and (76), therefore we find, c 2 > 0; c 0 2 > 0, and integrate (58) over X m for the case of decreasing porosity, C 1;m < 0, which proves (74), and this completes the proof. Lemma 5.5 For both cases of increasing or decreasing fracture porosity, namely, C 1;f > 0 or C 1;f < 0, and sufficiently large p h f , there exists a positive constant, Proof. Again, we have, where C l 4 ; C u 4 > 0. Therefore, for the case of increasing fracture porosity, C 1;f > 0, and sufficiently large p h f , Therefore, As stated above, for the case of the decreasing porosity, C 1;f < 0, the coefficient ð/ r;f þ C 1;f Þ remains positive, then the above inequality holds. This completes the lemma proof.
Lemma 5.6 For both cases of increasing or decreasing fracture porosity, namely, C 1;f > 0 or C 1;f < 0, and sufficiently large p h f , there exists a positive constant, such that, Proof. Similar to the previous proofs.
Proof. Since p fe is bounded by p w and p 0 , i.e., p w p fe p 0 , then, 0 < p fe À p w p 0 À p w ¼ C p , such that C p is positive number. Holding the conditions (84), one may find that, Moreover, suppose that p fe ; p w 2 L 2 ð0; T ; L 2 ðXÞÞ. Then,