The Dual-Reciprocity Boundary Element Analysis for Hydraulically Fractured Shale Gas Reservoirs Considering Diffusion and Sorption Kinetics

Abstract

This work presents a new application of boundary element method (BEM) to model fluid transport in unconventional shale gas reservoirs with discrete hydraulic fractures considering diffusion, sorption kinetics and sorbed-phase surface diffusion. The fluid transport model consists of two governing partial differential equations (PDEs) written in terms of effective diffusivities for free and sorbed gases, respectively. Boundary integral formulations are analytically derived using the fundamental solution of the Laplace equation for the governing PDEs and Green’s second identity. The domain integrals arising due to the time-dependent function and nonlinear terms are transformed into boundary integrals employing the dual-reciprocity method. This transformation retains the domain-integral-free, boundary-integral-only character of standard BEM approaches. In the proposed solution, the free- and sorbed-gas flow in the shale matrix is solved simultaneously after coupling the fracture flow equation of free gas. Well production performance under the effect of relaxation phenomenon due to delayed responses of sorbed gas under nonequilibrium sorption condition is rigorously captured by imposing the zero-flux condition at fracture–matrix interface for the sorbed-gas transport equation. The validity of proposed solution is verified using several case studies through comparison against a commercial finite-element numerical simulator.

Appendices

Appendix 1: Derivation of the Dimensionless Equation and Variables

The continuity equation for free-gas flow in matrix, as flux to fracture counted as source term, is expressed as:

$$\nabla \cdot \left( {D_{{{\text{eff}}}} \nabla C} \right) - Q_{m} = \frac{{\partial \left( {\phi C} \right)}}{\partial t} - R_{s} .$$

(33)

Q_m is the strength distribution of the source in mass per unit volume of the porous medium per unit time. Following similar development in Liggett and Liu (1983) and Pecher and Stanislav (1997), we treat N_p discrete fractures as plane sources (Γp_j) and rewrite Q_m as follows:

$$Q_{m} = \sum\limits_{j = 1}^{{n_{p} }} {\frac{1}{{\Gamma_{p,j} }}\int\limits_{{\Gamma_{p,j} }} {Q_{p,j} \delta \left( {x - x_{p} } \right)\delta \left( {y - y_{p} } \right)} } d\tau$$

(34)

where strength Q_p,j is a function of position along plane sources Γ_p,j. Using dimensionless group provided in Table

Table 2 The dimensionless variables used in the proposed DRBEM formulation

2, the dimensionless form of Eq. 33, after linearization, is given as:

$$\nabla^{2} C_{{\text{D}}} - \sum\limits_{k = 1}^{{N_{p} }} {\frac{1}{{\Gamma_{p,k} }}\int\limits_{{\Gamma_{p,k} }} {Q_{pD,k} \delta \left( {x_{{\text{D}}} - x_{pD} } \right)\delta \left( {y_{{\text{D}}} - y_{{p_{{\text{D}}} }} } \right)} } d\tau_{{\text{D}}} = \frac{{\partial C_{{\text{D}}} }}{{\partial t_{{\text{D}}} }} - R_{sD} .$$

(35)

Similarly, for the sorbed-gas equation:

$$\nabla \cdot \left[ {\left( {1 - \phi } \right)D_{{\text{s}}} \nabla C_{{\text{s}}} } \right] = \frac{{\partial \left[ {\left( {1 - \phi } \right)C_{{\text{s}}} } \right]}}{\partial t} + R_{s} .$$

(36)

Substituting the dimensionless group in Table 2 leads to:

$$\nabla^{2} C_{sD} = \frac{{\partial C_{sD} }}{{\partial t_{sD} }}{ + }R_{sD} .$$

(37)

Appendix 2: The Dual-Reciprocity Boundary Element Method for Gas Flow in Matrix

Multiplication of Eq. 8 by 13 leads to:

$$G\nabla^{2} C_{{\text{D}}} - \sum\limits_{k = 1}^{{N_{p} }} {\frac{1}{{\Gamma_{p,k} }}\int\limits_{{\Gamma_{p,k} }} {GQ_{pD,j} \delta \left( {x_{{\text{D}}} - x_{pD} } \right)\delta \left( {y_{{\text{D}}} - y_{{p_{{\text{D}}} }} } \right)} } d\tau_{{\text{D}}} = \sum\limits_{i = 1}^{{N_{t} }} {Gf_{i} \left( {x_{{\text{D}}} ,y_{{\text{D}}} } \right)} \alpha_{i} \left( {t_{{\text{D}}} } \right).$$

(38)

Integration of Eq. 38 over the entire domain ${\Omega }$ leads to

$$\int\limits_{\Omega } {G\nabla^{2} C_{{\text{D}}} } d\Omega - \sum\limits_{k = 1}^{{N_{p} }} {\frac{1}{{\Gamma_{p,k} }}\int\limits_{{\Gamma_{p,k} }} {GQ_{pD,j} } } d\tau_{{\text{D}}} = \sum\limits_{i = 1}^{{N_{t} }} {\alpha_{i} \left( {t_{{\text{D}}} } \right)\int\limits_{\Omega } {G\nabla^{2} \hat{C}_{{{\text{D}},i}} } d\Omega } .$$

(39)

Recall Green’s second identity:

$$\int\limits_{\Omega } {\left( {G\nabla^{2} C_{{\text{D}}} - C_{{\text{D}}} \nabla^{2} G} \right)} d\Omega = \int\limits_{\Gamma } {\left( {G\frac{{\partial C_{{\text{D}}} }}{\partial n} - C_{{\text{D}}} \frac{\partial G}{{\partial n}}} \right)} d\tau$$

(40)

where n is the outward normal vector of the boundary Γ. Because G is complementary solution of the homogenous counterpart of Eq. 4, which represents the system response to the unit source strength at (ξ, η):

$$\nabla^{2} G = - \delta \left( {x_{{\text{D}}} - \xi } \right)\delta \left( {y_{{\text{D}}} - \eta } \right).$$

(41)

Using the shifting property of Dirac delta function, for a point C_D(x_D,y_D) that is inside domain ${\Omega }$, the integral of Dirac delta function over region ${\Omega }$ is evaluated over the integral region of $\Gamma_{ \in }$ where $\in \to 0$; if a point C_D(x_D,y_D) is on the smooth boundary Γ, the integral region of $\Gamma_{ \in }$ is divided into two parts, and only half of integration area ($\in^{\prime }$) should be evaluated. Therefore, the second term on the LHS of Eq. 40 is:

$$\int\limits_{\Omega } {\left( { - C_{{\text{D}}} \nabla^{2} G} \right)} d\Omega = \theta C_{{\text{D}}} .$$

(42)

θ equals to one if the C_D is inside domain, and 0.5 if it is on smooth boundary. Combining Eqs. 39, 40 and 42, the LHS of Eq. 39 is now transformed to boundary integral form as:

$$\begin{gathered} G\nabla^{2} C_{{\text{D}}} - \sum\limits_{k = 1}^{{N_{p} }} {\frac{1}{{\Gamma_{p,k} }}\int\limits_{{\Gamma_{p,k} }} {GQ_{pD,k} \delta \left( {x_{{\text{D}}} - x_{pD} } \right)\delta \left( {y_{{\text{D}}} - y_{{p_{{\text{D}}} }} } \right)} } \hfill \\ = \theta C_{{\text{D}}} + \int\limits_{\Gamma } {\left( {C_{{\text{D}}} \frac{\partial G}{{\partial n}} - G\frac{{\partial C_{{\text{D}}} }}{\partial n}} \right)} d\tau - \sum\limits_{k = 1}^{{N_{p} }} {\frac{1}{{\Gamma_{p,k} }}\int\limits_{{\Gamma_{p,k} }} {GQ_{pD,k} } } d\tau_{{\text{D}}} . \hfill \\ \end{gathered}$$

