A Dynamic Permeability Model in Shale Matrix after Hydraulic Fracturing: Considering Mineral and Pore Size Distribution, Dynamic Gas Entrapment and Variation in Poromechanics

Abstract

Traditional research on apparent permeability in shale reservoirs has mainly focussed on effects such as poromechanics and porosity-assisted adsorption layers. However, for a more realistic representation of field conditions, a comprehensive multi-scale and multi-flowing mechanism model, considering the fracturing process, has not been thoroughly explored. To address this research gap, this study introduces an innovative workflow for dynamic permeability assessment. Initially, an accurate description of the pore size distribution (PSD) within three major mineral types in shale is developed using focussed ion beam-scanning electron microscopy (FIB-SEM) and nuclear magnetic resonance (NMR) data. Subsequently, an apparent permeability model is established by combining the PSD data, leading to the derivation of dynamic permeability. Finally, the PSD-related dynamic permeability model is refined by incorporating the effects of imbibition resulting from the fracturing process preceding shale gas production. The developed dynamic permeability model varies with pore and fracture pressures in the shale reservoir. The fracturing process induces water blockage, water-film formation, and water-bridging phenomena in shale, requiring additional pressure inputs to counteract capillary effects in hydrophilic minerals in shale, But also increases the overall permeability from increasing permeability at larger scale pores. Unlike traditional reservoirs, the production process commences when the fracture is depleted to 1–2 MPa exceeds the pore pressure, facilitated by the high concentration of hydrophobic organic matter pores in shale, this phenomenon explains the gas production at the intial production stage. The reduction in adsorption-layer thickness resulting from fracturing impacts permeability on a nano-scale by diminishing surface diffusion and the corresponding slip flow of gas. this phenomenon increases viscous-flow permeability from enlarged flow spacing, but the increased viscous flow does not fully offset the reduction caused by adsorbed-gas diffusion and slip flow. In addition to the phenomena arising from various field conditions, PSD in shale emerges as a crucial factor in determining dynamic permeability. Furthermore, considering the same PSD in shale, under identical pore spacing, the shape factor of slit-like clay minerals significantly influences overall permeability characteristics, much more slit-shaped pores(higher shape factor) reduce the overall permeability. The dynamic permeability-assisted embedded discrete fracture model (EDFM) showed higher accuracy in predicting shale gas production compared to the original model.

1. Introduction

Shale, as an unconventional reservoir, necessitates distinct exploitation techniques. The production of shale gas has played a pivotal role in global energy supply [1]. The sedimentary structure of shale amalgamates various minerals, including organic matter, clay, carbonates, feldspars, quartz, and pyrite [2,3,4]. Pores within the shale matrix exist across multiple scales, from nano- to micro-scale, rendering shale formation difficult to characterise [5,6,7]. Additionally, shale gas exploration using the fracturing technique has had a significant impact on the poro-mechanism of shale gas reservoirs (SGR). To enhance the ultimate recovery of SGR, a comprehensive understanding of reservoir characterisation, flowing mechanism, and dynamic petrophysics play key roles in shale gas production [3,8,9,10,11].

For shale characterisation, researchers have systematically classified the minerals found within shale into three predominant categories based on their petrophysical properties: brittle minerals, clay, and organic matter (Figure 1). Studies suggest that the pores within these three minerals contribute to shale permeability [12,13,14]. These include inter-granular pores formed between clay minerals, round pores caused by erosion and grains overlapping in brittle minerals, and intra-granular pores within organic matter [3,15,16,17,18,19,20]. For quantifying the pore structure and distribution within different minerals, researchers have applied different routines to quantitatively analyse the pore size distribution (PSD). The numerical description of PSD addresses the fractural theory and the log–Gaussian distribution. The fractural theory is more commonly used to describe the micro-fracture while the log–Gaussian distribution describes PSD directly from experiment results [10,21,22,23,24,25]. The experimental invasion approach includes mercury intrusion porosimetry (MIP), low-temperature nitrogen adsorption (LTNA), and carbon dioxide adsorption. The non-invasive methodologies include nuclear magnetic resonance (NMR) and optical measurements methods include SEM and focussed ion beam-scanning electron microscopy (FIB-SEM) [5,8,25,26] (Figure 1). However, the single experimental methodology data do not reflect whole scale PSD owing to the range limitations inherent to each method, and there is no whole scale methodology that can both reflect PSD and shape factor concurrently [17,27,28]. Recent research has attempted to describe whole-scale PSD using the integration of nano- and micro-scale approaches [23,29]. However, a simple integration method will engender an accuracy issue at the integration boundary. In this article, to describe PSD, we devised a novel coupling method to integrate the NMR and SEM data based on the mechanism of nuclear relaxivity.

As previously noted, SRG development requires a fracturing technique. Based on observation, a significant amount of hydrophobic organic matter surrounds the brittle minerals. Thus, substantial nano- and micro-pores escalate the capillary phenomenon during fracture, causing substantial fluid leakage during fracturing [13,14,30,31,32,33]. The water imbibition process caused by the fluid leakage compressed the original gas in place and displaced it into a deeper location within the pores under the fracturing pressure. the imbibed fracturing fluid will race to control the adsorption sites of the adsorbed gas owing to the stronger adsorption capacity of water molecules, resulting in a reduction in the adsorption layer thickness and water-blockage. As the pressure of the fracturing process decreases, the entrapped gas expands and drains the water at the pore entrance [4,13,34]. The imbibition shows a difference based on mineral petrophysics, and fracturing fluid leakage results in water blockage in pores, drastically impacting the production of shale gas [7,26,35,36]. Quantitative analysis of fracturing fluid leakage in pores and its corresponding impact on gas production has been explained through experiments conducted on capillary tubes with varying radii. Additionally, scholars have utilized molecular simulations to describe the impact on the reduction of adsorbed gas due to fracturing fluid leakage, a factor proven to greatly influence the gas’s apparent permeability [2,17,18,37].

Water and gas transportation within shale formations covers multiple scales and follows different flow regimes owing to the petrophysics of SGR [12]. Confinement in nano-micro pores can alter the phase behaviour of fluids in the pores, and, consequently, the flow regime owing to the distortion of vapour–liquid equilibrium conditions [19,21,38,39]. Based on the Knudsen number, the flow regime can be classified as adsorption gas diffusion, slip flow, and viscous flow. Recent research [35,40] involving accurate explanation has focussed on free gas adsorption. Further explanations related to the flow regime in SGR will be presented later in the article.

Conventional numerical description for dynamic permeability in SGR integrated stress-dependent pore radii variation and poro-elastic poromechanics. Figure 2 illustrates the experimental results from multiple studies. It is well-documented [17,41] that geological characteristics during SGR depletion and production are a deterministic index of shale gas capability in hydrocarbon generation and production process. Numerous researchers have investigated the impact on permeability from dynamic radius at varying stress levels [42,43,44,45,46]. Based on numerical simulations, researchers have concluded that effective stress on pores will change the pore formation, and, thus, significantly impact the apparent permeability by changing the flow mechanism. Additionally, the matrix undergoes a slight increase early during production and becomes larger at the later stage because of slip flow and Knudsen diffusion [30,38,47,48,49].

This study aimed to comprehensively address all key attributes of dynamic permeability in Shale Gas Reservoirs (SGRs) through a detailed workflow (Figure 3). Firstly, complex PSD characterization is conducted based on the stitched FIB-SEM data and NMR data, including imaging, segmentation, and statistical analysis to extract equivalent pore radius at different scales. The whole-scale PSD of equivalent granular-shaped pores is characterized by coupling both experimental results at this stage. After that, the PSD in shale is separated by three minerals with the help of Bayesian statistics for Log-Gaussian characterization. The characterized PSD includes probability density at various radii and the shape factor of clay pores. This Log-Gaussian-based PSD description in different minerals is then utilized for radius correlation considering water film thickness under field conditions and to extract the critical liquid bridging parameters. After this, the correlated PSD becomes a crucial factor influencing the fracturing-induced imbibition process. During this stage, the force balance at different imbibed conditions under varying pressures and their corresponding impact on the adsorbed gas in pores are characterized. Finally, with the probability density at each correlated pore scale, critical water blockage, and liquid bridging condition at different pressures under varying radii, the apparent permeability was calculated at each correlated equivalent radius in the capillary bundle model. Combined with probability density, liquid bridging, and water blockage, the dynamic permeability encompasses all the aforementioned flowing mechanisms and water–gas conditions to reflect real-world field conditions.

