**6. Conclusions**

The mechanisms by which circadian rhythm and sleep regulate systemic metabolism and nutritional stimuli from the intestine modulate circadian rhythm and sleep were summarized and discussed dietary therapies could be a novel treatment strategy for both metabolic and sleep disorders, although future studies are needed to validate these strategies.


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

**Reda Abdel Azim 1,\*, Abdulrahman Aljehani 2 and Saad Alatefi 3**

**\*** Correspondence<EMAIL_ADDRESS>
**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−1 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)*, ξl* 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βg* + *Wpβw* − *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*1 is the organic matrix volume change due to free gas expansion (the first term in Equation (26)), *x*2 is the organic matrix volume change due to deformation from matrix shrinkage and stress sensitivity (the second term in Equation (26)), and *x*3 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 *ρg 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 →→ *L* 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, *n* is the number of nodes, → *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*1 and *W*2 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 (R2 = 0.99), a slope of (1.22 × 106), and an intercept of (0.16 × 104). Thus, according to Equation (28) the total free gas in place is equal to 1.2216 × 106 m3, which includes a free gas volume in shale (organic matrix) of 1.22 × 106 m3 and a volume of 0.16 × 104 m3 of free gas in the inorganic matrix. Furthermore, the total free gas in place from the in-house simulator was equal to 1.16 × 106 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.

*x***1** = volume changes due to free gas expansion in Organic matrix, *x***2** = volume changes due to shale matrix shrinkage and stress sensitivity, *x***3** = volume changes due to desorbed gas volume. Note: *X* = 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 × 108 m3 while the intercept (i.e., the free gas volume in the inorganic matter) is equal to 2.25 × 105 m3; moreover, the absorbed gas volume is equal to 4.56 × 107 m3 using equation 29. Accordingly, the sum of these three volumes will give a total gas in place volume of 2.03 × 108 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*2 and *x*3) in Table 6) then the estimated total gas volume will be equal to 2.13 × 108 m3, 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.

*x***1** = volume changes due to free gas expansion in Organic matrix, *x***2** = volume changes due to shale matrix shrinkage and stress sensitivity, *x***3** = volume changes due to desorbed gas volume. Note: *X* = 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 × 108 as reported by Hu et al. 2018 [37]. Figure 5 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 × 108 m3, while the intercept (i.e., the free gas volume in the inorganic matter) is equal to 2.38 × 104 m3; moreover, the absorbed gas volume is equal to 3.72 × 106 m3 using equation 29. Therefore, the MBE calculated total gas in place volume of Well-B is equal to 1.32 × 108 m3 which is in close agreement with the total gas in place of 1.24 × 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 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.

 MBE equation calculation procedures for Well-B production data.

*x***1** = volume changes due to free gas expansion in Organic matrix, *x***2** = volume changes due to shale matrix shrinkage and stress sensitivity, *x***3** = volume changes due to desorbed gas volume. Note: *X* = 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:<PHONE_NUMBER>).

**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:<PHONE_NUMBER>). 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 σH = 5800 psi, σh = 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, σH = 5800 psi, σh = 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, σH = 5800 psi, σh = 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.


*Supporting open minds since 2005*

Biodegradation Technology of Organic and Inorganic Pollutants http://dx.doi.org/10.5772/intechopen.94650 Edited by Kassio Ferreira Mendes, Rodrigo Nogueira de Sousa and Kamila Cabral Mielke

#### Contributors

