*Article* **Estimation of Free and Adsorbed Gas Volumes in Shale Gas Reservoirs under a Poro-Elastic Environment**

**Reda Abdel Azim 1,\*, Abdulrahman Aljehani <sup>2</sup> and Saad Alatefi <sup>3</sup>**


**\*** Correspondence: reda.abdulrasoul@auk.edu.krd

**Abstract:** Unlike conventional gas reservoirs, fluid flow in shale gas reservoirs is characterized by complex interactions between various factors, such as stress sensitivity, matrix shrinkage, and critical desorption pressure. These factors play a crucial role in determining the behavior and productivity of shale gas reservoirs. Stress sensitivity refers to the stress changes caused by formation pressure decline during production, where the shale gas formation becomes more compressed and its porosity decreases. Matrix shrinkage, on the other hand, refers to the deformation of the shale matrix due to the gas desorption process once the reservoir pressure reaches the critical desorption pressure where absorbed gas molecules start to leave the matrix surface, causing an increase in shale matrix porosity. Therefore, the accurate estimation of gas reserves requires careful consideration of such unique and complex interactions of shale gas flow behavior when using a material balance equation (MBE). However, the existing MBEs either neglect some of these important parameters in shale gas reserve analysis or employ an iterative approach to incorporate them. Accordingly, this study introduces a straightforward modification to the material balance equation. This modification will enable more accurate estimation of shale gas reserves by considering stress sensitivity and variations in porosity during shale gas production and will also account for the effect of critical desorption pressure, water production, and water influx. By establishing a linear relationship between reservoir expansion and production terms, we eliminate the need for complex and iterative calculations. As a result, this approach offers a simpler yet effective means of estimating shale gas reserves without compromising accuracy. The proposed MBE was validated using an in-house finite element poro-elastic model which accounts for stress re-distribution and deformation effects during shale gas production. Moreover, the proposed MBE was tested using real-field data of a shale gas reservoir obtained from the literature. The results of this study demonstrate the reliability and usefulness of the modified MBE as a tool for accurately assessing free and adsorbed shale gas volumes.

**Keywords:** shale gas reservoirs; poro-elasticity is shale gas; gas absorption–desorption process

### **1. Introduction**

The complexities involved in shale gas systems present challenges when it comes to estimating the volume and production of shale gas. This is primarily due to factors such as low permeability and highly heterogeneous porous media, which make conventional estimation methods ineffective. Therefore, a crucial step in determining shale gas production is selecting the appropriate material balance equation. Gas not only exists in small pores and open fractures but is also adsorbed into the grains of shale. Previous research and field data indicate that over 50% of the total gas is adsorbed within the matrix of the shale formation. The Langmuir pressure and volume parameters are used to estimate this volume of adsorbed gas [1–8].

Research by Darishchev et al. 2013 [9] demonstrates how the gas desorption volume is measured using shale samples at different temperatures. The study findings provide

**Citation:** Azim, R.A.; Aljehani, A.; Alatefi, S. Estimation of Free and Adsorbed Gas Volumes in Shale Gas Reservoirs under a Poro-Elastic Environment. *Energies* **2023**, *16*, 5798. https://doi.org/10.3390/en16155798

Academic Editors: Xingguang Xu, Kun Xie, Yang Yang and Shu Tao

Received: 15 May 2023 Revised: 18 July 2023 Accepted: 27 July 2023 Published: 4 August 2023

**Copyright:** © 2023 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

clear insights into the requirement of maintaining reservoir pressure below the Langmuir pressure to effectively release adsorbed gas. In rich shale formations, the overall recovery of gas heavily relies on the quantity of adsorbed gas that can be liberated. Despite hydraulic fracturing being employed to enhance permeability in shale matrices, there remains a portion of residual absorbed gas that cannot be extracted without utilizing additional high-value stimulation techniques like thermal stimulation [10–12].

As the main gas production contribution from the adsorbed gas, neglecting the desorption process may lead to unrealistic estimations of well production rates [13–18]. The desorption of gas is a vital factor in the deformation of the solid matrix of a reservoir. Therefore, this study takes this process into account using the poro-elasticity theory based on finite element modeling. Unlike conventional reservoirs, shale gas reservoirs require extraordinary technologies to extract the gas at an economical rate due to their extremely low matrix permeability and porosity. Accordingly, modeling shale gas well production and estimating accurate gas in place volumes, both free and adsorbed, are essential to predict the economic viability of the gas production process.

Shale gas reservoirs are different than the conventional gas reservoirs, in that they contain not only free gas, but also adsorbed gas and dissolved gas [19–21]. King (1990, 1993) [22,23] developed an MBE for a shale gas reservoir with the condition of limited water influx. The author considered the adsorption/desorption characteristics of the adsorbed gas and free gas. In addition, the author constructed a linear relationship through a pseudodeviation factor to estimate the original gas in place (OGIP). However, the work of King required an intensive iteration process. Clarkson and Smith [24] (1997) considered gas in the adsorbed state only when deriving an MBE for a coalbed methane (CBM) reservoir. Liu et al. (2020) [25] considered free gas, formation compressibility, gas adsorption/desorption, and water expansion when proposing an MBE for a CBM reservoir. On the other hand, Moghadam et. al. (2011) [26] ignored water and formation compressibility when deriving the MBE for CBM and shale gas reservoirs. Meng et al. (2017) [27] used the same principle as King's to derive an MBE and presented a new definition of the Z factor. The original adsorbed gas in place and original free gas in place were estimated [28–37] using a developed MBE for a CBM reservoir. Moreover, Shi et al. (2018a) proposed a modified material balance equation that incorporates the effects of gas desorption and matrix shrinkage, using a pseudo-deviation factor. However, determining this pseudo-deviation factor is complex and requires an iterative process. More recently, Meng et al., Shi et al., and Yang et al. [38–40] also developed modified MBEs for shale gas reservoirs that account for critical desorption pressure, stress sensitivity, and matrix shrinkage effects. Nevertheless, these equations still rely on the use of a pseudo-deviation factor which requires an iteration process. In summary, it can be concluded that current MBE formulations either neglect important parameters in shale gas reservoir analysis (such as stress sensitivity, matrix shrinkage, and critical desorption pressure), or they include these factors but require iterative calculations. Therefore, the primary objective of this study was to introduce a straightforward adjustment to the material balance equation. This modification will enable more accurate estimation of free and adsorbed shale gas volumes by considering stress sensitivity and variations in porosity during shale production, along with accounting for the effect of critical desorption pressure. By establishing a linear relationship between reservoir expansion and production terms, we eliminate the need for complex and iterative calculations. As a result, this approach offers a simpler yet effective means of estimating shale gas reserves without compromising accuracy. The structure of this study is as follows: the derivation of the proposed modified shale gas material balance equation is presented in Section 2. Then, an introduction to the in-house finite element simulator which accounts for the coupled poro-elasticity (deformation-diffusion process) usually encountered in shale gas reservoirs flow behavior is presented in Section 3. Finally, verification of the modified MBE using the in-house poro-elastic model and real-field cases are presented in Section 4.

#### **2. General Material Balance Equation of the Shale Gas Reservoir**

The derivation of the modified shale gas material balance equation (MBE) is presented in this section. The proposed MBE takes into consideration three distinct forms in which gas is stored: free, absorbed, and dissolved gases. However, the amount of gas dissolved in formation water can be ignored as it is negligible when compared to the volumes of free and absorbed gases. When the pressure decreases below a certain point, known as the critical desorption pressure, the organic matter releases (desorbs) its stored gas into matrix pores causing deformation in the shale due to stress changes and matrix shrinkage. Additionally, its assumed that during production, there is an isothermal process occurring and natural fractures act primarily as pathways rather than storage spaces.

The material balance equation can be expressed as