To validate the dynamic permeability model, it was integrated into the EDFM simulator as an amendment permeability module. The proposed model demonstrated an increased precision in reservoir production prediction compared to the original simulated results.

2. Complex Pore Characterisation

2.1. Description of PSD within Different Minerals

For shale characterisation, numerous studies [21,48,49] have examined dynamic petrophysics from pores within different minerals. Drawing insights from visualisation technology, our proposed mineral-pore system is divided into three different pores: granular-shaped inter-pores between brittle minerals; inter-pores between tetrahedral sheets-shaped clay minerals; and intra-pores within organic matter for characterisation (see Figure 4).

Quantifying PSD involves a diverse set of implementations across the nano- and micro-scales. Numerous studies contribute to the realm of accurate PSD description by utilising cutting-edge visualisation techniques and characterising interactive systems within three types of minerals [8,13,50,51]. However, those pore size characterised methodologies with merits and limitations, the visualisation methods such as X-ray CT, FIB-SEM, and TEM are capable of describing the micro-nano pore and mineral shape from a small 2D planer or 3D structure, whereas the adsorption testing and NMR can reflect the PSD based on a shale sample scale. Additionally, different implementations reflect varied results at different scales because of the theory behind the methodology. In this study, for complex shale PSD characterisation, a coupled NMR and FIB-SEM methodology is used to accurately reflect the valid PSD based on the NMR methodology and the characterisation of pore geometry based on FIB-SEM.

In our study, nano-micro PSD is extracted from FIB-SEM. Sup-resolution image from the WY formation is created utilising the FIB-SEM line-by-line scanning image stitching technology [52]. To characterise the PSD in shale, we stitched 16,000 FIB-SEM images with 1024 × 1024 pixels (50 nm per pixel) via python codes, pores within sup-resolution images are extracted from the greyscale threshold segmentation software, the surface pore area and related probability density, shape factor statistics are extracted form the stiched FIB-SEM image. combined with the image segmentation software (Figure 5), the equivalent pore distribution with radius is shown in Figure 6. Based on the box-plot statistics, we conclude that there are two spikes in the boxplot statics above 0.1 um. Consequently, referring to different nano-micro test results [3,53] and different PSD models [10] (Figure 7), we assume the nano-micro scaled PSD of the WY formation follows the log–Gaussian distribution.

2.2. Coupling PSD by Comparing NMR and FIB-SEM Data

To overcome the limitations of the coupled PSD distribution in previous studies, this study utilised both FIB-SEM and NMR data to comprehensively clarify pore structure and PSD [2,48,54]. FIB-SEM can directly observe the PSD within different minerals in the shale matrix, while NMR transfers the signal decay to obtain the volume of pores. The transformation between pore volume and NMR relaxivity data was described based on the fundamental principles of the porous media NMR theory [55].

$$\frac{1}{T_{2,S}}={\rho }_2{S}_{\nu }$$

(1)

where, $S_{\nu }$ is the ratio of surface area to pore volume; $T_2$ is the surface relaxivity—a value that characterises the strength of relaxation induced by the solid/fluid interface.

As mentioned earlier in this study, pores in shale are considered to have granular and slit geometries. The NMR interpretation alone cannot directly illustrate the geometry of pores, and FIB-SEM has limitations in describing nano-scaled pores in shale due to resolution constraints. To comprehensively characterise the PSD and the corresponding shape factor, we integrated NMR data and FIB-SEM data using information from larger granular-shaped organic matter pores (OMPs) in the range of approximately 100 nm to 1 $\mathsf{\mu }$m, where pore scale can be characterised by both methodologies.

PSD in larger OMPs is directly described from visualized FIB-SEM data, while the original T₂ data interpretation relies on the pore structure. By adjusting the S_v, the PSD curve from interpreted NMR data aligns with FIB-SEM results, and the corresponding S_v ratio reflects the equivalent granular pore capillary bundle PSD. Consequently, an equivalent granular-shaped capillary bundle model is generated from NMR data with FIB-SEM coupled S_v to demonstrate the full-scaled distribution in the shale matrix.

Thus, the PSD at varying equivalent radii in shale is considered as:

$$P\left(r\right)=\sum \frac{1}{rlg\sqrt{2\pi }\sigma }\mathit{EXP}\left[-\frac{{\left(lgr-\mu \right)}^2}{2{\sigma }^2}\right]$$

(2)

where $\sigma$ and $\mu$ are standard derivations and mean; r is the equivalent radii.

2.3. Log–Gaussian Mixed Capillary Model for Describing Overlapped Pores Distribution in Clay and Brittle Minerals Based on Coupled Data

Based on multiple scholarly studies, it is evident that the PSD in various shale formations cannot be simplistically characterised as a mere combination of Gaussian distributions arising from organic matter and inorganic pores [23]. Notably, certain PSD tests conducted on the Barnett shale did not exhibit the conventional bimodal distribution with two distinct peaks in the curve. This underscores the complexity and unique characteristics of shale formations, challenging the application of a one-size-fits-all approach [27] in describing their PSD [25,28]. Fortunately, research into distinct pore types has revealed that the pore distribution at various equivalent radii can still be effectively described based on the probability density of the three predominant minerals. The portrayal of PSD within the shale matrix is described as follows:

$$P\left(r\right)=\sum _{}^K{c}_{k}P(r,{\mu }_k,{\sigma }_k)\left(k=b,c,om\right)$$

(3)

$$P\left(r|{\mu }_k,{\sigma }_k\right)=\frac{1}{\sqrt{2\pi }{\sigma }_k}{e}^{-\frac{{\left(r-{\mu }_k\right)}^2}{2{\sigma }_k}}$$

(4)

where, $r$ is the equivalent pore radius; ${\sigma }_k$ and ${\mu }_k$ are the standard deviation and mean, respectively, of the corresponding pore distribution of brittle mineral $b$, clay mineral $c$, and organic matter $om$. $c_k$ reflects the volume concertation of brittle, organic matter and clay, respectively.

As mentioned, the micro-scale distribution is considered a Gaussian mixed model (GMM) of three clusters. To accurately interpret the slit-like PSD, the Gaussian distribution of equivalent clay pores can be extracted, thanks to the Expectation Maximisation Algorithm. A latent variable, z is added to the algorithm to describe probability of categorisation. Considering the full range probability of categorising a random pore includes two results: clay mineral pores and brittle mineral pores. The latent variable, $z_k$= 1 reflects the correct categorisation. The full range probability can then be described with the help of latent variable $z\in$ {(0, 1) (1, 0)} by describing the random categorisation based on Bayesian interpretation. The joint probability distribution of z, as follows:

$$p\left(z\right)=p\left(z_1\right)p\left(z_2\right)\cdot \cdot \cdot p\left(z_K\right)=\prod _{k=1}^K{\pi }_k^{z_k}$$

(5)

$$P\left(r|z\right)=\prod _{}^{K}P\left(r|{\mu }_k,{\sigma }_k\right)$$

(6)

Because, $z_k\in \left(0,1\right)$. Thus, the conditional probability of categorising pores are described as:

$$P\left(r|z\right)=\sum _{}^K{c}_{k}P\left(r|{\mu }_k,{\sigma }_k\right)\left(k=b,c,om\right)$$

(7)

$c_k$ is the probability that certain mineral pores will be selected, which are described via volume concertation. Based on the algorithm mentioned above, we already know the prior probability $P\left(z\right)$ and likelihood $P\left(r\right|z)$. The posterior probability $P\left(z\right|r)$ is described as follows:

$$P\left(z|r\right)=P\left(z_k=1|r\right)=\frac{P\left(z_k=1\right)P\left(r|z_k=1\right)}{\sum _{j=1}^{K}P\left(r|{\mu }_k,{\sigma }_k\right)}$$

(8)

with prior probability, likelihood, and posterior probability from the Bayes decision theory, the perimeter ${\mu }_k$ and ${\sigma }_k$ of the GMM, which represents the two majority micro-scale pores, can be accurately predicted using the converged expectation–maximization algorithm with the help of a Python-coded iteration function. The schematic workflow is shown in Figure 8. Thus, the equivalent capillary bundle model reflecting the brittle mineral and entrapped clay pores will be characterised.