Shekh Md. Mamun Kabir, Joonseok Koh, Nibal Khaleel Mousa, Abdul-Jabbar Ali, Maha Hussien, Carlos M. Morales-Bautista, Selene Lagunas-Rivera, Sugey López-Martínez, Sebastian A. Ramos-Arcos, Edith G. González-Mondragón, Eduardo S. López-Hernández, Ana R. Rodríguez-Luna, Pedro Eulogio Cisterna Osorio, Barbara Faundez-Miño, Sandra Aparecida de Assis, Raquel Guimarães Benevides, Alison Borges Vitor, Geise Camila Ribeiro, Cleide Souza de Santana, Larissa Emanuelle Da Silva Almeida, Adelaide Perdigão, José Luís da Silva Pereira, Zwahruddeen Muhammad Salisu, Diya'uddeen Bashir Hasan, Yahaya Gambo Liman, Suleiman Mohammed Awal, Ukanah Suleiman Pendo, Mohammed Kabir Yakubu, Ishiaku Semo Umaru, Željko Savković, Miloš Stupar, Nikola Unković, Aleksandar Knežević, Jelena Vukojević, Milica Ljaljević Grbić, Ju-Chen Chia, Mahima Begum, Bhaswati Sarmah, Gayatri Goswami Kandali, Sontara Kalita, Ipsita Ojah, Raktim Bhagawati, Lipika Talukdar, Danushika C. Manatunga, Rohan S. S Dassanayake, Renuka N. Liyanage, Anita Verma, Madhuri Girdhar, Anand Mohan, Kopal Singh, Zeba Tabassum, Vaishali Arora, Babita Khosla, Aaltharackal Vikraman Meera, Manorama Thampatti KC, Jacob John, Bhadra Sudha, Abdulmajeed Sajeena, Selvanayaki Krishnasamy, Ramkumar Lakshmanan, Mythili Ravichandran, Muzafar Riyaz, Rauf Ahmad Shah, Kuppusamy Sivasankaran, Chandan Singh, Deepak Vyas, Eustace Fernando, Godfrey Kyazze, Pavithra Fernando, Ahmed Ahsan, Abubakar Lawal, Lukman Bola Bola Abdul'rauf, Melinda-Haydee Kovacs, Emoke Dalma Kovacs, Sushil Kumar Singh, Manish Singh Singh Rajput, Ravish Singh Rajput, Virendra Kumar, Ashutosh Singh Gautam, Sakshi Singh, Asmae Bouziani, Mohamed Yahya, Dibakar Roy, Dasari Sreekanth, Deepak Pawar, Himanshu Mahawar Himanshu Mahawar, Kamal K. Barman, Kassio Ferreira Mendes, Kamila Cabral Mielke, Rodrigo Nogueira de Sousa, Bruna Aparecida de Paula Medeiros

#### © The Editor(s) and the Author(s) 2022

The rights of the editor(s) and the author(s) have been asserted in accordance with the Copyright, Designs and Patents Act 1988. All rights to the book as a whole are reserved by INTECHOPEN LIMITED. The book as a whole (compilation) cannot be reproduced, distributed or used for commercial or non-commercial purposes without INTECHOPEN LIMITED's written permission. Enquiries concerning the use of the book should be directed to INTECHOPEN LIMITED rights and permissions department (permissions@intechopen.com).

Violations are liable to prosecution under the governing Copyright Law.

Individual chapters of this publication are distributed under the terms of the Creative Commons Attribution 3.0 Unported License which permits commercial use, distribution and reproduction of the individual chapters, provided the original author(s) and source publication are appropriately acknowledged. If so indicated, certain images may not be included under the Creative Commons license. In such cases users will need to obtain permission from the license holder to reproduce the material. More details and guidelines concerning content reuse and adaptation can be found at http://www.intechopen.com/copyright-policy.html.

#### Notice

Statements and opinions expressed in the chapters are these of the individual contributors and not necessarily those of the editors or publisher. No responsibility is accepted for the accuracy of information contained in the published chapters. The publisher assumes no responsibility for any damage or injury to persons or property arising out of the use of any materials, instructions, methods or ideas contained in the book.

First published in London, United Kingdom, 2022 by IntechOpen IntechOpen is the global imprint of INTECHOPEN LIMITED, registered in England and Wales, registration number: 11086078, 5 Princes Gate Court, London, SW7 2QJ, United Kingdom Printed in Croatia

British Library Cataloguing-in-Publication Data A catalogue record for this book is available from the British Library