$$\begin{array}{l}\text{Production} = \text{Expansion (including expansion in both Organcic and Inor} \\ \text{gonic matters)} \end{array} \tag{1}$$

The production term in Equation (1) is stated as

$$G\_{\mathcal{P}}\beta\_{\mathcal{S}} + \mathcal{W}\_{\mathcal{P}}\beta\_{\mathcal{W}} - \mathcal{W}\_{\mathcal{i}} \tag{2}$$

where *Gp* and *Wp* are the cumulative gas and water production respectively, while *Wi* is the injection volume. In Equation (1), all production and injection volumes are in cubic meters unit. *Bg* and *Bw* are the volume factors of gas and water respectively, in (m3/m3) unit. Moreover, the expansion term can be stated as:

Expansion = volume change of inorganic matrix + volume change of organic matrix + desorbed gas volume + pore volume occupied by free gas—produced water. (3)

> All the volume changes included in the expansion term in Equation (3) are presented in the next sections (Sections 2.1–2.3).

#### *2.1. The Inorganic Matrix Volume Changes*

As the pressure in shale gas reservoirs decreases, there is a change in the volume of water and the pore spaces within the inorganic matrix. This change can be attributed to two factors: elastic variation in pore space within the reservoir, as well as expansion of formation water. The combined effect of these changes accounts for alterations observed in the inorganic matrix volume.

First, the variation in the inorganic matrix pore volume is equal to:

$$x\_p = \frac{1}{v\_p} \frac{\Delta v\_p}{\Delta p} \tag{4}$$

where *cp* is the pore space compressibility in MPa−<sup>1</sup> and *vp* is the pore volume in m3.

$$c\_w = \frac{1}{v\_w} \frac{\Delta v\_w}{\Delta p} \tag{5}$$

Similar to Equation (4) *cw* and *vw* are the water compressibility volumes, respectively. Re-arranging Equations (4) and (5) will give the change in reservoir pore and water volumes as follows:

$$
\Delta v\_p = c\_p \,\,\Delta p \,\, v\_p \tag{6}
$$

$$
\Delta v\_w = c\_w \Delta p \ v\_w \tag{7}
$$

Adding Equations (6) and (7), will present the total inorganic matrix volume change due to expansion as stated in Equation (8):

$$
\Delta v\_{exp} = c\_w \Delta p \, v\_w + c\_p \, \Delta p \, v\_p \tag{8}
$$

Using the definition of water saturation (*sw*), the water volume can be expressed in terms of pore volume, as follows:

$$
\upsilon\_w = s\_w \upsilon\_p \tag{9}
$$

Thus, Equation (8) can be re-written as follows:

$$
\Delta v\_{exp} = c\_w \Delta p \ s\_w v\_p + c\_p \Delta p \ v\_p = \Delta p \ v\_p \ (c\_w \ s\_w + c\_p \ )\tag{10}
$$

Moreover, the pore volume (*vp*) can be written in terms of the inorganic matrix free gas in place (*Ginorg*) as follows:

$$v\_p = \frac{v\_{\mathcal{g}}}{s\_{\mathcal{g}i}} = \frac{G\_{\text{inorg }\mathcal{B}\_{\text{\mathcal{g}i}}}}{1 - s\_{wi}} \tag{11}$$

where (*Ginorg*) is the inorganic matrix free gas in-place (m3) and *Bgi* is the initial gas volume factor, m3/m3. Therefore, the total expansion of pore space and formation is equal to:

$$
\Delta v\_{\rm exp} = c\_{\rm w} \Delta p \ s\_{\rm w} v\_p + c\_{\rm p} \Delta p \ v\_p = \Delta p \, \frac{G\_{\rm imag} \, \mathcal{J}\_{\rm gi}}{1 - s\_{\rm wi}} \left( c\_{\rm w} \, s\_{\rm w} + c\_{\rm p} \right) \tag{12}
$$

In addition, the free gas expansion in the inorganic matrix pores can be stated in terms of the inorganic matrix free gas and gas volume factor as follows:

$$G\_{inv\,\%} \left(\boldsymbol{\beta}\_{\mathcal{S}} - \boldsymbol{\beta}\_{\mathcal{S}^i}\right) \tag{13}$$

where *Bgi* is the initial gas volume factor and *Bg* is the gas volume factor at any pressure P and both are measured in (m3/m3) units.

#### *2.2. The Organic Matrix Volume Changes*

Shale deformation can be characterized by considering the impact of compression and matrix shrinkage, as explained by Plamer et al. (1998) [41], as follows:

$$\frac{\phi\_{sh}}{\phi\_{sh}} = 1 + \frac{c\_m}{\phi\_{sh}}(P - P\_i) + \frac{\xi\_l}{\phi\_{sh}} \left(\frac{K}{M} - 1\right) \left(\frac{P}{P\_L + P} - \frac{P\_i}{P\_L + P\_i}\right) \tag{14}$$

where *Φsh* and *Φish* are the initial shale porosity at the initial reservoir pressure and the shale porosity at any reservoir pressure, both dimensionless. *Cm* is the expansion factor of organic matter, *K* is the bulk modulus (MPa)*, ξ<sup>l</sup>* is the volumetric strain change (dimensionless), *M* is the axial strain modulus (MPa), *PL* is Langmuir's pressure (MPa), and *VL* is the maximum amount of adsorbed gas (Langmuir's volume, m3).

$$
\omega\_{\rm ff} = \frac{1}{M} - \left(\frac{K}{M} + 0.5\right)\gamma \tag{15}
$$

Here, *γ* is the solid compressibility of shale; *K/M* in Equation (15) is a function of Poisson's ratio (*ν*) as follows:

$$\frac{K}{M} = \frac{1}{3} \left( \frac{1+\upsilon}{1-\upsilon} \right) \tag{16}$$

The solid compressibility of shale (*γ*), can be neglected, so Equation (15) is re-written as:

$$\mathbb{C}\_{\text{ll}} = \frac{1}{M} \tag{17}$$

Equation (14) could be further simplified as follows:

$$\mathcal{Q}\_{\rm sh} - \mathcal{Q}\_{\rm shi} = \frac{1}{M} (P - P\_i) + \xi\_I \left( \frac{1}{3} \left( \frac{1+\upsilon}{1-\upsilon} \right) - 1 \right) \left( \frac{P}{P\_L + P} - \frac{P\_i}{P\_L + P\_i} \right) \tag{18}$$

Therefore, the organic matrix volume change is written as:

$$
\Delta V\_{\rm sh} = \frac{G\_{\rm forg} \mathcal{\beta}\_{\rm gi}}{\mathcal{Q}\_{\rm sh}} \left\{ \frac{1}{M} (P - P\_i) + \mathfrak{f}\_l \left( \frac{1}{3} \left( \frac{1+\upsilon}{1-\upsilon} \right) - 1 \right) \left( \frac{P}{P\_L + P} - \frac{P\_i}{P\_L + P\_i} \right) \right\} \tag{19}
$$

where *Gforg* is the organic matrix free gas volume in m3 and *Bgi* is the initial gas volume factor.

#### *2.3. Volume Change of the Adsorbed Shale Gas*

The desorbed gas volume in the organic matrix can be calculated as follows:

$$V\_{\rm des} = V\_L V\_b \rho\_{\rm slv} \frac{P}{P + P\_L} \tag{20}$$

where *VL* and *PL* are Langmuir's volume and pressure respectively; *P* is the reservoir pressure *Vb* is the reservoir bulk volume (m3), *ρsh* is the shale density.

Using Equation (20), the gas desorption rate into the matrix pore space can be written as follows:

$$-\frac{\partial V\_{des}}{\partial t} = -\left. V\_L V\_b \rho\_{sh} \frac{1}{\left(P + P\_L\right)^2} \frac{\partial P}{\partial t} \right. \tag{21}$$