2.4. Capillary Bundle Model for Describing Shale Matrix

The process of constructing a shale matrix model based on PSD extraction is illustrated in Figure 9. Initially, the FIB-SEM image was analysed using a segmentation algorithm to extract pores ranging from micro- to meso-scale, including slit-like pores between clay minerals. Thereafter, the data were processed to eliminate the impact of micro fractures. By combining the nano-scale NMR data, an equivalent capillary tube could be established, facilitating a comprehensive characterisation of the shale matrix.

2.5. Tortuosity of Capillary Tube Model

However, the GMM-based capillary tube model cannot directly reflect the flow through porous media and the equivalent capillary length cannot be simply explained as the length of shale matrix [45]. The dimensionless parameter, $\tau$, which is called hydraulic tortuosity, is used to describe the average length of all particle path lines that pass through a cross-section—in this case, the shale matrix. Hydraulic tortuosity was first proposed by Kozeny via the well-known Kozeny–Carman equation [56]:

$$K=\frac{{\varphi }^3}{{\beta \tau }_h^2{S}^2}$$

(9)

where, ϕ is the porosity defined as the fraction of pore space in the porous medium; S is the specific surface area equal to the ratio of interstitial surface area to bulk volume; and β is the shape factor (a constant depending on the type of granular material).

However, this method is not a perfect fit for the capillary bundle model in shale owing to its high heterogeneity. Another representative theoretical model was derived from a fixed bed of randomly packed identical particles [39]:

$$\tau =1.23\frac{{\left(1-\varphi \right)}^{4/3}}{{\xi }^2\varphi }$$

(10)

where, ξ is the shape factor (sphericity) of the particle; ξ = 1 for sphere particle and ξ < 1 for non-spherical particles. In this study, the length of the capillary bundle model considering hydraulic tortuosity for each mineral pore is described as follows [57]:

$$\Omega \left(L\right)=1.23\frac{L{\left(1-{\varphi }_i{c}_i\right)}^{4/3}}{{{\xi }_i}^2{\varphi }_i}$$

(11)

2.6. Stress Dependent Pore Radius Correlation

2.6.1. Pores with Brittle Mineral and Clay Mineral

Shale is viewed as a naturally fractured reservoir sensitive to changes in stress conditions and increased effective stress decreases the fracture aperture. Based on laboratory research, the size changes for non-organic matter (brittle minerals and clay) are not as strong as for organic matter [58,59]. In this study, the equation applied for describing brittle minerals and clay minerals is as follows:

$$\varphi ={\varphi }_r+\left({\varphi }_i-{\varphi }_r\right)e^{\left(a{\sigma }_m^{\prime }\right)}$$

(12)

where, $\varphi$ represents the porosity of the medium, ${\varphi }_r$ signifies the residual porosity under high-stress conditions, ${\varphi }_i$ denotes the initial porosity, a is a media property determined through experimental methods, and ${\sigma }_m^{\prime }$ prime stands for the mean effective stress. The numerical values for parameters $a$ have been extracted from an earlier investigation [60].

Thus, assuming the length in the equivalent capillary bundle model in the shale matrix remains constant during production, the radius in the equivalent capillary bundle model is described as follows:

$$\varnothing \propto R^2\Omega \left(L\right)$$

(13)

Considering the tortuosity-compensated capillary length remains constant during reservoir depletion, substituting (13) into (12) the radius changes in the circular-shaped equivalent capillary bundle model could be described as follows:

$${R_{iom}}^2\Omega \left(L\right)={R_r}^2\Omega \left(L\right)+\left({R_i}^2\Omega \left(L\right)-{R_r}^2\Omega \left(L\right)\right)e^{\left(a{\sigma }_m^{\prime }\right)}$$

(14)

$$\Lambda {\left(r\right)}_{iom}=\sqrt{{R_r}^2\left(1-e^{a{\sigma }_m^{\prime }}\right)+r^2{e}^{a{\sigma }_m^{\prime }}}$$

(15)

2.6.2. Water Film Thickness Considering Adsorption

The stability and thickness of adsorbed films on substrates are primarily influenced by force interactions between the thin film and the solid surface, which can be measured by the disjoining pressure. Many studies have investigated the thickness of water films in geological media [1,61,62]. In the thermodynamic system of water film adsorption on a flat substrate, the relationship between disjoining pressure and relative humidity can be described [62]. To consider the structural force caused by OMP, Li [62] proposed a novel model that takes into account both hydrophobic and hydrophilic minerals. They consider the water film with slit-like pores as the combination of surface force on both slit-like pores. Hence, the water films adhering to slit surfaces maintain a flat configuration, and due consideration has been given to the capillary pressure arising from the gas–water meniscus interface.The disjoining pressure for capillary tube and slit-like pores can be described as follows [16,63]:

For equivalent sphere pores:

$${\Pi }_{capillary}\left(h\right)=\frac{r}{r-h}·{\Pi }_{flat}\left(h\right)+P_c=\frac{r}{r-h}·{\Pi }_{flat}\left(h\right)+\frac{\gamma }{r-h}$$

(16)

$${\Pi }_{flat}\left(h\right)={\Pi }_m\left(h\right){+\Pi }_e\left(h\right)+{\Pi }_s\left(h\right)$$

(17)

$${\Pi }_m\left(h\right)=-\frac{A_{gws}\left(15.96\frac{h}{l}+2\right)}{12\pi h^3{\left(1+5.32\frac{h}{l}\right)}^2}$$

(18)

$${\Pi }_e\left(h\right)=2{\epsilon }_r{\epsilon }_0{\left(k_{B}T/ze\right)}^2·K^2$$

(19)

$${\Pi }_s=k{e}^{-\frac{h}{\lambda }}$$

(20)

where, ${\Pi }_1\left(h\right)$ denotes the force between the water film and the surface on the same side, accounting for both short- and long-range interactions, MPa; Similarly, ${\Pi }_2\left(h\right)$ represents the force between the water film and the opposing surface; ${\Pi }_3\left(h\right)$ signifies the force between two adsorbed water films, considering only molecular action, Mpa. It’s worth noting that Awgw in $\Pi$₃(h) represents the Hamaker constant for liquid-gas-liquid interactions caused by two adsorbed water films.

The relationship between water film thickness $h,r$ described by disjoining pressure $\Pi \left(h\right)$ and $P_v/P_0$ in the case of water film inside granular and slit-liked capillary can be described via the implicit function, the corresponding relationship between:

$${\Pi }_{capillary}\left(h\right)=-\frac{RT}{V_m}\mathrm{l}\mathrm{n}\frac{P_v}{P_0}$$

(21)

The water film thickness in the capillary bundle model at reservoir condition (23 Mpa) is shown in Figure 10. The designated gas flow radius based on the capillary bundle radius is described as follows:

$$\Phi \left(r\right)=r-h\left(T,V_m,P_v/P_0\right)$$

(22)

2.6.3. Pores in Organic Matter

Organic matter in shale reservoirs is considered as both source rock and porous media for gas. A high total organic carbon (TOC) implies greater adsorption and desorption capacity. As shale reservoir depletes, adsorbed gas begins to desorb from the pore surface. The matrix shrinkage process is relatively similar to that of coalbed methane with coal; recently, studies have referred to the model of coal to quantitatively analyse this phenomenon [64,65]. In this study, a well-known mathematical model [47] is applied to calculate the variation of permeability for the organic matter. Additionally, many studies regarding coalbed methane have revealed the relationship between porosity and permeability, based on which, researchers have explored the relationship between porosity and permeability of organic matter [24]. For shale reservoir under specific stress, the porosity change and porosity–permeability relationship before and after depletion are described as follows:

$$k=k_i{e}^{-3c_f\Delta \sigma }$$

(23)

$$\Delta \sigma =-\frac{v}{1-v}\left(p-p_i\right)+\frac{E}{3\left(1-v\right)}{\epsilon }_l\left(\frac{p}{p+p_{\epsilon }}-\frac{p_i}{p_i+p_{\epsilon }}\right)$$

(24)