Additional hard and PDF copies can be obtained from<EMAIL_ADDRESS>
Biodegradation Technology of Organic and Inorganic Pollutants Edited by Kassio Ferreira Mendes, Rodrigo Nogueira de Sousa and Kamila Cabral Mielke p. cm. Print ISBN 978-1-83968-895-9 Online ISBN 978-1-83968-896-6 eBook (PDF) ISBN 978-1-83968-897-3


1. Introduction

Various mathematical theories and simulation methods were developed in the past for describing gas flows in nonequilibrium, in particular, hypersonic rarefied regime. They range from the mesoscale models like the Boltzmann equation [1–6], the direct simulation Monte Carlo methods [7], and the high order hydrodynamic equations [1–6, 8–20]. Among these models, the kinetic Boltzmann equation plays a central role in the hierarchy of PDE-based mathematical models for gas kinetic theory. The kinetic Boltzmann equation can be transformed into the

© 2018 The Author(s). Licensee InTech. This chapter is distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

© The Author(s). Licensee InTech. This chapter is distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0), which permits unrestricted use, distribution, and eproduction in any medium, provided the original work is properly cited.

moment equations by introducing the statistical average in velocity space. Based on the Maxwell's equation of change and the so-called method of moments, Grad [8] in 1949 derived the constitutive equations of viscous shear stresses and heat fluxes from the kinetic Boltzmann equation of the distribution of monatomic gas particles. However, it was found by Grad [9] himself that, within the framework of his constitutive equations, there is a critical Mach number (1.65) beyond which no continuous shock wave solution in high compressive regime is possible.

2. The second-order constitutive model of the Boltzmann equation

The Boltzmann equation plays a central role in the hierarchy of mathematical models for gas kinetic theory. It was derived as an evolution equation for the singlet distribution function of a gas by considering the collision dynamics of two particles and combining it with a statistical molecular chaos assumption. Since the molecular chaos assumption is not of a mechanical nature, that is, the Boltzmann equation is based on the assumptions made to "arrive at it" from the reversible Liouville equations of motion, the Boltzmann equation should be regarded as a fundamental kinetic equation at the mesoscopic level of description of macroscopic processes. Thus, it is a postulate for dynamic evolution of singlet distribution functions f(t,r,v) in the

Numerical Simulation of Hypersonic Rarefied Flows Using the Second-Order Constitutive Model…

� �fð Þ¼ v;r; t C f ; f 2

which cannot be derived from the pure mechanical deterministic consideration. Although it is a first-order partial differential equation in space and time, its solution becomes very complicated because it is nonlinear owing to the collision integral C[ f,f2], which is made up of

The moment equations can be obtained by differentiating the statistical definition of the variable in question with time and later combining with the Boltzmann equation [2–5, 8]; it

The symbols c,u,〈〉,Λ(n) denote the peculiar velocity, the average bulk velocity, the integral in

The conservation laws of mass, momentum, and total energy can be derived directly from the kinetic Boltzmann equation. For example, in the case of momentum conservation law, differentiating the statistical definition of the momentum with time and combining with the

¼ �h i mð Þ v � ∇f v þ mvC f ; f 2

�h i mð Þ v � ∇f v ¼ �∇ � h i mvvf ¼ �∇ � ρuu þ h i mccf � �: (4)

d dthð Þ n � � � �, (1)

http://dx.doi.org/10.5772/intechopen.70657

5

D E � � � �: (2)

� � � � : (3)

� f c�∇hð Þ n D E¼Λð Þ n � hð Þ n C f ;f 2

2.1. The kinetic Boltzmann equation and the method of moments

∂ ∂t

yields for molecular expressions of general moment h(n)

D Eþ chð Þ n f � � D E � f

2.2. Exact derivation of the conservation laws

∂ ∂t

Then the first term on the right-hand side becomes

h i mvf ¼ mv

velocity space, and the dissipation (or production) terms, respectively.

∂f ∂t � �

þ v � ∇

phase space (time, position, velocity),

products of distribution functions.

D Eþ∇� u hð Þ n f

Boltzmann equation yield

∂ ∂t hð Þ n f