The volume of desorbed gas is written in term of the organic matrix free gas volume as follows:

$$V\_{dcs} = \frac{G\_{fory} \beta\_{gi}}{\phi\_{sh}} \rho\_{sh} \beta\_{\mathcal{S}} \left(\frac{V\_L P\_i}{P\_i + P\_L} - \frac{V\_L P}{P + P\_L}\right) \tag{22}$$

#### *2.4. The Shale Gas Material Balance Equation*

The general form of the shale gas material balance equation can be written by summing all the derived volume changes in organic and non-organic matrix (i.e.; the sum of Equations (12), (13), (19) and (22)), as follows:

*Gpβ<sup>g</sup>* + *Wpβ<sup>w</sup>* − *wi*

$$\begin{aligned} \mathcal{I} &= \mathsf{G}\_{\text{inorg}} \Big[ \begin{pmatrix} \beta\_{\mathcal{S}} - \beta\_{\mathcal{g}i} \end{pmatrix} + \Delta p \frac{\beta\_{\mathcal{S}i}}{1 - \mathsf{S}\_{\text{in}}} \Big( c\_{w} s\_{w} + c\_{p} \Big) \Big] \\ + \mathsf{G}\_{f \text{or}} & \begin{pmatrix} \beta\_{\mathcal{S}} - \beta\_{\mathcal{g}i} \end{pmatrix} - \frac{\beta\_{\mathcal{g}i}}{\phi\_{\mathcal{s}h}} \left\{ \begin{array}{l} \frac{1}{M} (P - P\_{i}) + \\ \frac{\beta\_{\mathcal{I}}}{\phi\_{\mathcal{I}}} \left( \frac{1}{3} \left( \frac{1+\upsilon}{1-\upsilon} \right) - 1 \right) \left( \frac{P}{P\_{\mathcal{L}} + P} - \frac{P\_{i}}{P\_{\mathcal{L}} + P\_{i}} \right) \end{pmatrix} \right\} \\ + \frac{\beta\_{\mathcal{S}i}}{\phi\_{\mathcal{A}h}} \rho\_{\mathcal{S}h} \beta\_{\mathcal{S}} \left( \frac{V\_{\mathcal{I}} P\_{i}}{P\_{\mathcal{I}} + P\_{\mathcal{L}}} - \frac{V\_{\mathcal{I}} P}{P + P\_{\mathcal{L}}} \right) \end{aligned} \tag{23}$$

To simplify the presentation of Equation (23), the new parameters (*F*, *Eg*, and *X*) are introduced. The parameter *F* represents the production term in Equation (23). The inorganic matrix expansion term in Equation (23) is represented by parameter *Eg*, while the organic matrix expansion term is represented by the parameter *X*. The three parameters are expressed as follows:

$$F = G\_p \beta\_\mathcal{S} + \mathcal{W}\_p \beta\_{\mathcal{W}} - \mathcal{W}\_i \tag{24}$$

$$E\_{\mathcal{S}} = \left[ \left( \beta\_{\mathcal{S}} - \beta\_{\mathcal{S}^i} \right) + \Delta p \, \frac{\beta\_{\mathcal{S}^i}}{1 - s\_{\text{tvi}}} \left( c\_{\text{uv}} s\_{\text{uv}} + c\_p \right) \right] \tag{25}$$

$$X = \left\{ \begin{array}{l} \left( \beta\_{\mathcal{S}} - \beta\_{\mathcal{S}^{\bar{i}}} \right) - \frac{\beta\_{\mathcal{S}^{\bar{i}}}}{\mathcal{Q}\_{\text{sh}}} \left\{ \frac{1}{\mathcal{M}} (P - P\_{\bar{i}}) + \tilde{\xi}\_{\text{I}} \left( \frac{1}{\mathfrak{Z}} \left( \frac{1 + \upsilon}{1 - \upsilon} \right) - 1 \right) \left( \frac{P}{P\_{\mathcal{L}} + P} - \frac{P\_{\bar{i}}}{P\_{\mathcal{L}} + P\_{\bar{i}}} \right) \right\} \right\} \\ \qquad + \frac{\beta\_{\mathcal{S}^{\bar{i}}}}{\mathcal{Q}\_{\text{sh}}} \rho\_{\text{sh}} \beta\_{\mathcal{S}} \left( \frac{V\_{\mathcal{I}} P\_{\mathcal{I}}}{P\_{\mathcal{I}} + P\_{\mathcal{L}}} - \frac{V\_{\mathcal{I}} P}{P + P\_{\mathcal{L}}} \right) \end{array} \tag{26}$$

The parameter *X* can be written in a shorter form as follows:

$$X = \{\mathbf{x}\_1 - \mathbf{x}\_2 + \mathbf{x}\_3\} \tag{27}$$

where *x*<sup>1</sup> is the organic matrix volume change due to free gas expansion (the first term in Equation (26)), *x*<sup>2</sup> is the organic matrix volume change due to deformation from matrix shrinkage and stress sensitivity (the second term in Equation (26)), and *x*<sup>3</sup> is the organic matrix volume change due to gas desorption process (the third term in Equation (26)).

In conclusion, Equation (23) can be re-arranged in terms of the new parameters (*F*, *Eg*, and *X*) as follow:

$$\frac{F}{E\_{\mathcal{S}}} = G\_{\text{inorg}} + G\_{\text{forg}} \frac{X}{E\_{\mathcal{S}}} \tag{28}$$

Equation (28) is a linear relationship between the production and expansion terms of the shale gas material balance equation; where the slope of the straight line in Equation (28) presents the free gas volume in the organic matrix *Gforg*, and the intercept is the free gas volume of the inorganic matrix *Ginorg*. Both volumes can be determined through linear fitting using the plot of *F/Eg* vs. *X/Eg*.

Finally, the volume of adsorbed gas (*GAbsd*) as a function of the free gas volume in the organic matrix *Gforg* and the Langmuir's parameters is as follow:

$$G\_{Absd} = \frac{G\_{for\text{\textdegree}} \beta\_{\text{\textdegree}i}}{\phi\_{sli}} \rho\_{sl} \left(\frac{V\_L P\_i}{P\_i + P\_L}\right) \tag{29}$$

where, *VL* is the Langmuir's volume and *PL* is the Langmuir's pressure; ρsh and ϕsh are the shale density and porosity, respectively.

The calculation procedures of the proposed MBE are summarized as follows:


In the following section, a brief description of the in-house poro-elastic finite element simulator which will be used to validate the proposed material balance equation is given.

#### **3. Poro-Elastic In-House Numerical Simulation Model**

The flow of fluid in shale gas reservoirs was modeled in this study using an in-house poro-elastic simulator. The finite element method was utilized to discretize the equations that govern fluid flow in two-dimensional space.

The 2-D continuity equation of gas flow is written as follows:

$$\frac{\partial}{\partial x} \left( \rho\_{\mathcal{S}} \mu\_{\mathcal{S}} \right) - a \frac{\partial}{\partial t} \left( \frac{\partial}{\partial x} \frac{\partial}{\partial y} \right) \vec{\mu} = -\frac{\partial}{\partial t} \left( \phi \rho\_{\mathcal{S}} \right) \tag{30}$$

where *ρ<sup>g</sup> and ug* are the gas density and velocity, respectively.

By incorporating the gas formation volume factor and the source/sink term (*qg*) into Equation (30), we obtain the following equation:

$$\frac{\partial}{\partial \boldsymbol{\alpha}} \left( \beta\_{\mathcal{S}} \boldsymbol{u}\_{\mathcal{S}} \right) - \boldsymbol{\alpha} \frac{\partial}{\partial t} \left( \frac{\partial}{\partial \boldsymbol{x}} \frac{\partial}{\partial \boldsymbol{y}} \right) - q\_{\mathcal{S}} = -\frac{\partial}{\partial t} \left( \phi \beta\_{\mathcal{S}} \right) \tag{31}$$