$$\frac{k}{k_i}={\left(\frac{\varphi }{{\varphi }_i}\right)}^3$$

(25)

where, $k$ represents media permeability; $k_i$ stands for the initial media permeability; $c_f$ denotes media compressibility, typically determined by fitting the correlation curve to laboratory test data; $\Delta \sigma$ represents the change in effective stress; ν is Poisson’s ratio; $p$ signifies media pressure; p_i is the initial media pressure; $E$ stands for Young’s modulus; and ε_l and pε denote the Langmuir-type organic matter shrinkage constants. Thus, the pore radius described by the pressure sensitivity are shown in:

$$\Lambda {\left(r\right)}_{om}=R_i\sqrt[6]{e^{-3c_f\frac{v}{1-v}\left(p-p_i\right)+\frac{E}{3\left(1-v\right)}{\epsilon }_l\left(\frac{p}{p+p\epsilon }-\frac{p_i}{p_i+p\epsilon }\right)}}$$

(26)

3. Dynamic Permeability Model

3.1. Flowing Mechanism Considering Hydraulic Fracturing Process

It is conceptualised that the micro-fracture system comprises three parts: the granular-shaped organic matter pores; brittle mineral pores; and slit-like entrapped pores within clay. The conceptual model for gas transportation after gas production is illustrated in Figure 11. Smaller pores in hydrophilic minerals (brittle, clay) in the shale with higher capillary forces are comparable easily to drain the water into the pore space, resulting in water blockage, while larger hydrophilic minerals with smaller capillary force are required smaller pressure differences to overcome water imbibition caused from hydraulic fracturing. organic matter pores suffered the least impact from water imbibition processes due to their hydrophobic characteristics. As the reservoir is depleted, the original gas in the pore will displace the imbibed water outside the pore, water film will form around the pore wall. Thus, the gas flow of pores within the shale matrix can be recognized in two categories, fully water-saturated pores due to imbibition, and water-film-affected pores. Brittle minerals and clay minerals pore with different shape factors displaced different shapes in water film.

To quantitatively analyse the pore, during the fracturing process, substantial fracturing fluid was pumped into the deeper shale matrix throughout the capillary tunnel [14]. Owing to differences in petrophysics, the flow mechanism for the capillary tube in different minerals varies [66]. If the water is leaking off into open system throats, gas will be displaced into nearby pores or throats, the closed system takes account of the gas entrapment in a single throat, recent studies have explored the imbibition process as in a closed system (Figure 11), based on its ultra-low permeability [13,26] and experiment results [14]. The concepts of open and closed systems are illustrated in Figure 12. In this article, to express the fracturing process, we assume the capillary bundle model in a closed system. However, the capillary bundle model is not applicable for super heterogeneity shale reservoir because the model considered the capillary bundle as in a single mineralogy, corresponding water-blockage phenomenon and in the proposed model will be drastically influenced at the muti-mineralogy capillary model.

For the interpretation of water saturation and the water-blockage phenomenon in the proposed model, we consider the force at both sides of the capillary tube (Figure 13). Fluid and gas saturation in the conceptualised model during production are analysed via the force balance equation at the water–gas interface in pores. The force balance can be described during the fracking and flow-back periods for the water–gas characterisation. Eventually, the water-blockage phenomenon at different pressures can be recognised. Following the force-balanced equation considering compressed OGIP, capillary pressure, fracturing liquid pressure, and fracturing pressure are investigated. The water blockage phenomenon at different pores and fracturing pressures can be characterised during the fracturing period in the equivalent capillary bundle model [38]. The force balance equation at different minerals is shown as follows.

For OMP:

$$\frac{P_{g0}{L}_0}{L_0-x_0}+\frac{4\sigma \mathit{cos}{\theta }_{omp}}{d}-P_f+\frac{4ZRT{n}_{omp}x\left(p\right)}{d\left(L_0-x\left(p\right)\right)}=0$$

(27)

For brittle minerals pores:

$$\frac{P_{g0}{L}_0}{L_0-x_0}+\frac{4\sigma \mathit{cos}{\theta }_b}{d}-P_f+\frac{4ZRT{n}_{b}x\left(p\right)}{d\left(L_0-x\left(p\right)\right)}=0$$

(28)

For entrapped clay minerals pores:

$$\frac{P_{g0}{L}_0}{L_0-x_0}+\frac{4\sigma \mathit{cos}{\theta }_c}{d}-P_f+\frac{4ZRT{n}_{c}x\left(p\right)}{d\left(L_0-x\left(p\right)\right)}+n{E}_{\pi }RT\left(C_0-C_f\right)=0$$

(29)

where, $x_0$ is the equivalent initial water saturated capillary length; $\sigma$ is the interfacial tension of water N/m; and $L_0$ reflects the tortuosity-considered equivalent capillary length. In this article, they are described as follows:

$$L_0=\Omega \left(L\right)L=1.23\frac{L{\left(1-{\varphi }_i{c}_i\right)}^{4/3}}{{{\xi }_i}^2{\varphi }_i}$$

(30)

3.2. Water Bridging Mechanism during Production

As previously mentioned, a water film is created because of force interactions. During the flow-back production period, as water is displaced from the capillaries through imbibition, hydrocarbon molecules begin to get adsorbed in the organic matter pores (OMP). Nevertheless, a thick water film persists within the capillary bundle, regardless of wettability. Surface forces become dominant in stabilising the liquid film. As the water film thickens, a liquid bridge forms, eventually spontaneously filling the entire capillary. This phenomenon of transition from the vapour phase to the liquid phase, known as ‘capillary condensation’, significantly impacts the apparent permeability.

Researchers [3,66] have discussed the characteristics of water-film thickness and stability within nano-pores and developed an approach considering fluid/pore-wall interactions to describe the phase behaviour of thin water-film transition into liquid condensation (Figure 14). They have identified that the instability mechanisms of adsorbed films differ inside slits and capillaries, with capillaries having higher surface interactions and easier condensation. The phase behaviour of adsorbed water film within nano-porous montmorillonite and shale was investigated, with water condensation found in hydrophilic clay samples, but only partial condensation in shale samples owing to the hydrophobic repulsion of organic minerals. This may impact the storage and transportation of gas in these nano-pores.

Generally, the instability condition of the thin wetting film could be determined by disjoining pressure isotherms, which was formulated as ∂Π(h)/∂h > 0 [67]. Thus, the critical film thickness, h^⁎ can be characterised as follows:

$${\left.\frac{\partial \Pi \left(h\right)}{\partial h}\right|}_{h=h^{\mathrm{*}}}=0$$

(31)

As per the previously derived water-film thickness equation, the corresponding critical water-film thickness within the capillary bundle model can be described as follows:

$$\Pi \left(h^{\mathrm{*}}\right)=-\frac{RT}{V_m}·\ln \frac{P_{\nu }^{\mathrm{*}}}{P_0}$$

(32)

The corresponding vapor pressure, $P_{\nu }^{\mathrm{*}}$ and saturated vapor pressure, $P_0$ are described as follows [68]:

$$P_0=7.0746-\left(1657.46/\left(T-273.15+227.02\right)\right)$$

(33)

$$P_{\nu }^{\mathrm{*}}=c_{water}P$$

(34)

3.3. Mixed Flow Regime

The capillary bundle model is considered stacked-up flow channels, each flow channel can be described as independent, such flows within different pore sizes show different flow regimes, have higher importance in engineering and are easy to get analytical solutions. The equivalent capillary bundle model adheres to the flow regimes within porous media, accommodating the coexistence of various flow regimes. Considering the pore-scale distribution in different minerals within the shale matrix, we categorise the flow regimes within a single capillary tube during shale gas flow-back production. This includes molecular flow from the adsorption layer and the Knudsen layer-affected viscous flow (Figure 15). Researchers have demonstrated that these two dominant flow regimes within nano-micro structures could be investigated independently [69,70,71]. For permeability calculations, both flow regimes could be considered independently. Subsequently, permeability could be calculated by adding free gas molecular flow and slip-correlated viscous flow.

Gas-flow regimes in unconventional reservoirs are divided by the K_n (ratio of a molecular mean free path (MFP) to the average channel size). Considering a dynamic pore aperture, the effective K_n is described as follows:

$$K_n=\frac{\lambda }{b_{eff}}$$