After Grad's pioneering work in developing gas kinetic theory and subsequent failure of his 13-moment method in describing hypersonic shock structure, there have been enormous efforts to resolve the problem from various perspectives, not only by physicists and mathematicians, but engineers and also chemists. Among such efforts, Eu's works [2–5] to develop the gas kinetic theory consistent with the second law of thermodynamics beyond the linear irreversible thermodynamics stand out. By recognizing the logarithmic form of the nonequilibrium entropy production, Eu [2] in 1980 proposed a canonical distribution function in the exponential form, instead of Grad's polynomial form. He also generalized the equilibrium Gibbs ensemble theory—providing the relationship between thermodynamic variables and the partition functions—to nonequilibrium processes. It turns out that such canonical exponential form assures the nonnegativity of the distribution function and satisfies the second law of thermodynamics in rigorous way, regardless of the level of approximations.

Recently, Myong [15] in 2014 developed a new closure theory which plays a critical role in the development of gas kinetic theory. The new closure was derived from a keen observation of the fact that, when closing open terms in the moment equations derived from the Boltzmann kinetic equation, the number of places to be closed is two (movement and interaction), rather than one (movement only) misled by the Maxwellian molecule assumption in previous theory. Therefore, the order of approximations in handling the two terms—kinematic (movement) and dissipation (interaction) terms—must be the same; for instance, second-order for both terms, leading to the name of the new closure as the balanced closure. Then, after applying the Eu's cumulant expansion based on the canonical distribution function to the explicit calculation of the dissipation term and the aforementioned new closure, Myong [15] derived the secondorder constitutive models from the Boltzmann kinetic equation and proved that the new models indeed remove the high Mach number shock structure singularity completely, which had remained unsolved for decades.

On the basis of these new theories, this chapter will first describe a recent development in theoretical models for numerical simulation of hypersonic rarefied flows from the viewpoint of the method of moments. It will focus on the detailed derivation of the second-order constitutive model from the original Boltzmann equation and the development of associated computational models for numerical simulation of hypersonic rarefied gas flows in simple geometry as well as complicated real vehicles. Finally, some practical applications of the second-order constitutive model to hypersonic rarefied flows are summarized.


**5. The new trends of early detection of glaucoma and disease progression**

#### **5.1 Detecting of RNFL thinning and disease progression**

Evidence shows that RNFL abnormalities can often be evident without detectable VF damage [32–35]. Therefore, measuring the change of RNFL is likely to be useful in detecting early disease progression. Spectral-domain OCT is now the invaluable investigation tool for glaucoma patients because it can measure RNFL thickness reliably [36] with high sensitivity and specificity to detect glaucoma [37, 38] and its progression [39]. Commercially available software event-based algorithm, such as the Guided Progression Analysis (GPA, Carl Zeiss Meditec), can detect progressive RNFL thinning using RNFL thickness maps. A study has demonstrated that GPA can detect and visualise different patterns of progressive RNFL thinning [39]. Trend-Based Progression Analysis (TPA) is another algorithm for detecting progressive RNFL thinning by measuring the rate of change in RNFL thickness for each superpixel of the RNFL thickness map (50 x 50 superpixels). A study that involved 139 POAG patients (240 eyes) followed up for ≥5 years showed that progressive RNFL thinning determined by GPA and TPA was predictive of detectable functional decline in glaucoma. The study showed that TPA outperformed GPA in detecting more eyes with progressive RNFL thinning at a similar level of specificity (84.2% vs. 81.7% for TPA and GPA, respectively) [40]. Furthermore, TPA also provides visualisation of the distribution of the rate of RNFL thinning. It was suggested that the detection of progressive RNFL loss can serve as a biomarker to reflect disease deterioration behaviour and hence guide glaucoma management [40]. However, TPA is not without its limitations. A minimum of four follow-up visits is required for the construction of the TPA, and performance can be undermined with fewer visits. In situations where there are abrupt RNFL changes or in eyes with large test-retest variability, the event-based analysis may be more useful [41]. The authors concluded that TPA enhances but may not replace GPA for topographic analysis of RNFL thinning.

### **5.2 The dynamic target IOP, disease progression and quality of life: The LiGHT trial**