Assuming off-diagonal components of the permeability tensor are zero →→ *k* = *kx* 0 0 *ky* , Equation (31) can be reformulated as Darcy's velocity and utilized as follows:

$$\frac{\partial}{\partial \mathbf{x}} \left( \frac{ck\_{\mathbf{x}} \boldsymbol{\beta}\_{\mathcal{S}}}{\mu\_{\mathcal{S}}} \frac{\partial p\_{\mathcal{S}}}{\partial \mathbf{x}} \right) + q\_{\mathcal{S}} = \boldsymbol{\phi} \frac{\partial p\_{\mathcal{S}}}{\partial t} \frac{\partial \boldsymbol{\beta}\_{\mathcal{S}}}{\partial p\_{\mathcal{S}}} \tag{32}$$

Writing Equation (32) in 2D form will give:

$$\frac{\partial}{\partial x}\left(\frac{ck\_{\mathfrak{x}}\beta\_{\mathfrak{S}}}{\mu\_{\mathfrak{S}}}\frac{\partial p\_{\mathfrak{S}}}{\partial x}\right) + \frac{\partial}{\partial z}\left(\frac{ck\_{\mathfrak{x}}\beta\_{\mathfrak{S}}}{\mu\_{\mathfrak{S}}}\frac{\partial p\_{\mathfrak{S}}}{\partial z}\right) + q\_{\mathfrak{S}} = \phi \frac{\partial p\_{\mathfrak{S}}}{\partial t} \frac{\partial \beta\_{\mathfrak{S}}}{\partial p\_{\mathfrak{S}}}\tag{33}$$

Modifying Equation (33) to include the amount of adsorbed gas as a source term (i.e., injection well) will give:

$$\frac{\partial}{\partial x}\left(\frac{ck\_{\mathfrak{f}}\beta\_{\mathfrak{f}}}{\mu\_{\mathfrak{f}}}\frac{\partial p\_{\mathfrak{f}}}{\partial x}\right) + \frac{\partial}{\partial y}\left(\frac{ck\_{\mathfrak{f}}\beta\_{\mathfrak{f}}}{\mu\_{\mathfrak{f}}}\frac{\partial p\_{\mathfrak{f}}}{\partial y}\right) + V\_{\mathbb{L}}\rho\_{\mathbb{R}}\frac{1}{(p+p\_{\mathfrak{f}})^{2}}\frac{\partial p}{\partial t} - a\frac{\partial}{\partial t}\left(\frac{\partial}{\partial x}\frac{\partial}{\partial y}\right)\overset{\rightarrow}{\mathcal{U}} = \phi\frac{\partial p\_{\mathfrak{f}}}{\partial t}\frac{\partial \beta\_{\mathfrak{f}}}{\partial p\_{\mathfrak{f}}}\tag{34}$$

Multiplying both sides of Equation (34) by a trial function *w* and integrating over the domain, Ω yields:

$$\begin{split} \int\_{\Omega} \left( w c\_t \boldsymbol{\uprho} \frac{\partial p\_\mathcal{E}}{\partial t} \frac{\partial \beta\_\mathcal{E}}{\partial p\_\mathcal{E}} \right) d\Omega & - V\_\mathcal{L} \rho\_\mathcal{R} \frac{1}{(p + p\_\mathcal{E})^2} \int\_{\Omega} \left( w \frac{\partial p\_\mathcal{E}}{\partial t} \right) d\Omega - \\ \int\_{\Omega} \left[ w u \boldsymbol{\uprho} \frac{\partial}{\partial t} \left( \frac{\partial}{\partial x} \frac{\partial}{\partial y} \right) \overrightarrow{u} \right] d\Omega & = \int\_{\Omega} w \left( \begin{array}{l} \frac{\partial}{\partial x} \left( \frac{c k\_\mathcal{E} \beta\_\mathcal{E}}{\mu\_\mathcal{E}} \frac{\partial p\_\mathcal{E}}{\partial x} \right) + \\ \frac{\partial}{\partial y} \left( \frac{c k\_\mathcal{E} \beta\_\mathcal{E}}{\mu\_\mathcal{E}} \frac{\partial p\_\mathcal{E}}{\partial y} \right) \end{array} \right) d\Omega \end{split} \tag{35}$$

Applying Green Formula, Equation (35) is written as:

$$\begin{array}{l} \int\_{\Omega} \left( wc\_{t} \rho \frac{\partial p\_{\mathcal{S}}}{\partial t} \frac{\partial \beta\_{\mathcal{S}}}{\partial p\_{\mathcal{S}}} \right) d\Omega - V\_{L} \rho\_{R} \frac{1}{(p+p\_{\mathcal{S}})^{2}} \int\_{\Omega} \left( w \frac{\partial p\_{\mathcal{S}}}{\partial t} \right) d\Omega = \int\_{\Omega} w \frac{\partial w}{\partial x} \frac{\partial}{\partial x} \left( \frac{ck\_{g} \beta\_{\mathcal{S}}}{\mu\_{\mathcal{S}}} \frac{\partial p\_{\mathcal{S}}}{\partial y} \right) d\Omega\\ \int\_{\Omega} \frac{\partial w}{\partial x} \left( \frac{ck\_{g} \beta\_{\mathcal{S}}}{\mu\_{\mathcal{S}}} \frac{\partial p\_{\mathcal{S}}}{\partial x} \right) d\Omega + \frac{\partial w}{\partial y} \left( \frac{ck\_{g} \beta\_{\mathcal{S}}}{\mu\_{\mathcal{S}}} \frac{\partial p\_{\mathcal{S}}}{\partial y} \right) d\Omega \\ - \int\_{\Omega} \oint\_{\Gamma} w \left( \frac{ck\_{g} \beta\_{\mathcal{S}}}{\mu\_{\mathcal{S}}} \frac{\partial p\_{\mathcal{S}}}{\partial x} n\_{\mathcal{S}} + \frac{ck\_{g} \beta\_{\mathcal{S}}}{\mu\_{\mathcal{S}}} \frac{\partial p\_{\mathcal{S}}}{\partial y} n\_{\mathcal{S}} \right) d\Gamma + \int\_{\Omega} \left[ w \alpha \frac{\partial}{\partial t} \left( \frac{\partial}{\partial x} \frac{\partial}{\partial y} \right) \stackrel{\rightarrow}{\mu} \right] d\Omega \end{array} \tag{36}$$

where (Γ) is the boundary and (*n*) is the outward normal to the boundary. The finite difference method is used to discretize the terms including differentiation with respect to time. Permeability and porosity remain constant with changes in time. Re-arranging Equation (36) gives:

$$\begin{split} \int\_{\Omega} \left( w c\_{t}^{i-1} \Phi^{j-1} \frac{p^{i} - p^{i-1}}{\Delta t^{i}} \frac{p^{i} - p^{i-1}}{\beta\_{\mathcal{S}} - \beta\_{\mathcal{S}}{t}^{i-1}} \right) d\Omega &+ V\_{L} \rho\_{\mathcal{R}} \frac{1}{(p + p\_{L})^{2}} \int\_{\Omega} \left( w \frac{p^{i} - p^{i-1}}{\Delta t^{i}} \right) d\Omega \\ &+ \int\_{\Omega} \frac{\text{div}}{\partial x} \left( \frac{c k\_{x}^{i-1} \beta\_{\mathcal{R}}}{\mu\_{\mathcal{S}}} \frac{\partial p^{i}}{\partial x} \right) d\Omega + \\ \int\_{\Omega} \frac{\text{div}}{\partial y} \left( \frac{c k\_{y}^{i-1} \beta\_{\mathcal{R}}}{\mu\_{\mathcal{S}}} \frac{\partial p^{i}}{\partial y} \right) d\Omega &= \int\_{\Omega} \left[ u u \left( \frac{\partial}{\partial x} \frac{\partial}{\partial y} \right) \frac{\overset{\rightarrow i}{u} - u^{i-1}}{\Delta t^{i}} \right] d\Omega \end{split} (37)$$

