A Fractional Decline Model Accounting for Complete Sequence of Regimes for Production from Fractured Unconventional Reservoirs

Abstract

On the basis of a well-based model (Model I) developed in a previous work (Liu and Valkó in SPE J 2019. https://doi.org/10.2118/197049-PA), in which a fractional production decline model is developed based on anomalous diffusion to structurally account for the heterogeneity related to the complex fracture network, we incorporate two more components, i.e., the tempered anomalous diffusion and a source term, in the present work to develop a generalized model, namely Model II. Because of the two new components, Model II is capable of physically considering both the scale lower bound of heterogeneity and the influx from the matrix into the conductive fracture system. It consequently enables the new model to describe the complete sequence of regimes for the production of the slightly compressible single-phase fluid from the fractured unconventional reservoirs. Then, in the case studies, the synthetic data used in the previous work are better fitted to the new type curves in all observed periods, which demonstrates the advantage of Model II over Model I regarding the late-time flow regimes. In addition, the good fit of other sets of synthetic data to the Model II type curves exhibits the claimed capabilities of the new model and indicates the production of this type could be well characterized by seven parameters, i.e., $\alpha$, $c$, ${\lambda }_{\mathrm{D}}$, $\omega$, $\sigma$, $\tau$, and $\mathrm{EUR}$. Finally, the characteristics of Model II are further investigated by analyzing the effects of several parameters on the production performance, based on which some insights about the practice of developing unconventional reservoirs using massive hydraulic fracturing are provided.

Appendices

Appendix 1: Derivation of the Hemispherical Steady State Flow

The extra pressure drop in the convergence region is modeled by the stabilized hemispherical Darcy’s flow of radius ${r}_{{\rm c}}$. The hemisphere is chopped by two parallel cross sections with radius of ${r}_{{\rm c}}$ and ${r}_{{\rm p}}$(Fig. 3c). The coordinates shown in Fig. 3d are utilized to solve the problem.

In the coordinates, we assume constant pressure condition $P={P}_{{\rm e}}$ at $X=0$ and $P={P}_{{\rm f}}$ at $X=\sqrt{{r}_{{\rm c}}^{2}-{r}_{{\rm p}}^{2}}$. Given the condition of stabilized flow, the flow rate across an arbitrary intermediate cross section of $X=x$ also keeps a constant of ${q}_{{\rm p}}$, as shown in Eq. (34):

$${q}_{{\rm p}}=-\frac{{k}_{{\rm f}}}{\mu }\frac{\mathrm{d}P}{\mathrm{d}x} \pi \left({r}_{{\rm c}}^{2}-{x}^{2}\right).$$

(34)

After separating the variables in Eq. (34), Eq. (35) is obtained

$$-{\rm d}P=\frac{{q}_{{\rm p}}\mu }{ \pi {k}_{{\rm f}}}\frac{\mathrm{d}x}{\left({r}_{{\rm c}}^{2}-{x}^{2}\right)}.$$

(35)

Integrating both sides of Eq. (35) over the whole range, we can get Eq. (36):

$$-{\int }_{{P}_{{\rm e}}}^{{P}_{{\rm f}}}\mathrm{d}P=\frac{{q}_{{\rm p}}\mu }{ \pi {k}_{{\rm f}}}{\int }_{0}^{\sqrt{{r}_{{\rm c}}^{2}-{r}_{{\rm p}}^{2}}}\frac{\mathrm{d}x}{\left({r}_{{\rm c}}^{2}-{x}^{2}\right)}.$$

(36)

Finally, the relationship between the extra pressure drop and the flow rate can be obtained as Eq. (37):

$$\Delta {P}_{{\rm c}}={P}_{{\rm e}}-{P}_{{\rm f}}={q}_{{\rm p}}\frac{\mu }{\pi {k}_{{\rm f}}{r}_{{\rm c}}}\mathrm{arctanh}\left(\sqrt{1-\frac{{r}_{{\rm p}}^{2}}{{r}_{{\rm c}}^{2}}}\right).$$

(37)

Appendix 2: Solution of the Boundary Value Problem in Laplace Space for the System of Model II

In this appendix, the solution of the system of Model II is obtained by solving the boundary value problem of Eq. (20).

Regarding the equations of the matrix continuum, the general solution of Eq. (20d) is

$$\widetilde{{P}_{{\rm mD}}}={C}_{1}\mathrm{exp}\left(\sqrt{\frac{s}{\sigma }}{\xi }_{{\rm D}}\right)+{C}_{2}\mathrm{exp}\left(-\sqrt{\frac{s}{\sigma }}{\xi }_{{\rm D}}\right)$$

(38)

where ${C}_{1}$ and ${C}_{2}$ are the constants of integration with arbitrary values.

Taking the derivative of Eq. (38) with respect to ${\xi }_{{\rm D}}$, we can get

$$\frac{\partial \widetilde{{P}_{{\rm mD}}}}{\partial {\xi }_{{\rm D}}}={C}_{1}\sqrt{\frac{s}{\sigma }}\mathrm{exp}\left(\sqrt{\frac{s}{\sigma }}{\xi }_{{\rm D}}\right)-{C}_{2}\sqrt{\frac{s}{\sigma }}\mathrm{exp}\left(-\sqrt{\frac{s}{\sigma }}{\xi }_{{\rm D}}\right).$$

(39)

Substituting Eqs. (38) and (39) into the boundary conditions of the matrix continuum, i.e., Eqs. (20e) and (20f), the equations with respect to ${C}_{1}$ and ${C}_{2}$ are obtained, as shown in Eq. (40):

$$\left\{\begin{array}{l}\mathrm{exp}\left(\sqrt{\frac{s}{\sigma }}\right){C}_{1}-\mathrm{exp}\left(-\sqrt{\frac{s}{\sigma }}\right){C}_{2}=0\\ {C}_{1}+{C}_{2}=\frac{1}{s}\end{array}\right.$$

(40)

${C}_{1}$and ${C}_{2}$ are obtained as Eq. (41) from Eq. (40):

$$\left\{\begin{array}{l}{C}_{1}=\frac{1}{s+ s\mathrm{exp}\left(2 \sqrt{\frac{s}{\sigma }}\right)}\\ {C}_{2}=\frac{\mathrm{exp}\left(2 \sqrt{\frac{s}{\sigma }}\right)}{s+ s\mathrm{exp}\left(2 \sqrt{\frac{s}{\sigma }}\right)}\end{array}\right.$$

(41)

Substituting Eq. (41) into Eq. (39), we can get

$$\frac{\partial \widetilde{{P}_{{\rm mD}}}}{\partial {\xi }_{{\rm D}}}=\frac{\mathrm{exp}\left(\sqrt{\frac{s}{\sigma }}{\xi }_{{\rm D}}\right)-\mathrm{exp}\left(\sqrt{\frac{s}{\sigma }}\left(2-{\xi }_{{\rm D}}\right) \right)}{\sqrt{s\sigma }+ \sqrt{s\sigma }\mathrm{exp}\left(2 \sqrt{\frac{s}{\sigma }}\right)}$$

(42)

Thus, the gradient at the interface is obtained by evaluating Eq. (42) at ${\xi }_{{\rm D}}=0$:

$${\left(\frac{\partial \widetilde{{P}_{{\rm mD}}}}{\partial {\xi }_{{\rm D}}}\right)}_{{\xi }_{{\rm D}}=0}=\frac{1-\mathrm{exp}\left(2\sqrt{\frac{s}{\sigma }} \right)}{\sqrt{s\sigma }+ \sqrt{s\sigma }\mathrm{exp}\left(2 \sqrt{\frac{s}{\sigma }}\right)}=-\frac{\mathrm{tanh}\left(\sqrt{\frac{s}{\sigma }} \right)}{\sqrt{s\sigma }}.$$

(43)

Then, substituting the expression in Eq. (43) into Eq. (20a) and rearranging the terms yield

$$\left(\omega +\left(1-\omega \right)\mathrm{ tanh}\left(\sqrt{s/\sigma }\right)/\sqrt{s/\sigma }\right)\left({\left(s+{\lambda }_{{\rm D}}\right)}^{\alpha }-{{\lambda }_{{\rm D}}}^{\alpha }\right)\widetilde{{P}_{{\rm fD}}}=\frac{{\partial }^{2}}{\partial {y}_{{\rm D}}^{2}}\widetilde{{P}_{{\rm fD}}}$$

(44)

Let $f\left(\alpha ,{\lambda }_{{\rm D}},s\right)=\sqrt{{\left(s+{\lambda }_{{\rm D}}\right)}^{\alpha }-{{\lambda }_{{\rm D}}}^{\alpha }}$ and $g\left(\omega , \sigma ,s\right)=\sqrt{\omega +\left(1-\omega \right)\mathrm{ tanh}\left(\sqrt{s/\sigma }\right)/\sqrt{s/\sigma }}$. Therefore, Eq. (44) can be written as Eq. (45):

$${\left(f\left(\alpha ,{\lambda }_{{\rm D}},s\right)g\left(\omega , \sigma ,s\right)\right)}^{2}\widetilde{{P}_{{\rm fD}}}=\frac{{\partial }^{2}}{\partial {y}_{{\rm D}}^{2}}\widetilde{{P}_{{\rm fD}}}$$

(45)

With the same form as Eqs. (20d), (45) has the general solution as shown in Eq. (46):

$$\widetilde{{P}_{{\rm fD}}}={C}_{3}\mathrm{exp}\left(fg{y}_{{\rm D}}\right)+{C}_{4}\mathrm{exp}\left(-fg{y}_{{\rm D}}\right)$$

(46)

where ${C}_{3}$ and ${C}_{4}$ are the another set of constants of integration with arbitrary.

Taking the derivative of Eq. (46) with respect to ${y}_{{\rm D}}$, we can get

$$\frac{\partial \widetilde{{P}_{{\rm fD}}}}{\partial {y}_{{\rm D}}}={C}_{3}fg\mathrm{exp}\left(fg{y}_{{\rm D}}\right)-{C}_{4}fg\mathrm{exp}\left(-fg{y}_{{\rm D}}\right)$$

(47)

Substituting Eqs. (46) and (47) into the boundary conditions of the fracture continuum, i.e., Eqs. (20b) and (20c), a system of equations with respect to ${C}_{3}$ and ${C}_{4}$ is obtained, as shown in Eq. (48):

$$\left\{\begin{array}{l}\mathrm{exp}\left(fg\right){C}_{3}-\mathrm{exp}\left(-fg\right){C}_{4}=0\\ \left(1-c\,sg/f\right){C}_{3}+\left(1+c\,sg/f\right){C}_{4}=\frac{1}{s}\end{array}\right.$$

(48)

By solving Eq. (48), ${C}_{3}$ and ${C}_{4}$ are expressed as Eq. (49)

$$\left\{\begin{array}{l}{C}_{3}=\frac{1}{s\left(1+\mathrm{exp}\left(2fg\right)\right)+c {s}^{2}g/f\left(\mathrm{exp}\left(2fg\right)-1\right)}\\ {C}_{4}=\frac{\mathrm{exp}\left(2 \sqrt{fg}\right)}{s\left(1+\mathrm{exp}\left(2fg\right)\right)+c {s}^{2}g/f\left(\mathrm{exp}\left(2fg\right)-1\right)}\end{array}\right.$$

(49)

Substituting Eq. (49) into Eq. (46), we can get

$$\widetilde{{P}_{{\rm fD}}}=\frac{\mathrm{exp}\left(fg{y}_{{\rm D}}\right)+\mathrm{exp}\left(2 fg\right)\mathrm{exp}\left(-fg{y}_{{\rm D}}\right)}{s\left(1+\mathrm{exp}\left(2fg\right)\right)+c {s}^{2}g/f\left(\mathrm{exp}\left(2fg\right)-1\right)}$$

(50)

After doing some further simplifications to Eq. (50), finally we can obtain $\widetilde{{P}_{{\rm fD}}}$ [Eq. (51)]:

$$\widetilde{{P}_{{\rm fD}}}\left({y}_{{\rm D}}, s\right)=\frac{\mathrm{cosh}\left(\left({y}_{{\rm D}}-1\right) fg\right)}{s \mathrm{cosh}\left(fg\right)+c {s}^{2}\mathrm{sinh}\left(fg\right)g/f}.$$

(51)

Appendix 3: Derivation of the Relationship Between Rate and Pressure

In this appendix, the relationship between the production rate and the pressure in the setting of the tempered fractional diffusivity equation is derived to fill the gap between Eqs. (21) and (23) as well as (24).

Within the dimensional domain, the production rate is expressed as Eq. (52):

$$q=2{h}_{{\rm f}}{N}_{{\rm p}}{x}_{{\rm c}}{\left|{v}_{\alpha ,\lambda }\right|}_{y=0}$$

(52)

The factor 2 stands for the wellbore receiving production from both perforations for a typical perforation pair.

Performing Laplace transform on both sides of Eq. (52), we can get

$$\widetilde{q}=2{h}_{{\rm f}}{N}_{{\rm p}}{x}_{{\rm c}}{\left|-\frac{{k}_{{\rm f}}^{*}}{\mu }\frac{s^{\prime}}{{\left(s^{\prime}+\lambda \right)}^{\alpha }-{\lambda }^{\alpha }}\widetilde{\frac{\partial {P}_{{\rm f}}}{\partial y}}\right|}_{y=0}$$

(53)

according to the definition of the tempered fractional flux law in Eq. (10). $s^{\prime}$ is the Laplace variable with respect to the dimensional time. Given the fact that Laplace transform is only performed with respect to the temporal coordinate, we have $\widetilde{\frac{\partial {P}_{{\rm f}}}{\partial y}}=\frac{\partial \widetilde{{P}_{{\rm f}}}}{\partial y}$. Therefore,

$$\widetilde{q}=2{h}_{{\rm f}}{N}_{{\rm p}}{x}_{{\rm c}}\frac{{k}_{{\rm f}}^{*}}{\mu }\frac{s^{\prime}}{{\left(s^{\prime}+\lambda \right)}^{\alpha }-{\lambda }^{\alpha }}{\left(\frac{\partial \widetilde{{P}_{{\rm f}}}}{\partial y}\right)}_{y=0}.$$

(54)

Due to the zero Neumann boundary condition on the other boundary of the problem domain [Eq. (13c)], Eq. (54) could be rewritten as Eq. (55):

$$\widetilde{q}=-2{h}_{{\rm f}}{N}_{{\rm p}}{x}_{{\rm c}}\frac{{k}_{{\rm f}}^{*}}{\mu }\frac{s^{\prime}}{{\left(s^{\prime}+\lambda \right)}^{\alpha }-{\lambda }^{\alpha }}\left[{\left(\frac{\partial \widetilde{{P}_{{\rm f}}}}{\partial y}\right)}_{y={y}_{{\rm e}}}-{\left(\frac{\partial \widetilde{{P}_{{\rm f}}}}{\partial y}\right)}_{y=0}\right]$$

(55)

Applying the fundamental theorem of calculus to Eq. (55), we can get

$$\widetilde{q}=-2{h}_{{\rm f}}{N}_{{\rm p}}{x}_{{\rm c}}\frac{{k}_{{\rm f}}^{*}}{\mu }\frac{s^{\prime}}{{\left(s^{\prime}+\lambda \right)}^{\alpha }-{\lambda }^{\alpha }}{\int }_{0}^{{y}_{{\rm e}}}\frac{\partial }{\partial y}\frac{\partial \widetilde{{P}_{{\rm f}}}}{\partial y}\mathrm{d}y.$$

(56)

After substituting the dimensionless groups into Eq. (56), the corresponding dimensionless form is obtained;

$$\widetilde{{q}_{{\rm D}}}=\frac{s}{{\left(s+{\lambda }_{{\rm D}}\right)}^{\alpha }-{\lambda }_{{\rm D}}^{\alpha }}{\int }_{0}^{1}\frac{{\partial }^{2}\widetilde{{P}_{{\rm fD}}}}{\partial {y}_{{\rm D}}^{2}}{\rm d}{y}_{{\rm D}}$$

(57)

where $s$ is the Laplace variable associated with the dimensionless time.

Then, substituting Eq. (44) into Eq. (57), we can get

$$\widetilde{{q}_{{\rm D}}}=s\left(\omega +\left(1-\omega \right)\mathrm{ tanh}\left(\sqrt{s/\sigma }\right)/\sqrt{s/\sigma }\right){\int }_{0}^{1}\widetilde{{P}_{{\rm fD}}}{\rm d}{y}_{{\rm D}}$$

(58)

which is the relationship we need for obtaining Eq. (24) from Eq. (21).

Additionally, with the relationship between production rate and cumulative production, we have

$$Q={\int }_{0}^{T}q\mathrm{d}t$$

(59)

Hence, in Laplace space Eq. (59) would be written as Eq. (60):

$$\widetilde{Q}=\frac{\widetilde{q}}{s^{\prime}}$$

(60)

Then, in dimensionless form, we have

$$\widetilde{{Q}_{{\rm D}}}=\frac{\widetilde{{q}_{{\rm D}}}}{s}=\left(\omega +\left(1-\omega \right)\mathrm{ tanh}\left(\sqrt{s/\sigma }\right)/\sqrt{s/\sigma }\right){\int }_{0}^{1}\widetilde{{P}_{{\rm fD}}}{\rm d}{y}_{{\rm D}}$$

(61)

after substituting Eq. (58) for $\widetilde{{q}_{{\rm D}}}$. Equation (61) is the expression we need for obtaining Eq. (23) from Eq. (21).

Appendix 4: Equivalence Between the Proposed and the Reference Models

According to Eq. (32), the proposed and the reference models have different time scales, which would lead to some factors between their respective Laplace variables as well as the functions in Laplace space. At first, these factors will be derived. For simplicity, we denote the group ${y}_{{\rm e}}^{2}/{A}_{{\rm cw}}$ in Eq. (32) as $C$.

Regarding a function $F$ with respect to ${t}_{{\rm D}}$, its Laplace transform could be performed by

$$\widetilde{F}\left(s\right)={\int }_{0}^{\infty }F\left({t}_{{\rm D}}\right){e}^{-s{t}_{{\rm D}}}{\rm d}{t}_{{\rm D}}=\frac{1}{C}{\int }_{0}^{\infty }F\left({t}_{{\rm D}}\left({t}_{{\rm D}}^{*}\right)\right){e}^{-\frac{s}{C}{t}_{{\rm D}}^{*}}\mathrm{d}{t}_{{\rm D}}^{*}$$

(62)

If we denote the Laplace variable with respect to ${t}_{{\rm D}}^{*}$ as ${s}^{*}$, then from Eq. (62) we get

$${s}^{*}=\frac{s}{C}$$

(63)

$$\widetilde{F}\left(s\right)=\frac{\widetilde{F}\left({s}^{*}\right)}{C}.$$

(64)

Next substituting $\left(\alpha , {\lambda }_{{\rm D}}, c\right)=\left(1.0, 0.0, 0.0\right)$ into Eq. (24), we could get

$$\frac{1}{\widetilde{{q}_{{\rm D}}}\left(s\right)}=\frac{{s}^{1/2}}{\sqrt{\omega +\left(1-\omega \right)\mathrm{ tanh}\left(\sqrt{s/\sigma }\right)/\sqrt{s/\sigma }}} \mathrm{coth}\left({s}^{1/2} \sqrt{\omega +\left(1-\omega \right)\mathrm{ tanh}\left(\sqrt{s/\sigma }\right)/\sqrt{s/\sigma }}\right)$$

(65)

In the current scenario, Eq. (18 g) has the form of Eq. (66):

$$\sigma =\frac{{\left({y}_{{\rm e}}^{2}/{\eta }_{{\rm f}}^{*}\right)}^{\frac{1}{\alpha }}}{{L}_{{\rm m}}^{2}/{\eta }_{{\rm m}}}=\frac{{y}_{{\rm e}}^{2}}{{L}_{{\rm m}}^{2}}\frac{{k}_{{\rm m}}/\left(\mu {\left(\phi {c}_{{\rm t}}\right)}_{{\rm m}}\right)}{{k}_{{\rm f}}/\left(\mu {\left(\phi {c}_{{\rm t}}\right)}_{f+m}\right)}=\frac{{y}_{{\rm e}}^{2}}{{L}_{{\rm m}}^{2}}\frac{{k}_{{\rm m}}}{{k}_{{\rm f}}}\frac{1}{1-\omega }=\frac{4}{{L}^{2}}\frac{{k}_{{\rm m}}}{{k}_{{\rm f}}}\frac{{y}_{{\rm e}}^{2}}{1-\omega }$$

(66)

The definition of the interporosity flow coefficient ${\lambda }_{{\rm Ac}}$ in the reference model (Bello and Wattenbarger 2010) is

$${\lambda }_{{\rm Ac}}=\frac{12}{{L}^{2}}\frac{{k}_{{\rm m}}}{{k}_{{\rm f}}}{A}_{{\rm cw}}$$

(67)

Combining Eq. (66) and (67), the relationship between the classic interporosity flow coefficient and the proposed characteristic time ratio is

$$\sigma =\frac{{y}_{{\rm e}}^{2}}{3\left(1-\omega \right){A}_{{\rm cw}}}{\lambda }_{{\rm Ac}}.$$

(68)

By Eq. (68) and (63), we know that

$$\sqrt{\frac{s}{\sigma }}=\sqrt{\frac{3\left(1-\omega \right){A}_{{\rm cw}}}{{y}_{{\rm e}}^{2}{\lambda }_{{\rm Ac}}}C{s}^{*}}=\sqrt{\frac{3\left(1-\omega \right){A}_{{\rm cw}}}{{y}_{{\rm e}}^{2}{\lambda }_{{\rm Ac}}}\frac{{y}_{{\rm e}}^{2}}{{A}_{{\rm cw}}}{s}^{*}}=\sqrt{\frac{3\left(1-\omega \right){s}^{*}}{{\lambda }_{{\rm Ac}}}}$$

(69)

After substituting Eq. (69), the argument of function $\mathrm{coth}$ in Eq. (65) becomes

$${s}^{1/2} \sqrt{\omega +\left(1-\omega \right)\mathrm{ tanh}\left(\sqrt{s/\sigma }\right)/\sqrt{s/\sigma }}=\frac{{y}_{{\rm e}}}{\sqrt{{A}_{{\rm cw}}}}{{s}^{*}}^{1/2} \sqrt{\omega +\frac{\left(1-\omega \right)}{\sqrt{\frac{3\left(1-\omega \right){s}^{*}}{{\lambda }_{{\rm Ac}}}}}\mathrm{ tanh}\sqrt{\frac{3\left(1-\omega \right){s}^{*}}{{\lambda }_{{\rm Ac}}}}}=\frac{{y}_{{\rm e}}}{\sqrt{{A}_{{\rm cw}}}}{{s}^{*}}^{1/2} \sqrt{\omega +\sqrt{\frac{{\lambda }_{{\rm Ac}}\left(1-\omega \right)}{3{s}^{*}}}\mathrm{ tanh}\sqrt{\frac{3\left(1-\omega \right){s}^{*}}{{\lambda }_{{\rm Ac}}}}}$$

(70)

According to Eq. (2) in Bello and Wattenbarger (2010), Eq. (70) is directly related to the function $f$ of the reference model as in Eq. (71):

$${s}^{1/2} \sqrt{\omega +\left(1-\omega \right)\mathrm{ tanh}\left(\sqrt{s/\sigma }\right)/\sqrt{s/\sigma }}=\sqrt{{s}^{*} f\left({s}^{*}\right)}\frac{{y}_{{\rm e}}}{\sqrt{{A}_{{\rm cw}}}}.$$

(71)

Thus, Eq. (65) becomes

$$\frac{1}{\widetilde{{q}_{{\rm D}}}\left(s\right)}=\frac{s}{{s}^{1/2}\sqrt{\omega +\left(1-\omega \right)\mathrm{ tanh}\left(\sqrt{s/\sigma }\right)/\sqrt{s/\sigma }}} \mathrm{coth}\left({s}^{1/2} \sqrt{\omega +\left(1-\omega \right)\mathrm{ tanh}\left(\sqrt{s/\sigma }\right)/\sqrt{s/\sigma }}\right)$$

$$=\frac{{y}_{{\rm e}}^{2}/{A}_{{\rm cw}}{s}^{*}}{\sqrt{{s}^{*} f\left({s}^{*}\right)}{y}_{{\rm e}}/\sqrt{{A}_{{\rm cw}}}} \mathrm{coth}\left(\sqrt{{s}^{*} f\left({s}^{*}\right)}\frac{{y}_{{\rm e}}}{\sqrt{{A}_{{\rm cw}}}}\right)$$

$$=\frac{{y}_{{\rm e}}}{2\pi \sqrt{{A}_{{\rm cw}}}}\frac{2\pi {s}^{*}}{\sqrt{{s}^{*} f\left({s}^{*}\right)}} \mathrm{coth}\left(\sqrt{{s}^{*} f\left({s}^{*}\right)}\frac{{y}_{{\rm e}}}{\sqrt{{A}_{{\rm cw}}}}\right)$$

(72)

Then, we could observe that the right-hand side of Eq. (72) includes the rate solution of the reference model [Eq. (2A) in Bello and Wattenbarger (2010)]. Hence, we obtain

$$\frac{1}{\widetilde{{q}_{{\rm D}}}\left(s\right)}=\frac{{y}_{{\rm e}}}{2\pi \sqrt{{A}_{{\rm cw}}}}\frac{1}{\widetilde{{q}_{{\rm D}}^{*}}\left({s}^{*}\right)}.$$

(73)

According to Eq. (64), we could use a single Laplace variable on both sides in Eq. (73):

$$\frac{1}{\widetilde{{q}_{{\rm D}}}\left(s\right)}=\frac{{y}_{{\rm e}}}{2\pi \sqrt{{A}_{{\rm cw}}}}\frac{1}{{y}_{{\rm e}}^{2}/{A}_{{\rm cw}}\widetilde{{q}_{{\rm D}}^{*}}\left(s\right)}=\frac{\sqrt{{A}_{{\rm cw}}}}{2\pi {y}_{{\rm e}}}\frac{1}{\widetilde{{q}_{{\rm D}}^{*}}\left(s\right)}.$$

(74)

Finally, we get

$$\widetilde{{q}_{{\rm D}}}\left(s\right)=\frac{2\pi {y}_{{\rm e}}}{\sqrt{{A}_{{\rm cw}}}}\widetilde{{q}_{{\rm D}}^{*}}\left(s\right)$$

(75)

or in physical space,

$${q}_{{\rm D}}\left(t\right)=\frac{2\pi {y}_{{\rm e}}}{\sqrt{{A}_{{\rm cw}}}}{q}_{{\rm D}}^{*}\left(t\right)$$

(76)

Consequently, we get back to the dimensionless-rate relationship [Eq. (31)] by the derivation involving the rate solutions of both models, which proves the equivalence between Eq. (24) with $\left(\alpha , {\lambda }_{{\rm D}}, c\right)=\left(1.0, 0.0, 0.0\right)$ and Eq. (2A) in Bello and Wattenbarger (2010).

Bello (2009) extends the above reference rate solution to account for convergence skin. After a bit reorganization, Eq. (6.18) in Bello (2009) could be written as

$$\frac{1}{\widetilde{{q}_{{\rm D}}^{*}}\left(s\right)}=2\pi {s}^{*}\left({s}_{{\rm Ac}}+\frac{\mathrm{coth}\left(\sqrt{{s}^{*} f\left({s}^{*}\right)}\frac{{y}_{{\rm e}}}{\sqrt{{A}_{{\rm cw}}}}\right)}{\sqrt{{s}^{*} f\left({s}^{*}\right)}}\right)$$

(77)

Based on the previous equivalence, it is straightforward to obtain that Eq. (77) is equivalent to Eq. (24) with $\left(\alpha , {\lambda }_{{\rm D}}\right)=\left(1.0, 0.0\right)$ if we take ${s}_{{\rm Ac}}=c$. That is, even though the formulas of convergence skin factor are different due to the different assumptions of the fracture continuum shape, this factor has exactly the same effect in both the proposed and reference models.