Clinical trials usually define a treatment IOP-lowering target. For instance, the OHT study aimed for an IOP lowering by 20% from baseline for patients in the treatment arm [3], whereas the Collaborative Normal Tension Glaucoma (CNTG) study targeted an IOP lowering by 30% from baseline [42]. The Laser in Glaucoma and Ocular Hypertension (LiGHT) trial is a multicentre RCT that compared eye drops versus selective laser trabeculoplasty as first-line treatment for POAG or OHT [43]. The study is unique with its well-constructed algorithm for detecting disease progression and guiding treatment escalation [44]. It has a novel approach to defining target IOP. Firstly, the target IOP is specific for each patient at baseline, based on disease severity and lifetime risk of loss of vision at recruitment (e.g. different target pressure and percentage IOP reduction according to the disease stratification suggested by Mills et al. [45]). Secondly, the IOP was adjustable based on IOP control and disease progression [44]. The disease progression (either glaucoma deterioration or conversion of OHT to POAG) was determined by a decision support software based on objective visual field and optic disc imaging criteria. Disease progression was defined as 'strong evidence', with the Humphrey GPA software showing 'likely progression' and/or HRT rim area > 1% per year (at *P* < 0.001), and 'less strong evidence' with GPA showing 'possible progression' and/or HRT rim area > 1% per year (at *P* < 0.01). 'Likely visual field progression' is the presence of three or more points on the GPA at <0.05 probability for change on three consecutive occasions, while 'possible visual field progression' is the same criterion but on only two consecutive occasions [44]. Optic disc progression was defined as the rate of neuroretinal rim loss exceeding 1% of baseline rim area/year on a minimum of five repeat HRT images, where this is equivalent to approximately twice the value of normal age-related rim area loss [22]. Following treatment escalation, there is a resetting of both the target IOP and visual field and optic disc baselines against which future assessments will be compared with.

Although the LiGHT trial is probably more complex in its target IOP setting algorithm compared with other glaucoma trials which defined treatment success based on the proportion of patients achieving a particular target percentage reduction of IOP, it resembled closer to our clinical practice. For instance, further IOP lowering beyond the 'target IOP' is probably required for patients with progressive disease. In some cases, patients might request to reduce medication use even when the target IOP is not achieved (e.g. OHT patients with IOP at 24 mmHg who do not want treatment and show no signs of POAG conversion). Furthermore, the LiGHT trial is also unique in that it included the evaluation of quality of life as an outcome measure. These are all novel features that might become important components for future glaucoma study design.


**Yuriy Shablya 1,***∗* **, Dmitry Kruchinin 1 and Vladimir Kruchinin 2**

Received: 6 May 2020; Accepted: 9 June 2020; Published: 12 June 2020

**Abstract:** In this paper, we study the problem of developing new combinatorial generation algorithms. The main purpose of our research is to derive and improve general methods for developing combinatorial generation algorithms. We present basic general methods for solving this task and consider one of these methods, which is based on AND/OR trees. This method is extended by using the mathematical apparatus of the theory of generating functions since it is one of the basic approaches in combinatorics (we propose to use the method of compositae for obtaining explicit expression of the coefficients of generating functions). As a result, we also apply this method and develop new ranking and unranking algorithms for the following combinatorial sets: permutations, permutations with ascents, combinations, Dyck paths with return steps, labeled Dyck paths with ascents on return steps. For each of them, we construct an AND/OR tree structure, find a bijection between the elements of the combinatorial set and the set of variants of the AND/OR tree, and develop algorithms for ranking and unranking the variants of the AND/OR tree.

**Keywords:** combinatorial generation; method; algorithm; AND/OR tree; Euler–Catalan's triangle; labeled Dyck path; ranking algorithm; unranking algorithm

**MSC:** 68R05; 05C05; 05A15


*Article* **An Improved Pedestrian Navigation Method Based on the Combination of Indoor Map Assistance and Adaptive Particle Filter**

**Zhengchun Wang 1,2, Li Xing 3,\*, Zhi Xiong 1,2, Yiming Ding 1,2, Yinshou Sun 1,2 and Chenfa Shi 1,2**