In which (*i*) and (*i* − 1) are current and previous times respectively. Using Galerkin approach [42], Equation (37) is re-written as follow:

$$
\begin{split}
\left[\int\_{\Omega} c\_{l} \boldsymbol{\phi}^{j} \overset{\rightarrow}{N}\_{p} \overset{\rightarrow}{N}\_{p} d\Omega \right] \frac{\overset{\rightarrow}{p}^{i} - \stackrel{\rightarrow}{P}^{i-1}}{\Delta t^{i}} \frac{\overset{\rightarrow}{p}^{i} - \stackrel{\rightarrow}{P}^{i-1}}{\hat{\mu}\_{\text{x}}^{i} - \hat{\mu}\_{\text{x}}^{i-1}} + V\_{l} \rho\_{R} \frac{1}{(p+p\_{L})} \begin{bmatrix} \int\_{\Omega} \stackrel{\rightarrow}{N}\_{p} \overset{\rightarrow}{N}\_{p} d\Omega \end{bmatrix} \frac{\overset{\rightarrow}{p}^{i} - \stackrel{\rightarrow}{P}^{i-1}}{\Delta t^{i}} \\ + \left[ \int\_{\Omega} \left( \frac{\mathrm{c} k\_{x}^{i-1} \hat{\mu}\_{x}}{\mu\_{\text{g}}^{i-1}} \frac{\overset{\rightarrow}{\partial \hat{N}\_{x}}^{\hat{\mu}}}{\mathrm{d}x} + c \oint\_{\mathcal{S}} \frac{k\_{y}^{i-1}}{\mu\_{\text{x}}^{i-1}} \frac{\overset{\rightarrow}{\partial \hat{N}\_{y}}^{\hat{\mu}}}{\mathrm{d}y} \frac{\overset{\rightarrow}{\partial \hat{N}\_{x}}}{\mathrm{d}y} \right) d\Omega \right] \overset{\rightarrow}{P}^{i} + \left( \int\_{\Omega} \stackrel{\rightarrow}{N}\_{p} \overset{\rightarrow}{a} \left( \frac{\partial}{\partial x} \frac{\partial}{\partial y} \right) d\Omega \right) \frac{\overset{\rightarrow}{\partial \hat{I}} - \hat{\mu}^{i}}{\Delta t^{i}} = 0
\end{split} \tag{38}
$$

where

$$\begin{array}{ll} \stackrel{\rightarrow}{P}^{T} &= (p\_1 p\_2 \cdots \cdot p\_n) \\ \stackrel{\rightarrow}{N}\_p^{T} &= (N\_1 N\_2 \cdots N n) \\ \stackrel{\rightarrow}{N\_u} &= \begin{bmatrix} N\_1 0 N\_2 \cdots \cdot 0 \\ 0 N\_1 0 \cdots N\_n \end{bmatrix} \\ \stackrel{\rightarrow}{U}^{T} &= \begin{pmatrix} \mu\_{x1} \mu\_{y1} \mu\_{x2} \cdots \mu\_{yu} \end{pmatrix} \end{array} \tag{39}$$

(*n*) is the total number of nodes.

Equation (37) can be further re-arranged as follow:

$$
\stackrel{\rightarrow}{L} \left( \stackrel{\rightarrow}{p}^i - p^{i-1} \right) + \Delta t^i \stackrel{\rightarrow}{H} \stackrel{\rightarrow}{p}^i - \stackrel{\rightarrow}{Q} \left( \stackrel{\rightarrow}{\mathcal{U}}^i - \mathcal{U}^{i-1} \right) = \stackrel{\rightarrow}{0} \tag{40}
$$

where →→ *<sup>L</sup>* is as follow: →→

$$\stackrel{\rightarrow}{L} = \sum\_{c=1}^{nc} \stackrel{\rightarrow}{L^c} \tag{41}$$

Here, (*e*) denotes to element and (*ne*) is the number of elements.

→→

$$\stackrel{\rightarrow}{L^{\varepsilon}} = \int\_{\Omega} \left( \frac{P^{i} - P^{i-1}}{\beta\_{\mathcal{S}}^{\cdot} - \beta\_{\mathcal{S}}^{\cdot -1}} c\_{l}^{j-1} \phi^{\varepsilon^{i-1}} \stackrel{\rightarrow}{N\_{p}}^{T} \stackrel{\rightarrow}{N\_{p}} \right) d\Omega + V\_{L} \rho\_{R} \frac{1}{(p + p\_{L})^{2}} \int\_{\Omega} \left( \stackrel{\rightarrow}{N\_{p}}^{T} \stackrel{\rightarrow}{N\_{p}}\_{p} \right) \tag{42}$$

In which

$$\stackrel{\rightarrow}{H} = \stackrel{\leftarrow}{\sum} \stackrel{\rightarrow}{H^c} \tag{43}$$

$$\stackrel{\rightarrow}{H^{\xi}} = \int\_{\Omega} e \left( c \mathcal{J}\_{\mathcal{S}} \frac{k\_{\mathcal{X}}^{\epsilon^{i-1}}}{\mu^{\epsilon^{i-1}}} \frac{\stackrel{\rightarrow}{\partial N\_p^{\epsilon}}}{\delta x} \frac{\stackrel{\rightarrow}{\partial N\_p^{\epsilon}}}{\delta x^{\epsilon^{i-1}}} + c \mathcal{J}\_{\mathcal{S}} \frac{k\_{\mathcal{Y}}^{\epsilon^{i-1}}}{\mu^{\epsilon^{i-1}}} \frac{\stackrel{\rightarrow}{\partial N\_p^{\epsilon}}}{\delta y} \frac{\stackrel{\rightarrow}{\partial N\_p^{\epsilon}}}{\delta y} \right) d\Omega \tag{44}$$
 
$$\stackrel{\rightarrow}{Q} = \sum\_{\epsilon=1}^{nc} \stackrel{\rightarrow}{Q}^{\epsilon}$$

where

$$Qe = \alpha \int\_{\Omega} e \overset{\rightarrow}{N}\_p^T \left( \frac{\partial N^\varepsilon u}{\partial x} + \frac{\partial N^\varepsilon u}{\partial x} \right) d\Omega \tag{45}$$

Equation (40) can be written in an incremental form as follows:

$$
\begin{array}{ccc}
\stackrel{\rightarrow}{U} & \stackrel{\rightarrow}{\rightarrow} & \stackrel{\rightarrow}{\rightarrow} & \stackrel{\rightarrow}{\rightarrow} \\
\stackrel{\rightarrow}{L} & \stackrel{\rightarrow}{\rightarrow} & \stackrel{\rightarrow}{\rightarrow} & \stackrel{\rightarrow}{\rightarrow} \\
\stackrel{\rightarrow}{\Delta P} = P^{\dot{i}} - P^{\dot{i}-1}, \stackrel{\rightarrow}{\Delta U} = \mathcal{U}^{\dot{i}} - \mathcal{U}^{\dot{i}-1} \\
& \stackrel{\rightarrow}{f}\_{1} = -\Delta t^{\dot{i}} H P^{\dot{i}-1}
\end{array} \tag{46}
$$

$$
\stackrel{\rightarrow}{N}\_{\mathcal{P}}^{T} = (N\_1 N\_2 \cdots N n)
$$

$$
\stackrel{\rightarrow}{N\_{\mathcal{U}}} = \begin{bmatrix} N\_1 & 0 & N\_1 & \cdots & 0 \\ 0 & N\_1 & 0 & \cdots & N\_1 \end{bmatrix} \tag{47}
$$

