Matrix–Fracture Interactions During Flow in Organic Nanoporous Materials Under Loading

Abstract

The transport of natural gas in shale resources is of multiscale and multiphysics nature. The gas transport involves dealing with pores of different length scales and geometrical orientation. Understanding and modeling the interactions of those pores are crucial for proper description of shale matrix deliverability and dynamics. These interactions are less understood at a microscopic scale during the depletion of shale reservoir. In this study, we consider transient flow behavior of a compressible fluid in organic nanoporous material with micro-fractures, or cracks, using a pore network modeling approach. We present a multistep workflow where the transport in the organic nanopores is studied considering the advective–diffusive–adsorptive mechanisms then, coupled with the microcracks for better understanding and optimizing of natural gas organic-rich shale. The natural gas is initially stored in the pore network at high pressure as free and adsorbed fluids and its pressure-driven viscous flow includes additional diffusion mechanisms. The percolation theory is used to obtain some approximations to the organic matrix flow regime coefficients. The organic materials–microcracks coupling term is derived and validated relating the flow rate exchange to the pressure difference with some modifications in order to account for the organic matrix transport dependency on the pressure. The coupling flow exchange term is used to link the local nanoscale heterogeneity to the large scale continuum modeling of shale reservoir. The upscaled model captures the transport exchange during the transient and steady-state flow conditions.

Appendices

Appendix A: Discretization of the Mass Balance in the Pore Network Model:

The calculations around each node (pore) are done by the material balance equation given in Eq. 1. Then, the discretized form is illustrated by taking a single pore interconnected with two throats as shown in Fig. 12.

To avoid complexity, the fluid velocity in the discretized form is taken as the Hagen–Poiseuille velocity. The same procedure can be done with the modified Hagen–Poiseuille form.

$$\begin{aligned}&\displaystyle A_\mathrm{in} \left( {u\rho } \right) _\mathrm{in} -A_\mathrm{out} \left( {u\rho } \right) _\mathrm{out} +\frac{\partial \left( {\rho V^{*}} \right) _\mathrm{pore} }{\partial t}=0 \end{aligned}$$

(21)

$$\begin{aligned}&\displaystyle \rho =\frac{PM}{zRT} \end{aligned}$$

(22)

$$\begin{aligned}&\displaystyle u=\frac{R^{2}}{8\mu }\frac{\partial P}{\partial x} \end{aligned}$$

(23)

$$\begin{aligned}&\displaystyle A=\pi R^{2} \end{aligned}$$

(24)

Then, the mass flux is:

$$\begin{aligned} u\rho =\frac{MR^{2}}{8RT}\left( {\frac{P}{\mu z}} \right) \frac{\partial P}{\partial x}=\frac{MR^{2}}{16RT\left( {\mu z} \right) _\mathrm{avg} }\frac{P_1^2 -P_2^2 }{L} \end{aligned}$$

(25)

Mass accumulation:

$$\begin{aligned} \frac{\partial \left( {\rho V} \right) _\mathrm{pore} }{\partial t}=\frac{\frac{M}{RT}\left[ {\left( {\frac{PV^{*}}{z}} \right) _{t1} -\left( {\frac{PV^{*}}{z}} \right) _{t2} } \right] }{\Delta t} \end{aligned}$$

(26)

Continuity equation becomes:

$$\begin{aligned}&\frac{M\pi R_\mathrm{in}^{4}}{16RT\left( {\mu z} \right) _\mathrm{avg, in} }\frac{P_\mathrm{in}^2 -P_\mathrm{pore}^2 }{L}-\frac{M\pi R_\mathrm{out}^{4}}{16RT\left( {\mu z} \right) _\mathrm{avg,out} }\frac{P_\mathrm{pore}^2 -P_\mathrm{out}^2 }{L}\nonumber \\&\quad +\,\frac{M}{RT}\frac{\left[ {\left( {\frac{P_\mathrm{pore} V^{*}}{z}} \right) _{t1} -\left( {\frac{P_\mathrm{pore} V^{*}}{z}} \right) _{t2} } \right] }{\Delta t}=0 \end{aligned}$$

(27)

where $V^{*}$ is the pore volume corrected for the presence of the adsorbed layer using Langmuir adsorption model:

$$\begin{aligned} V^{*}=V\left( {1-\frac{V_\mathrm{s} }{V}} \right) ,\quad V_\mathrm{s} =V_\mathrm{sL} \frac{P}{P_\mathrm{L} +P} \end{aligned}$$

(28)

Rearrangement and simplification of the equation:

$$\begin{aligned}&\left[ {\frac{\pi R_\mathrm{out}^{4}}{16L\left( {\mu z} \right) _\mathrm{avg,out} }+\frac{\pi R_\mathrm{in}^{4}}{16L\left( {\mu z} \right) _\mathrm{avg,\,in} }} \right] P_{\mathrm{pore}_{t2}}^2 +\frac{1}{\Delta t}\left( {\frac{P_\mathrm{pore} V^{*}}{z}} \right) _{t2}\nonumber \\&\quad =\frac{\pi R_\mathrm{in}^{4}}{16L\left( {\mu z} \right) _\mathrm{avg,\,in} }P_\mathrm{in}^2 +\frac{\pi R_\mathrm{out}^{4}}{16L\left( {\mu z} \right) _\mathrm{avg,out} }P_\mathrm{out}^2 +\frac{1}{\Delta t}\left( {\frac{P_\mathrm{pore} V^{*}}{z}} \right) _{t1} \end{aligned}$$

(29)

System of nonlinear equations is to be solved for pressure at each pore as a function of time step $\Delta {t}$.

$$\begin{aligned}&aP_\mathrm{m}^{2}+\hbox {b}P_\mathrm{m} +cP_{vS}^{2}+dP_\mathrm{N}^{2}+eP_\mathrm{NE}^{2}+fP_\mathrm{NW} ^{2}+gP_\mathrm{SE}^{2}+hP_\mathrm{SW}^{2}+iP_\mathrm{UN}^{2}\nonumber \\&\quad +\,jP_\mathrm{DS}^{2}+kP_\mathrm{UNE} ^{2}+lP_\mathrm{UNW}^{2}+\cdots +c1=0 \end{aligned}$$

(30)

a, b, c, d, e ... are all function of pressure and that nonlinear system of equations can be solved using the iterative Newton–Raphson method:

$$\begin{aligned} P_{\mathrm{n},t+1} =P_{\mathrm{n},t} +\frac{f\left( {P_{\mathrm{n},t} } \right) }{{f}^{\prime }\left( {P_{\mathrm{n},t} } \right) } \end{aligned}$$

(31)

For the multiscale network, the velocity in the crack is given by Eq. (12).

Appendix B: The Use of Crank Sphere Model in Defining the Transient Factor

On our case, the Crank sphere model cannot be used directly as given in Zimmerman et al. (1993) because of the pressure dependence of the organic material permeability.

$$\begin{aligned} \frac{P_\mathrm{om} -P_i }{P_\mathrm{f} -P_\mathrm{i} }=\left[ {1-\exp \left( {-Ckt} \right) } \right] ^{1/2} \end{aligned}$$

(32)

One way to linearize the problem is to replace the pressure with a pressure function $\square$. In this case, the permeability k would be the intrinsic or liquid permeability.

$$\begin{aligned} \frac{\psi (P_\mathrm{om} )-\psi (P_\mathrm{i} )}{\psi \left( {P_\mathrm{f} } \right) -\psi (P_\mathrm{i} )}=\left[ {1-\exp \left( {-Ckt} \right) } \right] ^{1/2} \end{aligned}$$

(33)

Rearrangement:

$$\begin{aligned} \left( {1-\left( {\frac{\psi (P_\mathrm{om} )-\psi (P_\mathrm{i} )}{\psi \left( {P_\mathrm{f} } \right) -\psi (P_\mathrm{i} )}} \right) } \right) ^{2}=\exp \left( {-Ckt} \right) \end{aligned}$$

(34)

Taking the log of both sides, then the plot of the right-hand side versus the left-hand side in the Cartesian coordinate should be a straight line with y-intercept $=$ 0. That is validated using the data obtained from the network as shown in Figs. 13 and 14.

Fig. 13

Schematic of single pore interconnected with two throats to explain how the material balance is discretized for the pore network model

Fig. 14

Validation of the modified Crank sphere model using the results obtained from the network model