**Abstract:** At present, the traditional indoor pedestrian navigation methods mainly include pedestrian dead reckoning (PDR) and zero velocity update (ZUPT), but these methods have the problem of error divergence during long time navigation. To solve this problem, under the condition of not relying on the active sensing information, combined with the characteristics of particles "not going through the wall" in the indoor map building structure, an improved adaptive particle filter (PF) based on the particle "not going through the wall" method is proposed for pedestrian navigation in this paper. This method can restrain the error divergence of the navigation system for a long time. Compared to the traditional pedestrian navigation method, based on the combination of indoor map assistance (MA) and particle filter, a global search method based on indoor MA is used to solve the indoor positioning problem under the condition of the unknown initial position and heading. In order to solve the problem of low operation efficiency caused by the large number of particles in PF, a calculation method of adaptively adjusting the number of particles in the process of particle resampling is proposed. The results of the simulation data and actual test data show that the proposed indoor integrated positioning method can effectively suppress the error divergence problem of the navigation system. Under the condition that the total distance is more than 415.44 m in the indoor environment of about 2600 m2, the average error and the maximum error of the position are less than two meters relative to the reference point.

**Keywords:** map assistance; particle filter; global search algorithm; pedestrian navigation


3.4.6. Thermal Stability Studies of HSA and the HSA-QTP System

The thermal stability of HSA without and with QTP was investigated using fluorescence measurements. The fluorescence spectra of the HSA (5 µM) and HSA-QTP (50 µM) complex were recorded (300–400 nm upon excitation at 295 nm) in the temperature range 25–80 ◦C (with 5 ◦C intervals). The solution mixture (HSA-QTP) was incubated for 1 h at 25 ◦C before fluorescence measurements.


**5. Final remark**

Red/dark‐colored fruits are considered healthy and nutritious, the major potential health benefits being a reduced risk for cardiovascular and neurodegenerative diseases. Phytochemicals from red/dark‐colored fruits are also shown to prevent body weight gain, lower blood cholesterol, and reduce cancer risk. Nevertheless, further rigorous, prospective studies are needed in order to bet‐ ter understand the benefits included in red/dark‐colored fruits in our diet. There is also an emer‐ gent interest in the study of red/dark‐colored fruits astringency because of the healthy properties of astringent substances found in red/dark‐colored fruits including antibacterial, antiviral, anti‐ inflammatory, antioxidant, anticarcinogenic, antiallergenic, hepatoprotective, and vasodilating. The role of phenolic compounds and their metabolites as prebiotics, contributing to beneficial gastrointestinal health effects by modulating gut microbial balance with the simultaneous inhibi‐ tion of pathogens and stimulation of beneficial bacteria, should also be highlighted.


**5. Conclusions**

The therapeutic approach to CML has changed since the introduction of the TKI imatinib, followed by dasatinib, nilotinib, which was approved for use in children, and ponatinib, which became the first-line treatment in adults and expanded the therapeutic options, pushing allogeneic stem cell transplantation to a third-line treatment for most pediatric cases. Unfortunately, the selection of a TKI continues to rely on clinical experience in adults [108] without sufficient data on efficacy and safety specific to pediatric patients.

The availability of three TKIs is challenging to choose a first-line option, but it also gives clinicians additional treatment chances in case of a suboptimal response [108].

More experience has been gained about efficacy, toxicity profiles, and comorbidities of Imatinib compared to the other TKIs [22,37,109]. The appearance of resistance or intolerance phenomena has required the development of a second generation of TKI. Drug availability, ease of administration, and financial issues should also be considered.

Evidence-based recommendations have been established for treating CML in adults treated with TKIs. Appropriate guidelines are, however, difficult to extend to pediatric patients given the very rare occurrence in children and adolescents [105]. As a matter of fact, the CML incidence increases with age, from 0.09/100,000 at ≤ 15 years of age to 7.88/100,000 at ≥ 75 years of age. Needless to say, here are several biological and clinical differences between pediatric and adult CML. Markedly increased leukocyte count and a higher incidence of splenomegaly are characteristic features at diagnosis of pediatric patients.