→ *Np T* is the pressure shape function, *<sup>n</sup>* is the number of nodes, <sup>→</sup> *Nu* is the displacement function, and P represents the pressure nodal values.

Gaussian quadrature is applied for integration. Also, for the coupled equilibrium equation multiplying Equation (46) by trial function and integrating over the domain Ω leads to:

$$\begin{aligned} \int\_{\Omega} \begin{pmatrix} \stackrel{\rightarrow}{W} \stackrel{\rightarrow}{S^T} \stackrel{\rightarrow}{\Delta \sigma} \\ \stackrel{\rightarrow}{W} \stackrel{\rightarrow}{S^T} \end{pmatrix} d\Omega &= \stackrel{\rightarrow}{0} \\ \stackrel{\rightarrow}{W} = \begin{bmatrix} \,^{W\_1(\mathbf{x}, \mathbf{y})} \\ \,^{W\_2(\mathbf{x}, \mathbf{y})} \end{bmatrix} \end{aligned} \tag{48}$$

and *W*<sup>1</sup> and *W*<sup>2</sup> are trial functions. Using Green's identity, one obtains:

$$\begin{aligned} \int\_{\Gamma} \begin{pmatrix} \stackrel{\rightarrow}{W} \stackrel{\rightarrow}{M} \stackrel{\rightarrow}{M} \stackrel{\rightarrow}{\partial \sigma} \\ \end{pmatrix} d\Gamma - \int\_{\Omega} \stackrel{\rightarrow}{\Delta \sigma}^{T} \begin{pmatrix} \stackrel{\rightarrow}{\partial \dot{W}} \\ \dot{\bar{S}} \dot{W} \end{pmatrix} d\Omega &= \stackrel{\rightarrow}{0} \\ \stackrel{\rightarrow}{M} \stackrel{\rightarrow}{\partial \dot{T}} &= \begin{bmatrix} n\_{x} \ 0 \ n\_{y} \\ 0 \ n\_{y} n\_{x} \end{bmatrix} \end{aligned} \tag{49}$$

substituting Equation (49) into Equation (48) and rearranging, we will obtain:

$$\int\_{\Omega} \begin{pmatrix} \stackrel{\rightarrow}{\rightarrow} & \stackrel{\rightarrow}{\rightarrow} & i \\ \stackrel{\rightarrow}{\mathbf{S}} \stackrel{\rightarrow}{\mathbf{W}} & \stackrel{\rightarrow}{\mathbf{S}} \stackrel{\rightarrow}{\mathbf{M}} \stackrel{i}{\mathbf{I}} \, d\Omega + a \int\_{\Omega} \begin{pmatrix} \stackrel{\rightarrow}{\mathbf{S}} \stackrel{\rightarrow}{\mathbf{W}} \\ \stackrel{\rightarrow}{\mathbf{S}} \stackrel{\rightarrow}{\mathbf{W}} \end{pmatrix} \stackrel{i}{\Delta P} \, d\Omega = \int\_{\Gamma} \begin{pmatrix} \stackrel{\rightarrow}{\mathbf{W}} \stackrel{\rightarrow}{\mathbf{M}} \stackrel{i}{\mathbf{I}} \stackrel{\rightarrow}{\mathbf{}} \\ \end{pmatrix} d\Gamma \end{pmatrix} \tag{50}$$

Using Galerkin method Equation (50) is re-written as follows:

$$\begin{aligned} &a\left(\int\_{\Omega} \left(\left(\frac{\stackrel{\scriptstyle\cdot}{\partial\boldsymbol{N}}\_{\boldsymbol{\mu}}\stackrel{\scriptstyle\cdot}{\partial\boldsymbol{N}}\_{\boldsymbol{\mu}}}{\partial\boldsymbol{X}}\right)^{\scriptstyle T}\_{\boldsymbol{P}}\right)^{T} d\Omega\right) \stackrel{\scriptstyle\cdot}{\Delta} + V\_{L} \rho\_{R} \frac{1}{(p+p\_{\boldsymbol{L}})^{\mathcal{L}}} \int\_{\Omega} \left(\int\_{\Omega} \left(\left(\frac{\stackrel{\scriptstyle\cdot}{\partial\boldsymbol{N}}\_{\boldsymbol{\mu}}\stackrel{\scriptstyle\cdot}{\partial\boldsymbol{N}}\_{\boldsymbol{\mu}}}{\partial\boldsymbol{X}}\right)^{\scriptstyle T}\_{\boldsymbol{P}}\right)^{T} d\Omega\right)^{i} \\ &= \left(\int\_{\Omega} \left(\stackrel{\scriptstyle\cdot}{\partial\boldsymbol{N}}\_{\boldsymbol{\mu}}\right)^{\scriptstyle T} \stackrel{\scriptstyle\cdot}{\partial\boldsymbol{N}}\_{\boldsymbol{\mu}} \stackrel{\scriptstyle\cdot}{\partial\boldsymbol{N}}\_{\boldsymbol{\mu}} d\Omega\right) \stackrel{\scriptstyle\cdot}{\Delta}{\boldsymbol{M}} \end{aligned} \tag{51}$$

Equation (51) in compact form is as follow:

$$
\stackrel{\rightarrow}{K} \stackrel{\rightarrow}{\Delta U} + \stackrel{\rightarrow}{Q}^{\rightarrow} \stackrel{\rightarrow}{\Delta P} = \stackrel{\rightarrow}{f\_2} \tag{52}
$$

Finally, Equations (53) to (57) are the final finite element equations simultaneously solved as a system of linear equations as follows:

$$\stackrel{\rightarrow}{K} = \sum\_{\mathfrak{c}=1}^{\mathfrak{nc}} \stackrel{\rightarrow}{K^{\mathfrak{c}}} \tag{53}$$

$$\stackrel{\rightarrow}{K^{\varepsilon}} = \int\_{\Omega^{\varepsilon}} \left( \stackrel{\rightarrow}{S} \stackrel{\rightarrow}{N\_{u}} \stackrel{\rightarrow}{\iota} \right)^{T} \stackrel{\rightarrow}{D} \stackrel{\rightarrow}{S} \stackrel{\leftarrow}{N\_{u}} \stackrel{\leftarrow}{d} \Omega \tag{54}$$

$$\stackrel{\rightarrow}{\vec{f}}\_2 = \sum\_{c=1}^{nc} \stackrel{\rightarrow}{\vec{f}}\_2^c \tag{55}$$

$$\stackrel{\rightarrow}{f}\_2 = \int\_{\Gamma^\varepsilon} \left( \stackrel{\rightarrow}{N}\_u \stackrel{\rightarrow}{M}^T \stackrel{\rightarrow}{\Delta \sigma}^i \right) d\Gamma \tag{56}$$

$$
\begin{bmatrix}
\stackrel{\rightarrow}{\rightarrow} & \stackrel{\rightarrow}{\rightarrow} & \stackrel{\rightarrow}{\rightarrow} \\
L + \Delta t H - Q \\
\stackrel{\rightarrow}{\rightarrow} & \stackrel{\rightarrow}{\rightarrow} \\
Q^T & K \end{bmatrix} \begin{bmatrix}
\Delta \stackrel{\rightarrow}{P} \\
\Delta U
\end{bmatrix} = \begin{bmatrix}
\stackrel{\rightarrow}{f\_1} \\
\stackrel{\rightarrow}{f\_2} \\
\stackrel{\rightarrow}{f\_2}
\end{bmatrix} \tag{57}
$$

A finite element poro-elastic model is presented in this section which accounts for the deformation effects in shale gas reservoirs resulting from the combined matrix shrinkage, stress sensitivity, and gas desorption processes. The model is then used to validate the proposed material balance equation presented earlier in Section 1 of this study. Moreover, two published real-field cases are also used to further validate the material balance equation using real production data, as will be shown in the next section. Finally, it should be noted that a brief demonstration of the poro-elastic model's ability to capture stress and deformation changes around production wells is presented in Appendix A.