(35)

$$\lambda =\frac{{\mu }_g}{P_p}\sqrt{\frac{\pi ZRT}{2M}}$$

(36)

where,

$${\mu }_g={\mu }_{g0}\left[1+\frac{Y_1}{T_r^5}\left(\frac{P_r^4}{T_r^{20}+P_r^4}\right)+Y_2{\left(\frac{P_r}{T_r}\right)}^2+Y_3\left(\frac{P_r}{T_r}\right)\right]$$

(37)

$$Z=\left(0.702e^{-2.5T_r}\right)P_r^2-\left(5.524e^{-2.5T_r}\right)P_r+\left(0.044T_r^2-0.164T_r+1.15\right)$$

(38)

$$P_r=P_p/P_g\hspace{0.25em},\hspace{0.25em}{T}_r=T/T_g$$

(39)

M is the molecular mass, dimensionless; $b_{eff}$ is the effect pore radius, m; ${\mu }_{g0}$, is the viscosity at standard condition, Pa-s; and $P_r$ and $T_r$ are the relative pressure and temperature, respectively. MPa; $Y_1,Y_2\mathrm{a}\mathrm{n}\mathrm{d}{Y}_3$ are the fitting curve coefficient, dimensionless; corresponding data will be displayed in

Table 1. Basic parameters of the dynamic permeability model for shale gas reservoir.

Parameters	Value
The maximum pore spacing, $R_{fmax}$ (μm)	10 [54]
The minimum pore spacing, $R_{fmin}$ (μm)	0.001
Porosity for fracture networks, $\varphi$	0.05 [57]
TMAC, $f$	0.79 [24]
Temperature, $T$ (K)	393
Pore pressure, $p$ (MPa)	23
Confining pressure, $Pc$ $\left(\mathrm{M}\mathrm{P}\mathrm{a}\right)$	100.56
Critical temperature, $T_c,\mathrm{K}$	190.6
Critical pressure, $P_c,$ Pa	4.599 × 106
Fitting constant, Y1	7.9
Fitting constant, Y2	−9 × 10−6
Fitting constant, Y3	0.28
The Langmuir pressure, $P$ (MPa)	6.72 × 106
Adapted coefficient (Organic matter), a	5 × 10−8
Media coefficient (in-organic matter), $C_f$ (1/Mpa)	22
Molecular weight, $M$ (kg/mol)	16 × 10−3
Universal gas constant,$R$ (J/(mol·K))	8.314
Interfacial tension (water-gas), $\sigma ,\mathrm{N}{\mathrm{m}}^{-1}$	25 [38]
Rarefaction coefficient, $\alpha$	1.19 [24]
Diffusion coefficient, $D_a^0$	1 × 10−12 [24]
Adsorption isothermals,$∆H$, $\mathrm{J}\mathrm{m}\mathrm{o}{\mathrm{l}}^{-1}$	14 × 103
Electrical conductivity, e, ($\mathrm{F}{\mathrm{m}}^{-1}$)	8.85 × 10−12
Potential difference (solid-liquid), $\mathrm{m}\mathrm{V}$	50
CH4 molecular radius, $\mathrm{m}$	0.38 × 10−9
CH4 viscosity, $\mathrm{M}\mathrm{P}\mathrm{a}·\mathrm{s}$	0.0184
Water concentration in $\mathrm{C}{\mathrm{H}}_4$	0.001

and

Table 2. The parameters of different minerals.

Parameters	OMP	Brittle	Clay
Porosity concentration	0.9	0.03	0.07
STD derivation	0.32	0.17	0.23
Mean	0.97 × 10−6	0.03 × 10−6	0.062 × 10−6
Contact angle (°)	116	10	12
Adsorbed gas concentration (mol/m2)	7 × 10−6	1.92 × 10−6	1.8 × 10−6

.

The viscosity correlation factor considering the rarefaction effect in capillary structure is as follows:

$$u_g=\frac{u_{g0}}{1+\alpha Kn}$$

(40)

3.3.1. Free Gas Viscous Flow after Fracturing Process

In smaller Knudsen-numbered equivalent capillaries, viscous flow predominates compared to the previously mentioned flow regime. Unlike the earlier-discussed flow regimes, viscous flow is highly sensitive to the geometry of the flow region. Within the equivalent capillary bundle context mentioned earlier, viscous flow in the mixed-flow regimes is distinguished between rectangular and round shapes. For entrapped clay minerals pores, the velocity distribution within a constant wide-height ratio slit-like clay pore is described using Fourier’s equation:

$$v\left(x,y\right)=\frac{\mathsf{\nabla }P}{u_g}\frac{h^2}{12}\sum _{n=1,3,5\dots }^{\infty }\frac{48}{{\pi }^3{n}^3}\left(1-\frac{cosh\left(n\pi x/h\right)}{cosh\left(n\pi w/2\right)}\right)sin\left(n\pi y/h\right)$$

(41)

For the inter-pore within the shale matrix such as brittle minerals pores and larger-scaled OMP pores, the viscous flow velocity along the radius could be expressed as follows:

$$v\left(r\right)=\frac{{R_f}^2}{4{\mu }_g}\left(1-r^2\right)\mathsf{\nabla }P$$

(42)

3.3.2. Adsorption Gas Diffusion with Impact of Fracturing Process

As the reservoir depletes, adsorbed gas in the organic matter pore surface diffuses along the flow-back direction. This phenomenon is explained by the Max-Stefan theory, given that the gas molecules are adsorbed as a single layer within entrapped organic matter pores and the diffusion process is powered by the electrochemical potential gradient. Assuming the desorption-diffusion process follows the Langmuir iso-thermal adsorption curve and the spherical shape of organic matter pore, the shale gas diffusion velocity in the adsorption layer driven by electrochemical gradient is given as follows [45]:

$$v_d=D_d{C}_a\frac{ZRT}{P^2}\mathsf{\nabla }P$$

(43)

where, $D_d$ is the diffusion coefficient, ${\mathrm{m}}^2/\mathrm{s}$, $C_a$ is the molar concentration of adsorbed gas $\mathrm{m}\mathrm{o}\mathrm{l}/{\mathrm{m}}^3$.

Based on the molecular simulation results, during the fracturing process, water molecules will displace the originally adsorbed methane (CH₄). Substantially, CH₄ molecular concertation, $C_a$ at different pressures is described as follows:

$$C_a=\frac{4px}{\pi N_A{d}_m^3\left(P_L+P\right)\Omega \left(L\right)}$$

(44)

where, $N_A$ is the Avogadro constant, $6.02\times {10}^{23}/\mathrm{m}\mathrm{o}\mathrm{l}$, $P_L$ is the Langmuir pressure, Pa.

Corresponding diffusion coefficient with fracturing impact is described as follows:

$$D_a={D_a}^0\frac{px}{{(P}_L+P\left)\Omega \right(L)}$$

(45)

where, ${D_a}^0$ is the diffusion coefficient from experiment, x and $\Omega \left(L\right)$ are the imbibition length in capillary and tortuosity-considered capillary length, respectively.

3.3.3. Correlation of Free Gas Flux from Slip Flow and Knudsen Diffusion

The movement of free gas within a capillary tube involves continuous flow because of molecular collisions and Knudsen flow caused by molecular-wall collisions. At nano-scale, the flow regime is significantly influenced by the Knudsen diffusion phenomenon; researchers consider this among the first boundary conditions in fluid dynamics. The viscous flow of free gas is impacted by the corresponding Knudsen diffusion boundary [18]. The slip phenomenon within the Knudsen layer will indeed impact the viscous flow of free gas.

In this scenario, the additional free gas viscous flow influenced by the slip condition in an equivalent granular-shaped capillary bundle, acting as a substitute for OMP and brittle minerals pores, could be expressed as follows [72]:

$$N_S=\frac{p}{ZRT}\frac{F{K}_v}{{\mu }_g}\mathsf{\nabla }P$$

(46)

$$F=4\left(\frac{2}{f}-1\right)K_n=\sqrt{\frac{9\pi ZRT}{M}}\frac{u_g}{p{R}_f}\left(\frac{2}{f}-1\right)$$

(47)

where, $M$ is the molecular mass, g/mol;$f$ is the tangential momentum accommodation coefficient (TMAC), dimensionless; and $R_f$ is the pore spacing, m.

