Addressing the Influence of a Heterogeneous Matrix on Well Performance in Fractured Rocks

Abstract

We address the influences of heterogeneity and complex geology of the matrix of fractured rocks on transient flow. Fractional constitutive flux laws reflect the stochastic framework that we consider. We model transient diffusion in both the matrix and fracture systems in terms of a continuous time random walk. Our procedure is particularly suited to address a complex geology that may exist on a number of scales and which may include dead ends and discontinuities. The performance of a horizontal well produced through arbitrarily located, multiple hydraulic fractures with distinct properties (length, width, permeability) forms the basis for our thesis. The pressure distribution in a rectangular drainage region where the well may be placed arbitrarily is expressed in terms of the Laplace transformation. The required solutions are obtained numerically. The focus of our work is on long-term behaviors for production at a constant pressure. In addition to numerical solutions, asymptotic results that provide information on the structure of the solutions are presented. Agreement between the asymptotic and numerical solutions is excellent. We show that long-term responses are governed by two distinct, two-parameter Mittag–Leffler functions and are an outcome of the complexities we desire to model in both the matrix and fracture systems. As a consequence, power-law behaviors that reflect the heterogeneity inherent in the system define long-time expectations. We show that the new solutions we derive do reduce to those of classical diffusion; that is, results corresponding to classical diffusion are a subset of the new results obtained here. Our results are particularly suited to model transient flow in shale reservoirs.

Appendices

Appendix 1: Model Formulation

We first briefly document details pertaining to the finite-conductivity model for a single fracture in a rectangular drainage region and then sketch out the details of the multiple-fracture model in that drainage region.

1.1 Single-Fracture Model

If the fractured well were to be located in a rectangle of sides $x_{eD} (x_\mathrm{e}/x_\mathrm{F})$ and $y_{eD} (y_\mathrm{e}/x_\mathrm{F})$ where the left corner is taken as the origin, then as summarized in Raghavan and Ozkan (1994) we have