#### **4. Results and Discussion**

#### *4.1. Verification of the Modified MBE Using the In-House Numerical Poro-Elastic Model*

The in-house simulator was used to verify the proposed modified material balance equation by conducting several numerical cases to simulate the production of a shale gas reservoir under a poro-elastic environment that accounts for the coupled effects of stress sensitivity, matrix shrinkage, and critical desorption pressure. Figure 1 presents a demonstration of one of the in-house simulator cases used during the verification process; in this figure the in-house simulator was run for a shale gas reservoir with a drainage radius of 100 m and a matrix permeability of 0.001 md. Other rock and fluid properties used in this run are shown in Table 1. Figure 1A presents the pressure decline rate for two elements (element 1 and element 2): element 1 is close to the wellbore while the other is far from the wellbore. It can be seen from this figure that the mesh element which is close to the well vicinity goes through a faster rate of pressure decline while the far field element shows a slower decline rate due to the low matrix permeability. Moreover, Figure 1B presents the amount of desorbed gas volumes in element 1 and element 2 discussed in Figure 1A; it can be seen from Figure 1B that the element which is close to the wellbore vicinity (element 1) starts producing a desorbed gas volume at a rate faster than that seen in the element far from the wellbore due to the fact that the pressure decline rate is faster near the wellbore where the pressure will quickly reach the critical desorption pressure and the absorbed gas layer in the matrix wall will start to desorb gas into the matrix pore spaces and then to the wellbore. Moreover, the well was run to an abonnement pressure of 5 MPa; then, the proposed modified material balance was used to estimate the expansion/production parameters (i.e., the F, X, and Eg parameters from Equations 24,25, and 26) for the determination of total gas in place. Figure 2A presents the cumulative gas production versus the average reservoir pressure, while Figure 2B presents the change in the gas volume factor as a function of the average reservoir pressure; data in Figure 2A, B were used as an input to the MBE calculation. Moreover, Table 2 presents the calculation procedures of the proposed MBE, and Figure 3 presents the relationship between the production and reservoir expansion terms. Based on Figure 3, it is evident that a strong correlation exists between the production and reservoir expansion terms of the studied reservoir. The relationship between these two terms can be described as a linear relationship, with a coefficient of determination of 0.99 (R<sup>2</sup> = 0.99), a slope of (1.22 × 106), and an intercept of (0.16 × <sup>10</sup>4). Thus, according to Equation (28) the total free gas in place is equal to 1.2216 × 106 <sup>m</sup>3, which includes a free gas volume in shale (organic matrix) of 1.22 × 106 <sup>m</sup><sup>3</sup> and a volume of 0.16 × <sup>10</sup><sup>4</sup> <sup>m</sup><sup>3</sup> of free gas in the inorganic matrix. Furthermore, the total free gas in place from the in-house simulator was equal to 1.16 × <sup>10</sup><sup>6</sup> m3, hence the error between the total free gas volume from the modified MBE and that from the in-house poro-elastic numerical model is around 5%, which highlights the reliability of the presented MBE in calculating shale gas in place. In the following section, additional validation tests of the proposed MBE will be presented using real-field cases.

**Figure 1.** (**A**) Pressure decline rate in two elements (element-1 near-wellbore and element-2 far- from wellbore); and (**B**) desorbed gas volume in two elements (element-1 near-wellbore and element-2 farfrom wellbore), after one year of production.

**Table 1.** Reservoir and rock properties for the in-house numerical model case.


**Figure 2.** (**A**) Cumulative gas production versus average reservoir pressure; and (**B**) gas volume factor versus average reservoir pressure.

**Table 2.** MBE equation calculation procedures for the numerical in-house model run.


*<sup>x</sup>***<sup>1</sup>** = volume changes due to free gas expansion in Organic matrix, *<sup>x</sup>***<sup>2</sup>** = volume changes due to shale matrix shrinkage and stress sensitivity, *<sup>x</sup>***<sup>3</sup>** = volume changes due to desorbed gas volume. Note: *<sup>X</sup>* = is the summation of (x1, x2, and x3).

**Figure 3.** Linear relationship between the production and expansion parameters proposed in the modified material balance equation presented in this study (using the production data from the numerical case of the in-house model).

#### *4.2. Verification of the Modified MBE Using Real-Field Case Study*

To further validate the proposed modified material balance equation, two real-field cases of shale gas reservoirs located in China were used from the literature, Hu et al. 2019 [37]. Tables 3 and 4 presents the cumulative gas and water production versus the average reservoir pressure for two wells, Well-A and Well-B. Moreover, Table 5 presents the rock and fluid properties of the reservoir.

**Table 3.** Shale gas production history of Well-A. Reprinted with permission from reference [37] Copyright 2019 Elsevier.


**Table 4.** Shale gas production history of Well-B. Reprinted with permission from reference [37]. Copyright 2019 Elsevier.


**Table 5.** Reservoir and rock properties for the real-field case study. Reprinted with permission from reference [37]. Copyright 2019 Elsevier.


Table 6 presents the proposed MBE calculation procedure for Well-A, and Figure 4 presents the linear relationship between the production and expansion terms of Well-A dynamic production data. In Figure 4, the slope of the straight line (i.e., the free gas volume in the organic matrix) is equal to 1.57 × <sup>10</sup><sup>8</sup> <sup>m</sup><sup>3</sup> while the intercept (i.e., the free gas volume in the inorganic matter) is equal to 2.25 × <sup>10</sup><sup>5</sup> m3; moreover, the absorbed gas volume is equal to 4.56 × <sup>10</sup><sup>7</sup> m3 using equation 29. Accordingly, the sum of these three volumes will give a total gas in place volume of 2.03 × <sup>10</sup><sup>8</sup> m3 which is in close agreement with the total gas in place of 2.10 × 108, as reported by Hu et al. 2018 [37], (i.e., the difference between the reported total gas in place and the one calculated using the modified MBE presented herein is around 3%). Furthermore, it should be noted that if the matrix shrinkage and stress sensitivity effects as well as the change in volume due to gas adsorption are ignored (i.e., ignoring the values of columns (*x*<sup>2</sup> and *x*3) in Table 6) then the estimated total gas volume will be equal to 2.13 × 108 <sup>m</sup>3, which gives an over-estimation of 5% compared to the case where all effects are included.


**Table 6.** MBE equation calculation procedures for Well-A production data.

*<sup>x</sup>***<sup>1</sup>** = volume changes due to free gas expansion in Organic matrix, *<sup>x</sup>***<sup>2</sup>** = volume changes due to shale matrix shrinkage and stress sensitivity, *<sup>x</sup>***<sup>3</sup>** = volume changes due to desorbed gas volume. Note: *<sup>X</sup>* = is the summation of (x1, x2, and x3).

**Figure 4.** The linear relationship between the production and expansion parameters proposed in the modified material balance equation presented in this study (using the production data from the real-field production data of Well-A).

Table 7 presents the proposed MBE calculation procedure for Well-B. The total gas in place of Well-B is equal to 1.24 × <sup>10</sup><sup>8</sup> as reported by Hu et al. 2018 [37]. Figure <sup>5</sup> presents the linear relationship between the production and expansion terms of Well-B production data. In Figure 5, the slope of the straight line (i.e., the free gas volume in organic matrix) is equal to 1.28 × <sup>10</sup><sup>8</sup> <sup>m</sup>3, while the intercept (i.e., the free gas volume in the inorganic matter) is equal to 2.38 × <sup>10</sup><sup>4</sup> m3; moreover, the absorbed gas volume is equal to 3.72 × <sup>10</sup><sup>6</sup> m3 using equation 29. Therefore, the MBE calculated total gas in place volume of Well-B is equal to 1.32 × 108 <sup>m</sup><sup>3</sup> which is in close agreement with the total gas in place of 1.24 × <sup>10</sup><sup>8</sup> as reported by Hu et al. 2018 [37], (i.e., the difference between the reported total gas in place and the one calculated using the modified MBE presented herein is around 6%). Such encouraging results from the real-field data of Well-A and Well-B highlight the reliability of the presented modified material balance equation in estimating shale gas in place, accounting for the combined effects of stress change, matrix shrinkage, water production, and critical desorption pressure.