TKIs are designed to inhibit BCR-ABL1 kinase, but they have unfavorable effects, so-called "off-target" complications, such as growth impairment, especially important in children, because they are constantly and actively growing. Long-term morbidity due to TKIs is unknown. Furthermore, the adverse effects on growing children have not been clearly elucidated, even though the exposure period to Imatinib is relatively short. To establish the standard therapeutic management for pediatric

CML, it is important to prospectively confirm the attractive outcomes obtained in adult studies via pediatric clinical trials with a careful monitoring system for TKI-induced adverse effects, especially in growing children [110].

Limited experience with very young children, the transition of teenagers to adult medicine, and the goal of achieving treatment-free remission for this rare leukemia are more significant obstacles that require further clinical investigations [48,108].

Lastly, in order to carry out a possible and viable therapy with TKIs in the pediatric age, a key role is played by developing appropriate formulations specifically customized toward this kind of patients, since these drugs are often available in solid dosage forms, which are difficult to administer in children.

**Funding:** Authors acknowledge Project POR Puglia FESR-FSE 2014–2020 Azione 1.6. Innonetwork progetto D.I.V.A., cod. JD6EDJ7 and the University of Bari "Aldo Moro" (Italy) for their financial support.

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


**4. Discussion**

The results elucidate the microstructural impact of the tribological behavior between the burnishing ball and the steel surface. According to Kuznetsov et al. [12], friction constitutes a fundamental parameter to obtain significant improvements or unwanted effects on surface integrity. Therefore, the consequences of high friction involve the generation of uniaxial tensile stress in the rear zone behind the indenter [12], as well as the increase in the pile-up [21]. However, based on the tribo-contact effects on the selected microstructures, the definitions of high and low friction are ambiguous. The allowable tribo-interaction ranges within the process will be given by the limited plastic deformation of the microstructure. In fact, the stress state conditions after ball burnishing with a 270 N load on the austenitic steel show the presence of a tensile state in the burnishing path direction (x-axis), which is in agreement with the approach of Kuznetsov et al. [12]. A clearly detrimental austenitic surface (reflected in the increase in the macro-texture parameters; Figure 6) and a residual anisotropic state [21,22] beyond the compressive condition (Figure 7) were evidenced. As a consequence, a new finish distribution (skipping tolerances) with valleys and peaks defined by the ball track (Figure 5c) was displayed. This high plastic deformation capacity of the austenitic crystallographic lattice allowed for low surface integrity. Nevertheless, under the same conditions (270 N), the textured martensitic stainless-steel surface offered a contrasting microstructural response to the ball-burnishing process. The lower COF (Figure 4) as an effect of the martensitic matrix, determined the displacement of the peaks toward the milling valleys, conferring uniformity on the surface (Figure 5d), whereas the compressive surface state improved in both directions according to the initial trend established by the milling finish. A higher COF (0.17 after load increment) led to the compressive layer's relocation to the surface (as stated by Kuznetsov et al. [12]), a heightening anisotropy (Figure 7), and an onset of pile-up (Figure 5e, Figure 6e). As seen in Figure 7, the new pressure exerted on the martensitic surface was far from inducing a tensile state in the parallel direction on the burnishing path (which defines high friction [12]) so the hypothesis of high tribo-interactions on this material was limited to the pile-up initiation. Therefore, defining the burnished component's functionality is pertinent. The process configuration must be prioritized, either the contact interactions with other components (roughness) [3] or the exposure to the corrosion, wear, and fatigue conditions whose resistance improves through the generation or increment of the surface compressive residual state [6–8,10,11]. It should be noted that determining the compressive layer's thickness (sub-surface residual state) as a function of the tribo-interaction degree may modify the high-friction hypothesis of martensitic steels established in this study. Nevertheless, excessive tribo-interaction (overloading) within the process must be prevented in order to allow the treated microstructures to retain some degree of ductility, as cited by Kuznetsov et al. [12].