$$\begin{aligned} \begin{aligned} \overline{p}_D (x_D, y_D)&= \frac{\pi }{ x_{eD} s^{1-\alpha _\mathrm{f}}}\left( {\frac{\tilde{\eta }}{\ell ^2}}\right) ^{\frac{1 -\alpha _\mathrm{f}}{\alpha _\mathrm{f}}}\int ^{x_{wD} + {x}_\mathrm{FD2}}_{x_{wD} - {x}_\mathrm{FD1}} \mathrm{d}x' \tilde{q}(x')\\&\quad \ \Biggl [{\frac{ch \sqrt{u} y_{eD1} +ch \sqrt{u} y_{eD2}}{\sqrt{u} \ sh \sqrt{u} y_{eD}}} + 2 \sum _{k=1}^{\infty } \cos \ k \pi \frac{x'}{ x_{eD}} \cos \ k \pi \frac{x_D}{ x_{eD}}\\&\quad \ {\frac{ch \sqrt{\epsilon }_k y_{eD1}+ ch \sqrt{\epsilon }_k y_{eD2}}{ {\sqrt{\epsilon }_k sh \sqrt{\epsilon }_k y_{eD}}}} \Biggr ], \end{aligned} \end{aligned}$$

(64)

with $\tilde{q}(x)$ being the flux distribution along the fracture face, $(x_w, y_w)$ being the well center, and ${x}_{F} = x_{F1} +x_{F2}$; see Fig. 11. All dimensionless distances are expressed with respect to the reference length, $\ell$: $y_{eD1} = y_{eD}-|y_D-y_{wD}|$, $y_{eD2} = y_{eD}-(y_D+y_{wD})$, ${x}_\mathrm{FD1} = {x}_\mathrm{FD2}$. The symbols ch(x) and sh(x) are the hyperbolic cosine function and hyperbolic sine function, respectively, and $\epsilon _k = {u+{{k^2\pi ^2}}/{{x^2_{eD}}}}$ with u obtained by the principle of subordination explained in Raghavan and Chen (2015)

$$\begin{aligned} u = \hat{s} f(\tilde{s}). \end{aligned}$$

(65)

To obtain the finite-conductivity solution, Eq. (64) is coupled with the formulation given in Cinco-Ley and Meng (1988) which leads to the following expression for the dimensionless pressure, $p_{wD}$, at the fractured well

$$\begin{aligned} \bar{p}_{wD} - \bar{p}_D(x_D, 0) = \frac{\pi }{C_\mathrm{FD}s} \left[ \vert x_D - x_{wD}\vert + s \int ^{\vert x_D - x_{wD}\vert }_0 \int ^{x'}_0 \bar{\tilde{q}}(x'')\mathrm{d}x''\mathrm{d}x'\right] , \end{aligned}$$

(66)

with the well center assumed to be at $(x_{wD}, 0)$; see Fig. 2, and $p_D(x_D, y_D, t)$ is the pressure distribution in the area drained by the well.

Fig. 11

A fractured well in a rectangular drainage region. The fracture with identical wing lengths with its center at $(x_w,y_w)$ is located anywhere in the drainage region and is parallel to one of the boundaries. The length, width and height of the fracture are $x_\mathrm{F}$, $w_\mathrm{F}$, and h, respectively. The horizontal well (not shown here) is assumed to intersect the fracture at $(x_w,y_w)$

If on the other hand, the fractured well were to be located in a reservoir that is infinite in its areal extent, then $\overline{p}_{D}$ is given by

$$\begin{aligned} \overline{p}_{D}={\frac{1}{ 2s^{1-\alpha _\mathrm{f}}}} \left( {\frac{\tilde{\eta }}{\ell ^2}}\right) ^{ \frac{ 1 -\alpha _\mathrm{f}}{\alpha _\mathrm{f}}} \int _{x_{wD}- {x}_\mathrm{FD1}}^{x_{wD}+ {x}_\mathrm{FD2}} \tilde{q}(x') K_0\left( \sqrt{\tilde{\epsilon }_{0} u}\right) \mathrm{d}x', \end{aligned}$$

(67)

where

$$\begin{aligned} \tilde{\epsilon _{0}} = \left( x_D -x_{wD}-x' \right) ^2+\left( y_D-y_{wD}\right) ^2. \end{aligned}$$

(68)

Additional details may be found in Chen and Raghavan (1997, 2013). Equation (64) forms the basis for coupling the hydraulic fractures under the assumption that fractures produce against well pressure, $p_\mathrm{wf}(t)$, and this matter is explored below.

1.2 The Multiple-Fracture Model

The scheme outlined in Raghavan et al. (1997) may be summarized as follows. Let the fractured well be at Location j and let $p_{D_{jk}}(x_{D}, y_{D})$ represent the pressure response caused by this well at Location k. If all the fractures produce against a common pressure, $p_{wD}$, then the pressure response within the horizontal wellbore is given by

$$\begin{aligned} p_{wD} = \sum ^{n}_{j=1} \int \limits ^{t_D}_0 q_{D_j}(\tau ) p'_{D_{jk}} (t_D-\tau )\mathrm{d}\tau , \end{aligned}$$

(69)

where $p'_{D_{jk}} (t_D) = d p_{D_{jk}}/ d t_D$. Equation (69) represents the convolution equation to compute the pressure response at the wellbore for the multiple-fracture system. Equation (69) along with Eqs. (43) and (44) in terms of the Laplace transformation leads to Eq. (45) in the text of this paper. Inversion of Eq. (45) by the Stehfest algorithm forms the basis for the results presented in the text.

To obtain an expression for production at a constant pressure, we use Duhamel’s theorem to combine the expression for $\overline{p}_{wD}$ with $\overline{q}_{D}$ through the relation

$$\begin{aligned} \overline{p}_{wD} \overline{q}_{D} = \frac{1}{s^2}. \end{aligned}$$

(70)

Appendix 2: Development of Asymptotic Solutions

We examine the structure of the solutions primarily during the boundary-dominated period for flow toward a fractured well in a rectangle and for flow to a well located at the center of a closed cylindrical reservoir where the focus is on deriving solutions for Flow Regime 4 and Flow Regime 5.

1.1 Well in a Rectangular Drainage Region

To obtain the needed expressions, we seek an expression for the long-time approximation of the pressure distribution given in Eq. (64) as $s \rightarrow 0$. It may be obtained by assuming that the flux distribution, $\tilde{q}(x)$, is a constant and is derived as follows. Integrating the right-hand side of Eq. (64) with respect to $x$ we obtain the following expression for the pressure distribution, $\overline{p}_{D}(x_D,y_D)$; see Ozkan & Raghavan (1991):

$$\begin{aligned} \begin{aligned} \overline{p}_D (x_D, y_D) =&\, H +\frac{2}{ s^{2-\alpha _\mathrm{f}}} \left( \frac{\tilde{\eta }}{\ell ^{2}}\right) ^{\frac{1-\alpha _\mathrm{f}}{\alpha _\mathrm{f}}}\sum _{k=1}^\infty \sin k \pi \frac{1}{ x_{eD}} \cos k \pi \frac{{x_{wD}}}{{ x_{eD}}}\\&\cos k \pi \frac{x_D}{x_{eD}} \frac{ch\, \epsilon _k (y_{eD}-\widetilde{y}_{D1}) + ch\, \epsilon _k (y_{eD}-\widetilde{y}_{D2})}{ k\epsilon _k \ sh\ \epsilon _k y_{eD}}, \end{aligned} \end{aligned}$$

(71)

with $\tilde{\epsilon }_k = \sqrt{{u+{{k^2\pi ^2}}/{{x^2_{eD}}}}}$, $\widetilde{y}_{D1}=|y_D-y_{wD}|$, and $\widetilde{y}_{D2}=y_D+y_{wD}$; and where

$$\begin{aligned} H = \frac{\pi }{ x_{eD}s^{2-\alpha _\mathrm{f}}}\left( \frac{\tilde{\eta }}{\ell ^{2}}\right) ^{\frac{1-\alpha _\mathrm{f}}{\alpha _\mathrm{f}}} \Biggl [ {\frac{ch \sqrt{u} (y_{eD}-\widetilde{y}_{D1}) + ch \sqrt{u} (y_{eD}-\widetilde{y}_{D2}) }{\sqrt{u} \ sh \sqrt{u} y_{eD}}}\Biggr ], \end{aligned}$$

(72)

with the reference length, $\ell$, being $x_\mathrm{F}/2 \equiv L_\mathrm{F}$.

We now seek the small-s approximation of the expression given in Eq. (71). The small-s approximation for $H$ may be obtained by first recasting Eq. (72) using the result (1.445,2) of Gradshteyn and Ryzhik (1965), namely,

$$\begin{aligned} \sum _{k=1}^\infty \frac{\cos k x}{k^2 +a^2} = \frac{\pi }{2a}\frac{\cosh [a (\pi -x)]}{\sinh (a \pi )} -\frac{1}{2 a^2 },\, \, \, [0 \le x \le 2\pi ]. \end{aligned}$$

(73)

On rewriting the right-hand side of Eq. (72), we obtain the small-s approximation of $H\, \text {as}\, s \rightarrow 0$

$$\begin{aligned} \begin{aligned} \lim _{s\rightarrow 0} H&=\left( \frac{\tilde{\eta }}{\ell ^{2}}\right) ^{\frac{1-\alpha _\mathrm{f}}{\alpha _\mathrm{f}}}\frac{1}{s^{2-\alpha _\mathrm{f}}}\left[ \frac{2\pi }{x_{eD}y_{eD}u}\right. \\&+\left. \frac{2\pi y_{eD}}{ x_{eD}} \left( \frac{1}{ 3} - \frac{\widetilde{y}_{D1}+\widetilde{y}_{D2}}{ {2y_{eD}}} + \frac{{\widetilde{y}_{D1}^2 + \widetilde{y}^2_{D2}}}{ {4y_{eD}^2}}\right) \right] , \end{aligned} \end{aligned}$$

(74)

where we have used the result (1.445,2) of Gradshteyn and Ryzhik (1965), namely,

$$\begin{aligned} \sum _{ k=1}^\infty \frac{\cos kx}{ k^2} = \frac{ \pi ^2}{6} - \frac{\pi x}{2 } +\frac{x^2}{ 4}, \, \, \, [0\le x \le 2 \pi ]. \end{aligned}$$

(75)

Based on the above development, at long enough times, we may write

$$\begin{aligned} \overline{p}_{D}(x_D, y_D, s \rightarrow 0) = \left( \frac{\tilde{\eta }}{\ell ^{2}}\right) ^{\frac{1-\alpha _\mathrm{f}}{\alpha _\mathrm{f}}}\frac{1}{s^{2-\alpha _\mathrm{f}}}\left( \frac{\tilde{A}}{u}+ \tilde{B}\right) , \end{aligned}$$

(76)

where

$$\begin{aligned} \tilde{A} = \frac{2\pi }{x_{eD}y_{eD}}, \end{aligned}$$

(77)

and

$$\begin{aligned} \begin{aligned} \tilde{B}&= 2\pi \frac{{y_{eD}}}{{x_{eD}}}\left( \frac{1}{ 3}-{\frac{\widetilde{y}_{D1}+\widetilde{y}_{D2}}{2y_{eD}}}+{\frac{\widetilde{y}^2_{D1}+\widetilde{y}^2_{D2}}{4y^2_{e D}}}\right. \\&\quad +\frac{{x^2_{eD}}}{\pi ^2 y_{eD}}\sum ^\infty _{k=1}{\frac{1}{k^2}}\sin k\pi {\frac{1}{x_{eD}}}\cos k \pi {\frac{x_{wD}}{x_{eD}}} \cos k \pi {\frac{x_D}{x_{eD}}} \\&\quad \left. {}{\frac{ch\, k\pi {\frac{(y_{eD}-\widetilde{y}_{D1})}{x_{eD}}}+ch\, k \pi {\frac{(y_{eD}-\widetilde{y}_{D2})}{x_{eD}}}}{sh\, k\pi {\frac{y_{eD}}{x_{eD}}}}}\right) . \end{aligned} \end{aligned}$$

(78)

Equation (76) may be used to derive the needed asymptotic expressions as detailed below.

1.1.1 Flow Regime 4

Using the expression for $f(\tilde{s})$ in Eq. (26) and simplifying Eq. (65), we obtain

$$\begin{aligned} u = \sqrt{\frac{\lambda ^\prime \omega ^\prime }{3}} \left( \frac{\ell ^2}{\tilde{\eta }}\right) ^{(2\alpha _\mathrm{f} - \alpha _\mathrm{m})/2}s^{(2\alpha _\mathrm{f} - \alpha _\mathrm{m})/2}, \end{aligned}$$

(79)

and on combining Eq. (76) with Eq. (79) and inverting, we have

$$\begin{aligned} p_{wD}( t_D \rightarrow \infty ) = \frac{2\sqrt{3}\pi }{x_{eD}y_{eD}\sqrt{\lambda ^\prime \omega ^\prime } } \frac{t_D^{\frac{2 - \alpha _\mathrm{m}}{2\alpha _\mathrm{f}}}}{\Gamma \left( \frac{(4 - \alpha _\mathrm{m})}{2}\right) } +\tilde{B}_w \frac{t_D^{\frac{1 - \alpha _\mathrm{f}}{\alpha _\mathrm{f}}}}{\Gamma (2 - \alpha _\mathrm{f})}; \end{aligned}$$

(80)

the subscript w on $\tilde{B}$ implies wellbore conditions; that is, $\tilde{B}_w = B(x_{wD} +x^{\star }_D, y_{wD})$, where $x^{\star }_D$ is defined in Table 10.5 of Raghavan (1993). Use of $x^{\star }_D$ enables us to obtain responses for $C_\mathrm{FD} \ge 5$. Actually, $\tilde{B}_w$ is related to the shape factor, $C_A$, by

$$\begin{aligned} \tilde{B}_w = \frac{1}{2}\ln \left( \frac{4 A}{e^\gamma C_A \ell ^2}\right) . \end{aligned}$$

(81)

Similarly, if we were to consider production at a constant pressure, then on combining Eqs. (70), (76) and (79), we have

$$\begin{aligned} \overline{q}_{D}( s \rightarrow 0) = \frac{1}{\tilde{B}_w s^{\alpha _\mathrm{f}} \left( \frac{\tilde{\eta }}{\ell ^{2}}\right) ^{\frac{1-\alpha _\mathrm{f}}{\alpha _\mathrm{f}}}}\frac{1}{\left[ 1 + \frac{\tilde{A}}{\tilde{B}_w}\sqrt{\frac{3}{\lambda ^\prime \omega ^\prime }} \left( \frac{\tilde{\eta }}{\ell ^{2}}\right) ^{\frac{2 \alpha _\mathrm{f} - \alpha _\mathrm{m}}{2}}s^{-\nu }\right] }; \, \,\, \nu = \alpha _\mathrm{f} -\alpha _\mathrm{m}/2,\nonumber \\ \end{aligned}$$

(82)

or

$$\begin{aligned} \overline{q}_{D}( s \rightarrow 0) = \frac{1}{\tilde{B}_w \left( \frac{\tilde{\eta }}{\ell ^{2}}\right) ^{\frac{1-\alpha _\mathrm{f}}{\alpha _\mathrm{f}}}}\sum _{ j=0}^\infty (-\tilde{C})^j s^{-\gamma }, \, \, \gamma = j \nu -\alpha _\mathrm{m}/2, \end{aligned}$$

(83)

where

$$\begin{aligned} \tilde{C} = \frac{\tilde{A}}{\tilde{B}_w}\sqrt{\frac{3}{\lambda ^\prime \omega ^\prime }}\left( \frac{\tilde{\eta }}{\ell ^{2}}\right) ^{\frac{2\alpha _\mathrm{f} - \alpha _\mathrm{m}}{2}}. \end{aligned}$$

(84)

On noting that

$$\begin{aligned} {{{\mathrm{\mathcal {L}}}}}\left[ \frac{t^{\gamma }}{\Gamma (\gamma +1)}\right] = \frac{1}{s^{\gamma +1}}; \, \, \, \gamma \ > \ -1, \end{aligned}$$

(85)

we may invert term by term as in Uchaikin (2013) leading to the result

$$\begin{aligned} q_D (t_D ) = \frac{1}{ \tilde{B}_w t_D^{\frac{1-\alpha _\mathrm{f}}{\alpha _\mathrm{f}}}} E_{{(2\alpha _\mathrm{f}- \alpha _\mathrm{m})/2},\alpha _\mathrm{f}}\left( - \frac{\tilde{A}}{\tilde{B}_w}\sqrt{\frac{3}{\lambda ^\prime \omega ^\prime }} t_D^{\frac{2\alpha _\mathrm{f} -\alpha _\mathrm{m}}{2\alpha _\mathrm{f}}}\right) , \end{aligned}$$

(86)

where $E_{\alpha , \beta }\left( z\right)$ is the two-parameter Mittag–Leffler function defined by Eq. (57).

1.1.2 Flow Regime 5

We now consider the approximation corresponding to Eq. (27). Using Eq. (27), the expression for u in Eq. (65) is

$$\begin{aligned} u = (1 + \omega ^\prime )\left( \frac{\ell ^2 s^{\alpha _\mathrm{f}}}{\tilde{\eta }_\mathrm{f}}\right) , \end{aligned}$$

(87)

and the pressure distribution at long enough times is

$$\begin{aligned} p_{D}(x_D,y_D, t_{D}\rightarrow \infty ) =\frac{2\pi }{x_{eD}y_{eD}(1 + \omega ')} t_{D}^{\frac{1}{\alpha _\mathrm{f}}} + \tilde{B}_w \frac{t_D^{\frac{1-\alpha _\mathrm{f}}{\alpha _\mathrm{f}}}}{\Gamma (2-\alpha _\mathrm{f})}. \end{aligned}$$

(88)

In the case of cylindrical systems, it has been shown in the literature that technically Eq. (88) only applies if Flow Regime 3 determines the start of boundary-dominated flow; that is, the influence of the external boundary manifests itself after the influence of the matrix system boundary becomes evident. If, however, Flow Regime 4 determines the time for the start of boundary-dominated flow, then it was determined numerically by Chen et al. (1985) that $\tilde{B}_w$ should be replaced by $\tilde{\tilde{B}}_w$, where

$$\begin{aligned} \tilde{\tilde{B}}_w = \tilde{B}_w + \frac{f(\omega ')}{\lambda ' x_{eD} y_{eD}}, \end{aligned}$$

(89)

where according to Ozkan et al. (1987) $f(\omega ')$ is given by

$$\begin{aligned} f(\omega ') = \frac{\omega '^2}{\lambda ' (1+ \omega ')^2}. \end{aligned}$$

(90)

The fact that Eq. (88) determines the well response during the boundary-dominated period appears to be independent of matrix properties and the ability of the matrix system to transfer fluids when transient conditions prevail in the matrix system at the start of boundary-dominated flow is, at first glance, disconcerting. Thus, the expression in Eq. (89), which was determined through numerical experiments, is indeed satisfying as it resolves intuitive expectations and, moreover, appears to be physically consistent. It may be shown that the extra pressure drop represented by the second term in Eq. (89) reflects the amount of fluid delivered to the fracture system per unit volume of the matrix system. We recommend that Eq. (89) be used for rectangular drainage regions, when appropriate.

We now examine constant pressure production. On using Eqs. (70), (76) and (89), we obtain

$$\begin{aligned} \overline{q}_D( s \rightarrow 0) = \frac{1}{\tilde{B}_w \left( \frac{\tilde{\eta }}{\ell ^{2}}\right) ^{\frac{1-\alpha _\mathrm{f}}{\alpha _\mathrm{f}}}}\frac{1}{\left[ s^\alpha + \frac{A}{\tilde{B}_w(1 + \omega ') }\left( \frac{\tilde{\eta }}{\ell ^{2}}\right) \right] }. \end{aligned}$$

(91)

Because

$$\begin{aligned} {{{\mathrm{\mathcal {L}}}}}^{-1}\left( \frac{1}{s^\nu +a}\right) = t^{\nu -1}E_{\nu , \nu }(- a t^\nu ), \end{aligned}$$

(92)

on inverting the expression in Eq. (91), we obtain

$$\begin{aligned} q_D (t_D \rightarrow \infty ) = \frac{1}{\tilde{B}_w t_D^{\frac{1-\alpha _\mathrm{f}}{\alpha _\mathrm{f}}}} E_{\alpha _\mathrm{f},\alpha _\mathrm{f}}\left[ - \frac{\tilde{A}}{\tilde{B}_w (1 + \omega ')} t_D\right] . \end{aligned}$$

(93)

This result is equivalent to that given in Raghavan and Chen (2016) for single-porosity systems. If Flow Regime 4 determines the beginning of the boundary-dominated period, then $\tilde{\tilde{B}}_w$ needs to be used instead of $\tilde{B}_w$ in the above expression.

1.2 Flow in a Cylindrical System

To show that the results are broadly applicable, we consider flow to a well of finite radius, $r_w$, that is located at the center of a cylindrical region and produced at a constant pressure, $p_w$, with its boundary, $r_\mathrm{e}$, sealed. Initially, the pressure is $p_i$ everywhere. We consider the analog of the transient diffusion equation given in Eq. (15) in cylindrical coordinates.

The desired solution of Eq. (15) in cylindrical coordinates needs to satisfy the following auxiliary conditions. Initially, we require that

$$\begin{aligned} p(r, t) = p_i, \, \, \text { when}\, t =0, \end{aligned}$$

(94)

and because the outer boundary is sealed, we require

$$\begin{aligned} \frac{\partial }{ \partial r} { p}(r,t) = 0, \, \, \text { for}\, r =r_\mathrm{e}. \end{aligned}$$

(95)

1.2.1 Production at a Constant Pressure

In this case we desire a solution of Eq. (15) (in cylindrical coordinates), Eqs. (94), (95) and

$$\begin{aligned} { p}(r,t) = p_w, \, \, \text {at}\, r =r_w, \end{aligned}$$

(96)

because the well is produced at a constant pressure. The production rate, q, for this system is given by

$$\begin{aligned} q = 2 \pi \tilde{\lambda }_\mathrm{f} h_\mathrm{ft} r_w \frac{\partial ^{1-\alpha _\mathrm{f}}}{\partial t^{1-\alpha _\mathrm{f}}}\left[ \frac{\partial }{ {\partial r}}p( {r},t) \right] , \, \, \text { for}\, r =r_w. \end{aligned}$$

(97)

We may readily show that the dimensionless flow rate, $q_D$, in terms of its Laplace transformation, $\overline{q} _D$, is given by

$$\begin{aligned} \overline{q}_D= & {} \left( \frac{r_w^2}{\tilde{\eta }_\mathrm{f}}\right) ^{\frac{1-\alpha _\mathrm{f}}{\alpha _\mathrm{f}}}\left[ \frac{{\sqrt{\hat{s} f(\tilde{s})}}}{s^{\alpha _{f}}}\right] \nonumber \\&\left[ \frac{I_1 (\sqrt{\hat{s} f(\tilde{s})} r_{eD}) K_1 (\sqrt{\hat{s} f(\tilde{s})}) - K_1 (\sqrt{\hat{s} f(\tilde{s})} r_{De} ) I_1 (\sqrt{\hat{s} f(\tilde{s})})}{K_1 (\sqrt{\hat{s} f(\tilde{s})} r_{De}) I_0 ({\sqrt{\hat{s} f(\tilde{s})}}) + I_1 (\sqrt{\hat{s} f(\tilde{s})} r_{De}) K_0 (\sqrt{\hat{s} f(\tilde{s})})}\right] . \end{aligned}$$

(98)

The scaling of the solution in terms of the dimensionless rate, $q_D$, naturally given in terms of the reference length, $r_w$, is given by

$$\begin{aligned} q_{D} (t_D) = \frac{\mu }{ 2 \pi \tilde{k}_{\alpha _{f}} h_\mathrm{ft} (p_i - p_w)}\left( {\frac{r_w^2}{\tilde{\eta }_\mathrm{f}}}\right) ^{\frac{1-\alpha _\mathrm{f}}{\alpha _\mathrm{f}}}q ( t). \end{aligned}$$

(99)

The long-term approximation of Eq. (98) for the flow rate may be obtained by using the ideas of Van Everdingen & Hurst (1949) and as outlined in Ozkan et al. (1987). The result is given by

$$\begin{aligned} \overline{q}_D (s \rightarrow 0) = \left( {\frac{r_w^2}{\tilde{\eta }}}\right) ^{\frac{1-\alpha _\mathrm{f}}{\alpha _\mathrm{f}}} \frac{\hat{s} f(\tilde{s})}{s^{\alpha _{f}}} \frac{{1}}{{\left( \ln r_{De} - \frac{3}{4} \right) }}\frac{1}{\left[ \frac{2/r_{De}^2}{\left( \ln r_{De} - \frac{3}{4} \right) } + \hat{s} f(\tilde{s})\right] }. \end{aligned}$$

(100)

In the following, we shall use two of the expressions for $f(\tilde{s})$ to illustrate characteristic behaviors during boundary-dominated flow, namely Flow Regimes 4 and 5.

1.2.2 Flow Regime 4

On combining Eqs. (100) and (26) and simplifying, we obtain

$$\begin{aligned} \overline{q}_D (s \rightarrow 0) = \left( {\frac{r_w^2}{\tilde{\eta }}}\right) ^{\frac{1-\alpha _\mathrm{f}}{\alpha _\mathrm{f}}} \frac{1}{s^{\alpha _\mathrm{f}}} \frac{{1}}{{\left( \ln r_{De} - \frac{3}{4} \right) }}\frac{1}{ \left( 1 + a b s^{- \nu } \right) }, \end{aligned}$$

(101)

with $\nu = \alpha _\mathrm{f} - \alpha _\mathrm{m}/2$,

$$\begin{aligned} a = \frac{2\sqrt{3}}{\sqrt{\lambda ^\prime \omega ^\prime } r_{De}^2 \left( \ln r_{De} - \frac{3}{4} \right) }, \end{aligned}$$

(102)

and

$$\begin{aligned} b = \sqrt{g(\tilde{\eta }_\mathrm{f})\frac{\tilde{\eta }_\mathrm{f}}{r_w^2}} = \left( \frac{\tilde{\eta }_\mathrm{f}}{r_w^2}\right) ^{\frac{2\alpha _\mathrm{f} -\alpha _\mathrm{m}}{2}}. \end{aligned}$$

(103)

On rewriting Eq. (101), as

$$\begin{aligned} \overline{q}_D (s \rightarrow 0 ) = \left( {\frac{r_w^2}{\tilde{\eta }}}\right) ^{\frac{1-\alpha _\mathrm{f}}{\alpha _\mathrm{f}}} \frac{{1}}{{\left( \ln r_{De} - \frac{3}{4} \right) }} \sum _{ j=0}^\infty (-ab)^j s^{-\gamma }, \, \, \gamma = j \nu +\alpha _\mathrm{f}, \end{aligned}$$

(104)

and on using Eq. (85), we may invert the right-hand side of Eq. (104) term by term to obtain

$$\begin{aligned} q_D (t_D \rightarrow \infty ) = \left( \frac{1}{\ln r_{eD} - \frac{3}{4} }\right) \frac{1}{t_D^{\frac{1-\alpha _\mathrm{f}}{\alpha _\mathrm{f}}}} E_{{(\alpha _\mathrm{f}- \alpha _\mathrm{m}/2)},\alpha _\mathrm{f}}\left( - a t_D^{\frac{2\alpha _\mathrm{f} -\alpha _\mathrm{m}}{2\alpha _\mathrm{f}}}\right) , \end{aligned}$$

(105)

where $E_{\alpha ,\beta } (z)$ is the two-parameter Mittag–Leffler function. If both $\alpha _\mathrm{f}$ and $\alpha _\mathrm{m}$ are set equal to 1, Eq. (105) leads to the expression given in Ozkan et al. (1987) because $E_{1/2,1}(-\sqrt{x}) = \exp (x) \mathrm {erfc}\,(x)$.

1.2.3 Flow regime 5

We now consider the approximation corresponding to Eq. (27). On combining Eqs. (100) and (27) and simplifying, we obtain

$$\begin{aligned} a = \frac{2}{(1 +\omega ^\prime ) r_{De}^2 \left( \ln r_{De} - \frac{3}{4} \right) }, \end{aligned}$$

(106)

and

$$\begin{aligned} b = 1. \end{aligned}$$

(107)

Thus, in this case

$$\begin{aligned} q_D (t_D \rightarrow \infty ) = \left( \frac{1}{\ln r_{De} - \frac{3}{4} }\right) \frac{1}{t_D^{\frac{1-\alpha _\mathrm{f}}{\alpha _\mathrm{f}}}} E_{\alpha _\mathrm{f},\alpha _\mathrm{f}}\left[ - \frac{2 t_D}{(1 +\omega ^\prime ) r^2_{De}\left( {\ln r_{De} - \frac{3}{4} }\right) }\right] . \end{aligned}$$

(108)

In concluding this discussion we need to recall that if Flow Regime 4 determines the beginning of the boundary-dominated period, then $\tilde{\tilde{B}}_w$ needs to be used instead of $\tilde{B}_w$ in the above expression.

Appendix 3: Steady Flow, Boundaries at Constant (Initial) Pressure

We consider subdiffusive flow under steady conditions with a goal to simply show the general characteristics of responses under such conditions and consider accordingly a fractured well in a rectangular area where all four boundaries of the porous medium are at a constant (initial) pressure. Again, using the principle of subordination and the expression for the pressure distribution given in Raghavan and Ozkan (1994) for a similar system under classical diffusion, we obtain

$$\begin{aligned} \begin{aligned} \overline{p}_D(x_D,y_D)=&\frac{2}{s^{2 -\alpha }}\sum ^\infty _{k=1} \frac{1}{ k}\sin k\pi {{x_D}\over {x_{eD}}}\sin k\pi \frac{{x_{wD}}}{{x_{eD}}} \sin k\pi \frac{1}{{x_{eD}}}\\&{\frac{ch\, \epsilon _k (y_{eD}-\widetilde{y}_{D1}) - ch \,\epsilon _k (y_{eD}-\widetilde{y}_{D2}) }{{\epsilon _k \ sh\, \epsilon _k y_{eD}}}}. \end{aligned} \end{aligned}$$

(109)

All expressions in Eq. (109) have been defined previously. We have dropped the subscript on the exponent $\alpha$ as we are considering a homogenous system. At long enough times, if we consider the expression in Eq. (109) as $s \rightarrow 0$ ($u \rightarrow 0$), then the inversion of the expression on the right-hand side of Eq. (109) is a simple matter, and we obtain

$$\begin{aligned} \begin{aligned} p_D(x_D,y_D, t_D \rightarrow \infty )=&\frac{2x_{eD}}{\pi }\sum ^\infty _{k=1} \frac{1}{ k^2}\sin k\pi {{x_D}\over {x_{eD}}}\sin k\pi \frac{{x_{wD}}}{{x_{eD}}} \sin k\pi \frac{1}{{x_{eD}}}\\&{\frac{ch\, k \pi \frac{ (y_{eD}-\widetilde{y}_{D1})}{x_{eD}} - ch \, k \pi \frac{ (y_{eD}-\widetilde{y}_{D2})}{x_{eD}}}{{ \ sh\, k \pi \frac{ y_{eD}}{x_{eD}}}}}\frac{t_D^{\frac{1 -\alpha }{\alpha }}}{\Gamma (2 -\alpha )}. \end{aligned} \end{aligned}$$

(110)

Similarly, for production of the well at a constant pressure we obtain

$$\begin{aligned} q_D( t_D \rightarrow \infty )=\frac{\pi }{2 \hat{B}_w x_{eD}}\frac{1}{\Gamma (\alpha ) t_D^{\frac{1 -\alpha }{\alpha }}}. \end{aligned}$$

(111)

Here $\hat{B}_w$ is the value of the summation term in Eq. (110) corresponding to the wellbore; that is,

$$\begin{aligned} \begin{aligned} \hat{B}_w =&\sum ^\infty _{k=1} \frac{1}{k^2}\sin k\pi {{x_D}\over {x_{eD}}}\sin k\pi \frac{{x_{wD}}}{{x_{eD}}} \sin k\pi \frac{1}{{x_{eD}}} \\&{\frac{ch\, k \pi \frac{ (y_{eD}-\widetilde{y}_{D1})}{x_{eD}} - ch \, k \pi \frac{ (y_{eD}-\widetilde{y}_{D2})}{x_{eD}}}{{ \ sh\, k \pi \frac{ y_{eD}}{x_{eD}}}}}. \end{aligned} \end{aligned}$$

(112)

with $x_D = x_{wD} +x^{\star }_D, y_D = y_{wD}$. Again this expression may be used for $C_\mathrm{FD} \ge 5$ as noted for $\tilde{B}_w$. It is clear that both Eqs. (110) and (111) suggest counterintuitive consequences under steady flow for subdiffusive conditions.