**Pressure (MPa)** *F x***<sup>1</sup>** *<sup>x</sup>***<sup>2</sup>** *<sup>x</sup>***<sup>3</sup>** *X Eg* 17.2 2.416 <sup>×</sup> <sup>10</sup><sup>5</sup> 0.00134 <sup>−</sup>1.9345 <sup>×</sup> <sup>10</sup>−<sup>4</sup> 8.30205 <sup>×</sup> <sup>10</sup>−<sup>8</sup> 0.001533529 0.710815 13.8 4.444 <sup>×</sup> <sup>10</sup><sup>5</sup> 0.00294 <sup>−</sup>3.1549 <sup>×</sup> <sup>10</sup>−<sup>4</sup> 1.74589 <sup>×</sup> <sup>10</sup>−<sup>7</sup> 0.003255668 0.994941 12.4 5.774 <sup>×</sup> <sup>10</sup><sup>5</sup> 0.00394 <sup>−</sup>3.8584 <sup>×</sup> <sup>10</sup>−<sup>4</sup> 2.39559 <sup>×</sup> <sup>10</sup>−<sup>7</sup> 0.00432608 1.112276 <sup>11</sup> 7.304 <sup>×</sup> <sup>10</sup><sup>5</sup> 0.00514 <sup>−</sup>4.6934 <sup>×</sup> <sup>10</sup>−<sup>4</sup> 3.30384 <sup>×</sup> <sup>10</sup>−<sup>7</sup> 0.005609674 1.229810 8.5 1.168 <sup>×</sup> <sup>10</sup><sup>6</sup> 0.00854 <sup>−</sup>6.8753 <sup>×</sup> <sup>10</sup>−<sup>4</sup> 6.29435 <sup>×</sup> <sup>10</sup>−<sup>7</sup> 0.009228155 1.440951 6.3 1.845 <sup>×</sup> <sup>10</sup><sup>6</sup> 0.01354 <sup>−</sup>9.7975 <sup>×</sup> <sup>10</sup>−<sup>4</sup> 1.20748 <sup>×</sup> <sup>10</sup>−<sup>6</sup> 0.014520956 1.628762 4.8 2.730 <sup>×</sup> <sup>10</sup><sup>6</sup> 0.01984 <sup>−</sup>1.3018 <sup>×</sup> <sup>10</sup>−<sup>3</sup> 2.08990 <sup>×</sup> <sup>10</sup>−<sup>6</sup> 0.021143915 1.759706 <sup>4</sup> 3.529 <sup>×</sup> <sup>10</sup><sup>6</sup> 0.02524 <sup>−</sup>1.5469 <sup>×</sup> <sup>10</sup>−<sup>3</sup> 2.95677 <sup>×</sup> <sup>10</sup>−<sup>6</sup> 0.026789896 1.831583

**Table 7.** MBE equation calculation procedures for Well-B production data.

*<sup>x</sup>***<sup>1</sup>** = volume changes due to free gas expansion in Organic matrix, *<sup>x</sup>***<sup>2</sup>** = volume changes due to shale matrix shrinkage and stress sensitivity, *<sup>x</sup>***<sup>3</sup>** = volume changes due to desorbed gas volume. Note: *<sup>X</sup>* = is the summation of (x1, x2, and x3).

**Figure 5.** The linear relationship between the production and expansion parameters proposed in the modified material balance equation presented in this study (using the production data from the real-field production data of Well-B).

#### **5. Conclusions**

The accurate estimation of recoverable reserves and production forecasting is incredibly important in shale gas reservoirs. For unconventional reservoirs like shale gas, the material balance equation is a critical tool used to estimate these reserves and forecast production. With the global demand for energy increasing and the potential for shale gas resources to meet this demand, it becomes even more crucial to have precise estimates of recoverable reserves and accurate production forecasts. To address this need, our study

presents a comprehensive derivation of a modified material balance equation for estimating free and adsorbed gas in place within shale gas reservoirs. This modification will enable more accurate estimation of shale gas reserves by considering stress sensitivity and variations in porosity during shale gas production, along with accounting for the effect of critical desorption pressure, water production, and water influx. By establishing a linear relationship between reservoir expansion and production terms, the proposed MBE eliminates the need for complex and iterative calculations usually encountered in existing shale gas MBEs. As a result, this approach offers a simpler yet effective means of estimating shale gas reserves without compromising accuracy. The proposed MBE was validated using an in-house finite element poro-elastic model which accounts for stress redistribution and deformation effects during shale gas production. Moreover, the proposed MBE was tested using real-field data of a shale gas reservoir obtained from the literature. The results of this study demonstrate the reliability and usefulness of the modified MBE as a tool for accurately assessing free and adsorbed shale gas volumes under a coupled poro-elastic environment.

**Author Contributions:** Conceptualization, R.A.A. and S.A.; methodology, R.A.A., S.A. and A.A.; software, R.A.A.; validation, R.A.A., S.A. and A.A.; formal analysis, S.A.; investigation, R.A.A. and S.A.; resources, A.A.; data curation, S.A.; writing—original draft preparation, R.A.A.; writing—review and editing, S.A.; visualization, A.A.; supervision, R.A.A.; project administration, A.A.; funding acquisition, A.A. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research work was funded by Institutional Fund Projects under grant no. (IFPIP: 423-145-1443).

**Data Availability Statement:** The datasets used and/or analyzed during the current study are available from the corresponding author upon reasonable request.

**Acknowledgments:** This research work was funded by Institutional Fund Projects under grant no. (IFPIP: 423-145-1443). The authors gratefully acknowledge technical and financial support provided by the Ministry of Education and King Abdulaziz University, DSR, Jeddah, Saudi Arabia.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **Appendix A**

Figure A1 shows the changes of *X*-component of the effective stresses around the wellbore. It can be seen from Figure A1 that the values of the effective stress simulated numerically are in a good agreement with the analytical solution at different simulation time steps. Hence, the poro-elastic effect led to changes in pore pressure due to fluid production and, consequently, altered the state of the stresses.

Figure A2 shows the pressure distribution around the wellbore due to fluid production after one hour of starting the production process at Pwf = 1000 psi, while Figure A3 shows the contour map of the shear component of effective stress after 1 h of the production process. It can be seen from Figure A3 that the maximum shear stress is around the wellbore and increases away from the wellbore.

**Figure A1.** X component of effective stress analytical versus numerical as a function of radius and time along the x-axis, for σ<sup>H</sup> = 5800 psi, σ<sup>h</sup> = 5500 psi, Pi = 5500 psi, Pw = 1000 psi, kx = 0.01 md, ky = 0.01 md (solid lines are the results obtained from the in-house numerical simulator).

**Figure A2.** Pore pressure contour map after 1 h (numerical results, σ<sup>H</sup> = 5800 psi, σ<sup>h</sup> = 5500 psi, Pi = 5500 psi, Pw = 1000 psi, kx = 0.01 md, ky = 0.01 md).

**Figure A3.** Contour map of shear component of effective stress after 1 h (Numerical results, σ<sup>H</sup> = 5800 psi, σ<sup>h</sup> = 5500 psi, Pi = 5500 psi, Pw = 1000 psi, kx = 0.01 md, ky = 0.01 md).

#### **References**


**Disclaimer/Publisher's Note:** The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.