Additional free-gas flux velocities from rectangular equivalent slit-shaped clay minerals from X, Y axes directions are expressed as:

$$v\left(x\right)=\frac{\mathsf{\nabla }P}{u_g}\frac{h^2}{12}\frac{2-f}{f}\frac{6K_n}{1-b{K}_n}\sum _{n=1,3,5\dots }^{\infty }\frac{8}{{\pi }^2{n}^2}\left(1-\frac{\mathit{cosh}\left(n\pi x/h\right)}{\mathit{cosh}\left(n\pi \epsilon /2\right)}\right)$$

(48)

$$v\left(y\right)=\frac{\mathsf{\nabla }P}{u_g}\frac{h^2}{12}\frac{2-f}{f}\frac{6K_n}{1-b{K}_n}\sum _{n=1,3,5\dots }^{\infty }\frac{8tanh\left(n\pi \epsilon /2\right)sin\left(n\pi y/h\right)}{{\pi }^2{n}^2}$$

(49)

where, $\epsilon$ is the shape factor of clay minerals, width/height, dimensionless.

The Knudsen diffusion refers to the collision between the CH₄ molecules and the pore wall. The diffusion process is recognised as the secondary boundary, which will also affect the free-gas flux. To characterise this phenomenon, studies have used the diffusion coefficient, $D_k$ from a linear algorithm [45].

The Knudsen diffusion-induced free-gas flux is expressed as follows:

$$N_k=\frac{1}{ZRT}{D}_k=\frac{1}{ZRT}\frac{2R_f}{3}\sqrt{\frac{8ZRT}{\pi M}}\mathsf{\nabla }P$$

(50)

where, $N_k$ is the gas flux from the Knudsen diffusion.

The additional flux due to the Knudsen layer affecting slip-gas flow could be expressed as:

$$N=\frac{p}{Z{\mu }_{g0}RT}\frac{{R_f}^2}{8}\left(1+\frac{2-f}{f}\frac{4K_n}{1-b{K}_n}\right)\mathsf{\nabla }P$$

(51)

3.4. Development of Dynamic Apparent Permeability Model

Above mentioned work in this study have already delineated the equivalent capillary bundle model for describing PSD, the water-film thickness and corresponding water-bridging mechanism influenced by fracturing, dynamic poro-mechanism, and the flow-velocity model for each capillary. Subsequent work will consolidate these mechanisms and models to extract dynamic permeability at varying pressures.

Combining Equations (42), (46) and (50), the free-gas flux considering slip flow at different granular pores is described as follows:

$$\begin{gathered} N_f=\int \pi v\left(r\right)dr+N_k+N_s \\ =\pi \frac{p\left(1+\alpha K_n\right)}{Z{\mu }_{g0}RT}\left[\left(1+\frac{2-f}{f}\frac{4K_n}{1-b{K}_n}+D_a{C}_a\frac{16M}{\pi \ast {\mu }_{g0}{K}_n^2}\right)\frac{R_f^2}{8}\right]\mathsf{\nabla }P\frac{A}{L} \end{gathered}$$

(52)

$$N_{ads}=v_d\frac{A}{L}=D_d{C}_a\frac{ZRT}{P^2}\mathsf{\nabla }P\frac{A}{L}$$

(53)

Based on the Darcy equation, the permeability model for a single capillary is described as follows:

$$Q=\frac{p}{ZRT}\frac{K}{{\mu }_{g0}}A\frac{\Delta P}{L}$$

(54)

Thus, the apparent permeability at the corresponding radius is described as follows:

$$K_{app}=\left(1+\alpha K_n\right)\left[\left(1+\frac{2-f}{f}\frac{4K_n}{1-b{K}_n}+D_a{C}_a\frac{16M}{\pi \ast {\mu }_{g0}{K}_n^2}\right)\frac{R_f^2}{8}\right]$$

(55)

where, $R_f$ is the equivalent pore spacing. Combining Equation (12) on stress depends correlation and Equation (26) on water-film thickness, $R_f$ is express as follows:

$$\begin{gathered} R_f\left(brittle\right)=\Lambda {\left(r\right)}_{iom}\Phi \left(r\right) \\ =\sqrt{{\left[r-h\left(T,V_m,P_v/P_0\right)\right]}^2\left(1-e^{a{\sigma }_m^{\prime }}\right)+{\left[r-h\left(T,{\mathrm{V}}_{\mathrm{m}},P_v/P_0\right)\right]}^2{e}^{a{\sigma }_m^{\prime }}} \end{gathered}$$

(56)

$$\begin{gathered} R_f\left(OMP\right)=\Lambda {\left(r\right)}_{om}\Phi \left(r\right) \\ =\left[r-h\left(T,V_m,P_v/P_0\right)\right]\sqrt[6]{e^{-3c_f\frac{v}{1-v}\left(p-p_i\right)+\frac{E}{3\left(1-v\right)}{\epsilon }_l\left(\frac{p}{p+p_{\epsilon }}-\frac{p_i}{p_i+p_{\epsilon }}\right)}} \end{gathered}$$

(57)

Combining Equations (31) and (32), the water-bridging mechanism at varying pressures is expressed as follows:

$${\left.\frac{\partial \left(\frac{r}{r-h}\cdot {\Pi }_{flat}\left(h\right)+\frac{\gamma }{r-h}\right)}{\partial h}\right|}_{h=h^{*}}=0$$

(58)

The water-bridging mechanism at different pore pressures is described as follows:

$$\Gamma \left(r,P\right)=\left\{\begin{array}{c}1\hspace{1em}P>P_{thw}\\ 0\hspace{1em}P<P_{thw}\end{array}\right.$$

(59)

where, combining Equations (32), (33) and (34), the $P_{threshold}$ is described as follows:

$$P_{thw}=\frac{1}{c}7.0746-\left(1657.46/\left(T-273.15+227.02\right)\right)e^{\left(\frac{r}{r-h*}\cdot {\Pi }_{flat}\left(h*\right)+\frac{\gamma }{r-h*}\right)}$$

(60)

Based on Equations (27), (28) and (29), the water blockage from the fracturing process is expressed as follows:

$${\rm T}\left(r,P\right)=\left\{\begin{array}{c}1\hspace{1em}P>P_{thb}\\ 0\hspace{1em}P<P_{thb}\end{array}\right.$$

(61)

$$P_{thb}=\left\{\begin{array}{c}\frac{P_{g0}{L}_0}{L_0-x_0}+\frac{4\sigma \mathit{cos}{\theta }_{omp}}{d}+\frac{4ZRT{n}_{omp}x\left(p\right)}{d\left(L_0-x\left(p\right)\right)}\\ \frac{P_{g0}{L}_0}{L_0-x_0}+\frac{4\sigma \mathit{cos}{\theta }_b}{d}+\frac{4ZRT{n}_{b}x\left(p\right)}{d\left(L_0-x\left(p\right)\right)}\\ \frac{P_{g0}{L}_0}{L_0-x_0}+\frac{4\sigma \mathit{cos}{\theta }_c}{d}+\frac{4ZRT{n}_{c}x\left(p\right)}{d\left(L_0-x\left(p\right)\right)}+n{E}_{\pi }RT\left(C_0-C_f\right)\end{array}\right.$$

(62)

Thus, the dynamic permeability model of brittle minerals and OMP at varying pressures is derived as follows:

$$\begin{gathered} K\left(P\right)=\sum _{r=1,2\dots }\frac{1}{rlg\sqrt{2\pi }\sigma }\mathit{EXP}\left[-\frac{{\left(lgr-\mu \right)}^2}{2{\sigma }^2}\right] \\ \left(1+\alpha K_n\right)\left[\left(1+\frac{2-f}{f}\frac{4K_n}{1-b{K}_n}+D_a{C}_a\frac{16M}{\pi {\mu }_{g0}{K}_n^2}\right)\frac{R_f^2}{8}\right]\Gamma \left(r,P\right){\rm T}\left(r,P\right) \end{gathered}$$

(63)

Combining Equations (41), (48) and (49), the free-gas flux considering the slip flow at clay pores is described as follows:

$$N_{freegas}=\oint v\left(x,y\right)+\int v\left(x\right)+\int v\left(y\right)+N_{ads}$$

(64)

$$\begin{gathered} N_{freegas}=\oint \frac{\mathsf{\nabla }P}{u_g}\frac{h^2}{12}\sum _{n=1,3,5\dots }^{\infty }\frac{48}{{\pi }^3{n}^3}\left(1-\frac{cosh\left(n\pi x/h\right)}{cosh\left(n\pi w/2\right)}\right)sin\left(n\pi y/h\right) \\ +\int \frac{\mathsf{\nabla }P}{u_g}\frac{h^2}{12}\frac{2-f}{f}\frac{6K_n}{1-b{K}_n}\sum _{n=1,3,5\dots }^{\infty }\frac{8}{{\pi }^2{n}^2}\left(1-\frac{cosh\left(n\pi x/h\right)}{cosh\left(n\pi \epsilon /2\right)}\right) \\ +\int \frac{\mathsf{\nabla }P}{u_g}\frac{h^2}{12}\frac{2-f}{f}\frac{6K_n}{1-b{K}_n}\sum _{n=1,3,5\dots }^{\infty }\frac{8tanh\left(n\pi \epsilon /2\right)sin\left(n\pi y/h\right)}{{\pi }^2{n}^2}+D_d{C}_a\frac{ZRT}{P^2}\mathsf{\nabla }P \end{gathered}$$

(65)

$$N_{freegas}=wh\frac{\mathsf{\nabla }P{h}^2}{12{\mu }_g}\left(1+\alpha K_n\right)F$$

(66)

$$\begin{gathered} F=1-\frac{192}{{\pi }^5\epsilon }\sum _{n=1,3,5\dots }^{\infty }\frac{tanh\left(n\pi \epsilon /2\right)}{n^5}+\frac{2-f}{f}\frac{6K_n}{1-b{K}_n}\sum _{n=1,3,5\dots }^{\infty }\frac{8}{{\pi }^2{n}^2}\left[\left(1-\frac{2tanh\left(n\pi \epsilon /2\right)}{n\pi \epsilon }\right)\right] \\ +\frac{2-f}{f}\frac{3K_n}{1-b{K}_n}\sum _{n=1,3,5\dots }^{\infty }\frac{32tanh\left(n\pi \epsilon /2\right)}{{\pi }^3{n}^3}+D_a{C}_a\frac{48M}{\pi {\mu }_g}{K}_n^2 \end{gathered}$$

(67)

Corresponding permeability on clay minerals pores is described as follows:

$$\begin{gathered} K\left(P\right)=\sum _{r=1,2\dots }\frac{1}{rlg\sqrt{2\pi }\sigma }\mathit{EXP}\left[-\frac{{\left(lgr-\mathsf{\mu }\right)}^2}{2{\sigma }^2}\right]\frac{h^2}{12}\left(1+\alpha K_n\right)\{1 \\ -\frac{192}{{\pi }^5\epsilon }\sum _{n=1,3,5\dots }^{\infty }\frac{tanh\left(n\pi \epsilon /2\right)}{n^5} \\ +\frac{2-f}{f}\frac{6K_n}{1-b{K}_n}\sum _{n=1,3,5\dots }^{\infty }\frac{8}{{\pi }^2{n}^2}\left[\left(1-\frac{2tanh\left(n\pi \epsilon /2\right)}{n\pi \epsilon }\right)\right] \\ +\frac{2-f}{f}\frac{3K_n}{1-b{K}_n}\sum _{n=1,3,5\dots }^{\infty }\frac{32tanh\left(n\pi \epsilon /2\right)}{{\pi }^3{n}^3} \\ +D_a{C}_a\frac{48M}{\pi {\mu }_g}{K}_n^2\}\Gamma \left(r,P\right){\rm T}\left(r,P\right) \end{gathered}$$

(68)

4. Model Validation

It is difficult to validate the proposed model with the imbibition process via experiments because of the significant pressure difference during hydraulic fracturing or the nano-scaled pore structure in shale, in this study, the proposed model is separated into two parts for validation. First, we conducted a comparison with the conventional dynamic permeability model and a porous-meter experiment to examine the righteousness of the PSD-coupled dynamic permeability model. The proposed model was then applied in MRST-shalemaster—an open-sourced EDFM simulator, to assess its accuracy while using conventional intrinsic permeability and the proposed dynamic permeability model.

4.1. Validation of Apparent Permeability Model for Single Capillary

To validate the permeability model in an experiment-viable condition, this study compared (Figure 16a) the permeability against a single capillary tube based on Sheng (2019) and Wang (2022) [24,66]. These two methods explain the permeability at a single capillary tube and originate from Javadpour’s [47] study and the Poiseuille theory, respectively. Figure 16 illustrates that the proposed capillary model is in good agreement with both models. Sheng’s work [24] focused on the influence on the shape factor, ignored the wettability on the Knudsen flow while Wang’s work [66] consider the constant film thickness. Javadepour-based models, compared to the Poiseuille-based AP model showed a relatively higher number in permeability, both AP models showed that the shape factor and water film thickness can be ignored at larger-scaled pores. In addition, our model uses equivalent round pore radius to describe the clay minerals while Sheng [24] uses height and shape factors to describe PSD in ellipse pores, making the curve of eclipse mineral pores align with granular pores and ignoring the PSD in the real situation, this interpretation will underestimate the influenced of shape factor, from our work, using equivalent radius interpretation by NMR data for describing slit-liked clay pores, limited the total flow area, shape factor is then becoming a dominating parameters in AP model.

4.2. Validation of Dynamic Permeability with Experiment Data

Experimental validation of the PSD-integrated dynamic permeability is conducted under the same shale sample extracted from the WY formation after the NMR tests. The cylindrical core sample is 2.5 cm in diameter with 5 cm in length. Permeability was tested under the pulse gradient porous meter with CH₄ as injecting media to maintain pressure from 1.2 Mpa to 5 Mpa. Owing to the limitation of the experiment setup, this part of the numerical model neglected the imbibition-induced water blockage and desorption phenomenon; the input parameters of this part are provided in Table 1. The comparison of the model’s calculation and experimental data is shown in Figure 16b.

As previously noted (Figure 2), studies from different samples showed different trends as the pressure increased. Zeng [73] suggested that the two main contributions to the change in permeability include the stress-induced matrix shrinkage and gas-slippage weakening, samples with nature fractures trends to reflect a “U-shape” in the dynamic permeability while others showed a declining trend as the pressure increases [74]. To validate the proposed dynamic permeability model, this part takes the comparison with the stress-dependence poroelastic approach and the gas-slippage dominated AP model, which is recognised as the upper boundary and lower boundary of the dynamic permeability model, respectively [73]. Our numerical model is situated in between both models because the proposed model takes into account both dynamic flow regime and poroelastic properties. Comparison with numerical model and experiment data also matches well.

4.3. Validation of Dynamic Permeability under Field Conditions: Numerical Simulation

It is challenging to reflect fracturing-influenced field conditions. To validate the proposed dynamic permeability during production, we made amendments to the permeability module in an open-source EDFM simulator based on the Matlab Reservoir Simulation Toolkit (MRST) [75]. The results of the amended EDFM with the original EDFM using the WY1 field production data are used to validate the overall dynamic permeability, and the relevant SGR parameters are presented in Table 1.

The WY 1 platform employs a long-section, multi-cluster fracturing with segmented perforation mode. The main segments have lengths ranging from 110 m to 135 m, averaging 120 m. Following the principle of segmenting without crossing small layers and choosing segments with similar physical properties, a single segment consists of 10 perforation clusters, each cluster having a length of 0.35. For longer segments (exceeding 130 m), a single segment consists of 11 clusters. In segments with a risk of fractures, a single segment consists of 12 clusters, there are a total of 15 fracturing segments, with cluster spacing in the range of 4.75 to 10.65 m. According to the fracturing design, each segment can be recognised as isolated. Thus, we set the EDFM simulation into three cases, Case 1 includes 12 hydraulic fractures with 100 m fracture half-length; Case 2 includes 11 hydraulic fractures and 105 m fracture half-length; and Case 3 includes 10 hydraulic fractures with 130 m fracture half-length, related SGR parameters are shown in

Table 3. The main parameters from EDFM and Amended EDFM model.

Parameters	EDFM	Amended EDFM
Matrix permeability (nD)	600	Dynamic
Fracturing pressure (Mpa)	-	50
Porosity	0.03
Fracture permeability (D)	0.1
Natural fracture perm (D)	0.01
Bottom hole pressure (MPa)	15
Initial pore pressure (MPa)	23
Case1 contribution	0.2
Case2 contribution	0.7
Case3 contribution	0.1

. The simulation takes into account the randomly generated natural fracture mapping from MRST toolbox, for comparison purposes, the random seed was set as 61,320 to get the extract same natural fracture pattern

Figure 17 illustrates the pressure distribution in three cases, to reflect 15 segments in the field, we use the contribution factor from three cases to reflect the cumulative production (Table 3). The results, in terms of gas production rate and cumulative production, are plotted in Figure 18, including field data that incorporates the shut-in period. The amended EDFM result initially underestimated the cumulative gas production, followed by an overestimation. However, when compared with the EDFM result, it demonstrated higher accuracy in predicting cumulative gas production over time.

5. Result and Discussion

5.1. Dynamic Permeability under Different Reservoir Conditions

Figure 19 illustrates the apparent permeability at varying fracture and pore pressure. Figure 19a presents the dynamic permeability when pore pressure is smaller than the fracture pressure. Figure 19b demonstrates the dynamic permeability when the pore pressure is bigger than the fracture pressure, owing to the capillary phenomenon. First, the graph shows that dynamic permeability decreased gradually, and unlike the experimental result, the trend at high pore pressure did not reach a plateau because the water-bridging phenomenon resulted in water blockage at high pore pressure. Second, at the same pore pressure, the effect of the imbibition process is mostly on the adsorption layer distribution on the pore wall and the corresponding pore-size shrinkage. Third, dynamic porous parameters also play a key role in dynamic permeability, including PSD distribution, water-film thickness, pressure, and temperature during production. The above-mentioned phenomenon will be discussed in detail in the following section to address the field suggestions based on dynamic permeability.

5.2. Fracturing-Induced Water Blockage and Water-Film Thickness during Production

Reservoir parameters used in this section are listed in Table 1 to describe the contribution along the pore radius. A schematic demonstration of permeability contribution at varying pore radii is presented in Figure 20, which is, the PSD ratio times the corresponding permeability. Four reservoir pressure conditions (29 Mpa, 23 Mpa, 22 Mpa, and 21 MPa) were compared to describe the fracturing-induced water-blockage phenomenon. We can conclude that, after the hydraulic fracturing process, as the reservoir depleted, hydrophobic OMP started to exhaust and contribute to the permeability owing to its ultra-strong capillary forces. As for hydrophilic minerals, extra capillary forces are required for water-blockage removal because of its wettability. Additionally, the water film mechanism inside pores will escalate the capillary force by changing the equivalent pore radius and substantially impact the permeability owing to its high concertation in OMP. However, the water-bridging mechanism will occur only under 50 nm pores under reservoir conditions.

Figure 20 demonstrates the relationship of three different minerals in shale at varying pressures and temperatures, within three different types of minerals. Note that in the permeability system, smaller pores (1 nm–200 nm) are more affected by the water-film phenomenon. Compared to brittle minerals and OMP, pores within clay minerals are more affected because of their slit-shaped structure. Additionally, Figure 21 illustrates the influences on pore pressure and reservoir temperature. Comparison of the permeability model conducted under 20 MPa and 23 Mpa at 383 K showed that greater pore pressure will have a heavier impact on permeability. The comparison between 440 K and 383 K at 23 Mpa pore pressure concluded that higher temperature will have a greater impact on permeability. Also, for a different type of shale rocks, the overall wettability of brittle minerals, clay and OMP at different TOC will have an influence on the waterblockage phenomenon at different pressure, smaller contact angles will escalate the waterblockage phenomenon, consequentially reducing the overall permeability.

5.3. Fracturing Effected Adsorption Layer Thickness and Consequentially Impacted Permeability

The fracturing response has a significant impact on the poro-mechanical of the pore. Figure 22a demonstrates the change of diffusion coefficient under different fracturing pressures. Figure 22a utilizes the 1 minus diffusion coefficient to at different fracturing pressure to illustrate the impact on the diffusion phenomenon. From the curve, we can see that the fracturing process will significantly impact the diffusion coefficient. The fracturing impact ratio increases as the fracturing pressure increases. For hydrophilic brittle and clay minerals, under identical fracturing pressure, small pores undergo greater diffusion coefficient damage than bigger pores owing to their wettability. For hydrophobic OMPs, smaller pores undergo less diffusion coefficient damage owing to strong capillary at smaller pores.

Apart from the impact on the diffusion coefficient, fracturing-induced thinning of the desorption portion layer also contributes to the desorption-induced matrix shrinkage, directly increasing the permeability contributed by free gas [43,67,72]. Figure 22b illustrates the permeability changes with different fracturing pressures for both free gas and adsorbed gas. The triangle marks with the corresponding colours depict the derived permeability at different radii. Higher fracturing pressure resulted in a reduction of adsorbed gas permeability by around 30–60%, while the incremental in free gas permeability was only at 0.02%.

5.4. Influence of Clay-Shaped Factors

The dynamics of the proposed model are significantly influenced by clay minerals due to their high concentration at a larger pore scale. Besides the concentration of clay minerals, the shape factor (SF) also plays a vital role in overall permeability.

The FIB-SEM description can reflect the shape factor distribution at a larger scale in the matrix (Figure 23). Neglecting the granular-shaped brittle mineral pores contributing at a small SF, from observations, the real SF distribution in slit-like pores also follows a Gaussian distribution. Thus, we present three cases with smaller, medium, and higher SF (7, 20, 43) dominated distributions with the same porosity to discuss the impact on permeability (Figure 24a).

Figure 24b describes the contributed permeability at different radii and overall permeability. From the contributed permeability curve, the permeability at different radii with a small shape factor has a larger number. Overall permeability is reduced by 70% when comparing Case 1 (mean SF = 7) and Case 2 (mean SF = 20), but reduced by 83% when compared with Case 3 (mean SF = 43). The increment of SF and the corresponding overall permeability follows a logarithmic relationship.

6. Conclusions

In this study, we developed the entire workflow to extract dynamic permeability under the field conditions impacted by the fracturing process, The dynamic permeability model explains the water-gas interaction in the shale matrix, which is perfect for explaining the gas flow during the reservoir depletion and well-testing period to prevent water blockage with the influence of hydraulic fracturing. For field application, this methodology could be used to assist in setting a development plan for enhancing EUR during the early period by controlling the bottom-hole pressure to optimise the matrix permeability. Also, this model demonstrated dynamic permeability in the shale matrix at different reservoir pressures. This article also demonstrated that the proposed dynamic permeability model, combining hydraulic fracture and natural fracture with the EDFM showed higher accuracy in predicting the production rate.

Key findings in this study include:

• Bayesian-assisted Gaussian description for the three majority minerals derived from coupled FIB-SEM and NMR data proves to be viable for accurately describing the PSD in shale. The corresponding dynamic permeability model demonstrates an intimate association with experimental data.
• The fracturing-induced imbibition process results in water blockage and water bridging mechanisms during shale gas production, impacting the dynamic permeability in the matrix. The water blockage phenomenon significantly reduces permeability in nano-scaled brittle minerals and clay. Substantial water blockage requires a larger pressure gradient to overcome. The impact on total permeability, however, depends on the subsequent PSD in the micro-scale. Nano-scaled dominated OMP and micro-scaled dominated brittle and clay pores reduce the impact on permeability from water blockage. Water bridging occurs only in nano-scale OMP below 50nm at high pressure and temperatures. Due to the high concentration of OMP, the permeability contribution cannot be neglected.
• Reservoir depletion has a substantial impact on permeability, showing a declining trend as pore pressure increases. In addition to poromechanics, the fracturing-induced imbibition phenomenon reduces the thickness of the water film inside pores, significantly impacting adsorbed gas permeability, and only marginally boosting the contribution of free gas to permeability.
• Due to the high permeability of large-scale clay mineral pores, the SF of clay minerals significantly influences dynamic permeability. With constant pore spacing, higher SF clay mineral pores reduce permeability in the shale matrix.