(43)

Applying Green’s second identity and following similar development from Eqs. 40 to 42, the RHS of is rewritten in terms of boundary integral as:

$$\sum\limits_{i = 1}^{{N_{t} }} {\alpha_{i} \left( {t_{{\text{D}}} } \right)\int\limits_{\Omega } {G\nabla^{2} \hat{C}_{{{\text{D}},i}} } d\Omega } = \sum\limits_{i = 1}^{{N_{t} }} {\left[ {\alpha_{i} \left( {\theta \widehat{C}_{{{\text{D}},i}} - \int\limits_{\Gamma } {\left( {\widehat{C}_{{{\text{D}}_{,i} }} \frac{\partial G}{{\partial n}} - G\frac{{\partial \widehat{C}_{{{\text{D}},i}} }}{\partial n}} \right)} d\tau } \right)} \right]} .$$

(44)

Combination of Eqs. 43 and 44 results in:

$$\theta C_{{\text{D}}} + \int\limits_{\Gamma } {\left( {C_{{\text{D}}} \frac{\partial G}{{\partial n}} - G\frac{{\partial C_{{\text{D}}} }}{\partial n}} \right)} d\tau - \sum\limits_{k = 1}^{{N_{p} }} {\frac{1}{{\Gamma_{p,k} }}\int\limits_{{\Gamma_{p,k} }} {GQ_{pD,k} } } d\tau_{{\text{D}}} = \sum\limits_{i = 1}^{{N_{t} }} {\left[ {\alpha_{i} \left( {\theta \widehat{C}_{{{\text{D}},i}} - \int\limits_{\Gamma } {\left( {\widehat{C}_{{{\text{D}}_{,i} }} \frac{\partial G}{{\partial n}} - G\frac{{\partial \widehat{C}_{{{\text{D}},i}} }}{\partial n}} \right)} d\tau } \right)} \right]} .$$

(45)

Appendix 3: Fracture Flow Model

Assuming fracture has no storage capacity, gas flow in one-dimensional fracture can be written as:

$$\frac{\partial }{\partial l}\left( {C\frac{{k_{{\text{f}}} }}{{\mu_{{\text{g}}} }}\frac{\partial p}{{\partial l}}} \right) - Q_{f}^{*} = 0$$

(46)

where l represents flow direction along 1D fracture, ${ }Q_{f}^{*}$ is the line source strength of the fracture and is opposite in sign to that of matrix. Notably, only viscous (Darcian) flow is considered in fracture domain. Substituting definition of gas compressibility c_g, Eq. 46 can be further approximated as following:

$$\frac{\partial }{\partial l}\left( {\frac{{k_{{\text{f}}} }}{{\mu_{{\text{g}}} c_{{\text{g}}} }}\frac{\partial C}{{\partial l}}} \right) - Q_{f}^{*} = 0.$$

(47)

The dimensionless form of Eq. 47, after employing definitions in "Appendix 1", is written as:

$$\frac{{\text{d}}}{{{\text{d}}l_{D} }}\left( {\lambda_{g} \frac{{dC_{{\text{D}}} }}{{dl_{{\text{D}}} }}} \right) - Q_{fD}^{*} = 0$$

(48)

where

$$\lambda_{g} = \frac{{\mu_{{{\text{g}},{\text{ini}}}} c_{{{\text{g}},{\text{ini}}}} }}{{\mu_{{\text{g}}} c_{{\text{g}}} }}$$

$$l_{fD} = l/\sqrt A$$

$$Q_{fD}^{*} = \frac{{Q_{f}^{*} \mu_{{g,{\text{ini}}}} c_{{g,{\text{ini}}}} A}}{{k_{f} \left( {C_{{{\text{ini}}}} - C_{{{\text{ref}}}} } \right)}}.$$

(49)

$\mu_{{{\text{g}},{\text{ini}}}} c_{{{\text{g}},{\text{ini}}}}$ is the gas viscosity–compressibility product at initial condition. Fracture domain can be discretized into a series of control volumes, and Eq. 48 is integrated over each control volume $\Delta V$:

$$\int\limits_{\Delta V} {\frac{{\text{d}}}{{{\text{d}}l_{D} }}\left( {\lambda_{g} \frac{{{\text{d}}C_{D} }}{{{\text{d}}l_{D} }}} \right)} dV - \int\limits_{\Delta V} {Q_{fD}^{*} } dV = w_{f} hl_{{\text{D}}} \int\limits_{{ - \frac{{ - l_{fD} }}{2}}}^{{\frac{{ + l_{fD} }}{2}}} {\frac{{\text{d}}}{{{\text{d}}l_{{\text{D}}} }}\left( {\lambda_{g} \frac{{{\text{d}}C_{{\text{D}}} }}{{{\text{d}}l_{{\text{D}}} }}} \right)} dl - \frac{{\overline{Q}_{fD}^{*} }}{{h\Delta l_{fD} }} = 0$$

(50)

where

$$\frac{{\overline{Q}_{fD}^{*} }}{{\Delta l_{fD} }} = \int\limits_{{\Delta l_{fD} }} {Q_{fD}^{*} } dl = \frac{{\mu_{{{\text{g}},{\text{ini}}}} c_{{{\text{g}},{\text{ini}}}} A}}{{k_{f} \left( {C_{{{\text{ini}}}} - C_{{{\text{ref}}}} } \right)}}\int\limits_{{\Delta l_{fD} }} {Q_{f}^{*} } dl.$$

(51)

Relating Eq. 50 to the source strength used in matrix flow, Eq. 50 is readily simplified using center difference approximation as:

$$\left. {\lambda_{g} } \right|_{{l_{fD} + \frac{{\Delta l_{fD} }}{2}}} \left( {\frac{{dC_{{\text{D}}} }}{{dl_{fD} }}} \right)_{{l_{fD} + \frac{{\Delta l_{fD} }}{2}}} - \left. {\lambda_{g} } \right|_{{l_{fD} - \frac{{\Delta l_{fD} }}{2}}} \left( {\frac{{dC_{{\text{D}}} }}{{dl_{fD} }}} \right)_{{l_{fD} - \frac{{\Delta l_{fD} }}{2}}} - \frac{{\mu_{{{\text{g}},{\text{ini}}}} c_{{{\text{g}},{\text{ini}}}} D_{{{\text{eff}}}} }}{{k_{{\text{f}}} w_{{\text{f}}} }}q_{fD} = 0.$$

(52)

The generalized form (Eq. 52) that is written for i-th cell is expressed as:

$$T_{i,i - 1} \left( {C_{{{\text{D}},i - 1}} - C_{{{\text{D}},i}} } \right) - T_{i + 1,i} \left( {C_{{{\text{D}},i}} - C_{{{\text{D}},i + 1}} } \right) - \frac{{\mu_{{g,{\text{ini}}}} c_{{g,{\text{ini}}}} D_{{{\text{eff}}}} }}{{k_{f} w_{f} }}q_{fD,i} = 0$$

(53)

where

$$T_{i,i - 1} = \frac{1}{{1/\gamma_{i} + 1/\gamma_{i - 1} }} = \frac{{\gamma_{i} \gamma_{j} }}{{\gamma_{i} + \gamma_{j} }};$$

$$\gamma_{i} = \frac{{2\lambda_{g,i} }}{{\Delta l_{fD,i} }}.$$

(54)

For special blocks including grid blocks neighboring to intersections and blocks connecting to wells, expressions for transmissibilities can be derived following similar development in Zhang and Ayala (2018).

Appendix 4: System of Equations Solved in DRBEM-FVM

Given the free-gas DRBEM and FVM equations (Eqs. 27 and 32), they are rearranged to move any known information (e.g., known boundary specifications and results from previous time level) to the RHS of the equations. Taking the example of fluxes being specified at all outer boundary elements (q_B,sp), the unknowns to be solved are: outer boundary nodal concentrations (C_B), internal nodal concentration (C_I), fracture element source strengths (q_F) and concentrations (C_F). The resulting coefficient matrix and right-hand-side of system of equations are shown below:

(55)

where

$$\begin{gathered} {\mathbf{B}}_{{\mathbf{B}}}^{{{\mathbf{(n}}_{{\mathbf{t}}} {\mathbf{ + 1)}}}} {\mathbf{ = }}\frac{{{\mathbf{2X}}_{{{\mathbf{BB}}}} }}{{{\mathbf{\Delta t}}_{{\mathbf{D}}} }}{\mathbf{ + H}}_{{{\mathbf{BB}}}} , \, {\mathbf{B}}_{{\mathbf{I}}}^{{{\mathbf{(n}}_{{\mathbf{t}}} {\mathbf{ + 1)}}}} {\mathbf{ = }}\frac{{{\mathbf{2X}}_{{{\mathbf{BI}}}} }}{{{\mathbf{\Delta t}}_{{\mathbf{D}}} }}, \, {\mathbf{B}}_{{\mathbf{F}}}^{{{\mathbf{(n}}_{{\mathbf{t}}} {\mathbf{ + 1)}}}} {\mathbf{ = }}\frac{{{\mathbf{2X}}_{{{\mathbf{BF}}}} }}{{{\mathbf{\Delta t}}_{{\mathbf{D}}} }} \hfill \\ {\mathbf{I}}_{{\mathbf{B}}}^{{{\mathbf{(n}}_{{\mathbf{t}}} {\mathbf{ + 1)}}}} {\mathbf{ = }}\frac{{{\mathbf{2X}}_{{{\mathbf{IB}}}} }}{{{\mathbf{\Delta t}}_{{\mathbf{D}}} }}{\mathbf{ + H}}_{{{\mathbf{IB}}}} , \, {\mathbf{I}}_{{\mathbf{I}}}^{{{\mathbf{(n}}_{{\mathbf{t}}} {\mathbf{ + 1)}}}} {\mathbf{ = }}\frac{{{\mathbf{2X}}_{{{\mathbf{II}}}} }}{{{\mathbf{\Delta t}}_{{\mathbf{D}}} }},{\mathbf{I}}_{{\mathbf{F}}}^{{{\mathbf{(n}}_{{\mathbf{t}}} {\mathbf{ + 1)}}}} {\mathbf{ = }}\frac{{{\mathbf{2X}}_{{{\mathbf{IB}}}} }}{{{\mathbf{\Delta t}}_{{\mathbf{D}}} }}. \hfill \\ \end{gathered}$$

(56)

The unknowns to be solved in the system of equations given by Eq. 28 include N_b unknowns at outer boundary nodes (C_sD,n or q_sD,n, depending on the prevailing boundary condition of sorbed gas at boundaries), and L unknowns at internal nodes (C_sD,i). Similarly, if the fluxes are specified at all outer boundary elements (q_sB,sp), the unknowns to be solved would be: outer boundary nodal concentrations (C_sB), internal nodal concentration (C_sI) and fracture nodal concentrations (C_sF). The resulting coefficient matrix and right-hand side of system equations are shown below:

(57)

After all unknowns of bulk- and sorbed-phase equations (Eqs. 55 and 57) are solved for in terms of nodal values in boundary, fracture and internal domain, concentrations of any position inside domain can be straightforwardly calculated via the DRBEM equations where λ = 1. Iteration is needed to arrive the final solution of Eqs. 55 and 57 due to the dependency of sorption kinetic term E. A straightforward successive substitution scheme is implemented in this paper. Following Zhang and Ayala (2020), the integrals that are time-invariant and dependent on geometry only in Eq. 56 (B, I, F and X) are carried out analytically and derived from the expression of G shown in Eq. 15. Detailed analytical expression of these integrals and interpolation functions are provided in Zhang and Ayala (2020) and are thus not repeated here.