Open Access
Numéro
Oil Gas Sci. Technol. – Rev. IFP Energies nouvelles
Volume 75, 2020
Numéro d'article 68
Nombre de pages 20
DOI https://doi.org/10.2516/ogst/2020062
Publié en ligne 9 octobre 2020

© R. Raghavan & C.C. Chen, published by IFP Energies nouvelles, 2020

Licence Creative CommonsThis is an Open Access article distributed under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Nomenclature

a : Exponent; see equation (A.39)

a ̃ $ \mathop{a}\limits^\tilde$ : See equation (A.41)

CFD : Dimensionless fracture conductivity

B : Formation volume factor [L3/L3]

c : Compressibility L T2/M

Eα,β(z): Mittag-Leffler function; see equation (24)

f ( s ̃ ) $ f(\mathop{s}\limits^\tilde)$ : See equation (15)

g ( η ̃ f ) $ g({\stackrel{\tilde }{\eta }}_f)$ : See equation (21)

H p , q m , n [ x ( a 1 , A 1 ) , , ( a p , A p ) ( b 1 , B 1 ) , , ( b q , B q ) ] $ {H}_{p,q}^{m,n}\left[x,\begin{array}{c}\left({a}_1,{A}_1\right),\dots,\left({a}_p,{A}_p\right)\\ \left({b}_1,{B}_1\right),\dots,\left({b}_q,{B}_q\right)\end{array}\right]$ : Fox function

ht : Thickness [L]

IRNP(te): Integral of rate-normalized pressure

k : Permeability [L2]

kα,β : See equation (3)

L : Width of linear system [L]

Mν(z): Mainardi or M-Wright function

p : Pressure [M/L/T2]

q : Flux [L3/T]

RNP(te): Rate-normalized pressure

T : Absolute temperature

t : Time [T]

te : Material balance time [T]

W(α, β; t): Wright function

w : Width of fracture [L]

xL : Distance to interface [L]

α : Exponent, time

β : Exponent, space

βf : See equation (A.5)

β : See equation (A.18)

β ⋆⋆ : See equation (A.29)

Γ(x): Gamma function

η : Diffusivity; various

η ̃ f $ {\stackrel{\tilde }{\eta }}_f$ : Equation (18) [ L β f + 1 / T f α $ {\mathrm{L}}^{{\beta }_f+1}/{\mathrm{T}}_f^{\alpha }$]

λ : Mobility [L T/M]

λ′: See equation (19)

μ : Viscosity [M/L/T]

ς : Similarity variable; see equation (64)

ϕ : Porosity [L3/L3]

ω′: See equation (20)

Ω: Domain

Subscripts

c: Constant rate

D : Dimensionless

f : Fracture system

F : Hydraulic fracture

N : The number of fractures

p : Constant pressure

m : Matrix system

t : Total

1, 2: Index

Superscript

–: Laplace transform

1 Introduction

The vast array of multiple mathematical models proposed to evaluate fractured reservoirs attest to their complexity. It is generally agreed that fluid is transported along preferred pathways with the pathways supported by a medium of lower conductivity. Both the pathways and lower conductivity system are complex geologically and in terms of geological gradients. Some of the complexity is usually better understood subsequent to the drilling of new (additional) wells; see Raghavan (2004). This view is a result of both mapping and a well-defined, spatiotemporal evolution of responses at wells. Broadly speaking, the types of models used today may be classified in two ways: (i) those that are known as equivalent porous media models (Yeh et al., 2015) as discussed in Pruess and Narasimhan (1985), Haggerty and Gorelick (1995), and Noetinger et al. (2001) which are extensions of the work of Barenblatt et al. (1960) and Warren and Root (1963) and, (ii) those that are known as Distributed Fracture Network (DFN) schemes where fractures are introduced as discontinuities with specific properties (orientation, aperture, size, position, etc.) in a uniform continuum and explicitly model pathways representing natural fractures (Artus, 2020; Cacas et al., 1990a, b; Dershowitz et al., 2002; Karimi-Fard et al., 2004; Noetinger, 2015). The equivalence of the two approaches is compared in Dong et al. (2019) and the circumstances in which the DFN network may be treated as a multiple, interacting continuum is noted there. Irrespective of the scheme adopted, one may not be indifferent to the manifold aspects of the transport of fluids and treat the system as an undifferentiated continuum concerning properties. Raghavan et al. (2001) elucidate the complexities and intrinsic connections of addressing the origins of reservoir heterogeneity and demonstrate the prodigious amount of work involved in dataset curation. Here as noted in Cacas et al. (2001) we must integrate of matters relating to fracture permeability, fracture aperture and fracture connectivity to mention only a few aspects; see Belayneh et al. (2006), Kang et al. (2015) and Bisdom et al. (2016). One may attempt to address these manifold features through either of the models noted above. Here our focus is on the use of equivalent-porous-media concepts based on anomalous diffusion by fractional calculus (Albinali et al., 2016; Fomin et al., 2011; Holy and Ozkan, 2016; Raghavan, 2011). In this scheme, rather than following the classical time-distance relationship, r 2 ~ t, with r being the distance and t being time; the transient follows the relation r 2 ~ t a where a is a constant with a < 1 (subdiffusion) or with a > 1 (superdiffusion), or even with a = 1 (combined influences). Chu et al. (2019, 2020) discuss the advantages of addressing atypical behaviors employing a subdiffusive model to evaluate wells of the Permian basin produced through multiple hydraulic fractures (for both single and multiple well-tests). Quantitative estimates of the exponent a for the tests discussed in Chu et al. (2019) are given in Raghavan and Chen (2018b). Our experience, confirmed by Chu (2018), and illustrated with an example here indicates that one must also consider superdiffusive flows in the fracture system. Significant challenges exist in addressing superdiffusive problems analytically; in fact, this is a fraught exercise. For example, Defterli et al. (2015) observe that this is an open problem and Metzler et al. (2007) note that the method of images may not be used to address reflection conditions; that is, closed boundary effects may not be modelled by image wells. One way to circumvent this issue is by considering approximate or asymptotic cases as we demonstrate here. Similarly, obtaining solutions of interest for 2D problems is difficult.

On the basis of O’Shaughnessy and Procaccia (1985a, b), many studies have proposed solutions for subdiffusive, transient flows based on rules for fractal objects (Chang and Yortsos, 1990). For vertical wells, such solutions lead to an expression where the pressure drop, Δp(r, t), may be expressed as a complementary incomplete gamma function: Γ ( a , x ) = x ζ a - 1 exp ( - ζ ) d ζ $ \mathrm{\Gamma }(a,x)={\int }_x^{\mathrm{\infty }} {\zeta }^{a-1}\mathrm{exp}(-\zeta )\mathrm{d}\zeta $ with the principal conclusion suggesting that at long enough times transients display characteristic power-law features, namely, Δp ~ t ϑ , where ϑ is a constant, contrary to the widely known result Δp ~ lnt. The finding that Δp ~ Γ(ax) may, however, be obtained in a number of ways without invoking the rules for fractal behavior, such as by assuming that the permeability, k, and porosity, ϕ, follow specific forms as in Carslaw and Jaeger (1959) or by assuming that the flow dimension is different from the Cartesian dimension as in Barker (1988). Incidentally, physical meaning may be ascribed to Barker’s work only in terms of a fractal model (Raghavan, 2008) or through a fractional differential equation. Conditions required by fractal models, specifically symmetry, are difficult to satisfy in certain circumstances for problems of interest; see Beier (1994); such difficulties, however, may be circumvented through fractional calculus. Also, Jourde et al. (2002) and Doe (1991) note that transient responses display characteristic power-law features when conduits are not characterized as fractals; a similar observation for fluvial deposits may be found in Thomas et al. (2005). Fractional laws proposed in this paper for both sub- and super-diffusion preclude conjuring counterfactual geological options.

Our objective is to address patterns of derivative responses shown in Figure 1, processed and presented by Doe et al. (2014). Early-time responses, as is usual in most wells producing unconventional reservoirs, reflect the lack of appropriate measurements particularly bottom-hole pressures; see Scott et al. (2015). The principal point to note here is that all responses at late times display a trend well above what would expect of models based on classical diffusion. None of the responses display slopes of 1/2, 1/4, or 1/8 and it is not appropriate to force a straight line fit with such slopes. Our point is that if we are trying to address heterogeneous systems, we are required to go beyond classical systems. We offer one choice to do so. Regardless of the approach taken, three elements, as incorporated here, are paramount: a matrix system, a fracture or fissure system and a hydraulic fracture.

thumbnail Fig. 1

Rate-normalized responses at five wells producing the same lease in the Eagle Ford Group; after Doe et al. (2014).

As in the classical case, we develop the model by combining the conservation of mass principle with the equation of state for a slightly compressible fluid and a fractional flux law. It is needless to emphasize the importance of a flux law as it is the marrow of reservoir-engineering predictions. Even in the classical case, the deterministic situation, a model derived by the three principles just noted, the basic driving force is attributed to random motions namely the random walk (Doyle and Snell, 1984; Klafter and Sokolov, 2011) that consists of independent and identically distributed steps, Xn; the first moment (the conservation of mass principle) and second moment (the time-distance relationship) being finite. Inherent in the classical Darcy law is the concept that the porous rock is a uniform body as the overall size of the porous body is much, much larger than the sizes the individual pores that constitute the porous rock and that any potential surface is a notional surface in that a particular cross-sectional element is not unlike any other and the normal gradient reflects the overall flux. If, however, we postulate that this is no longer the case, the system is heterogeneous and distant points connected, then modification is required; one such option being to consider anomalous diffusion with random jumps having a power law probability tail in which case we use an equation of the form t p = x β + 1 p $ {\mathrm{\partial }}_tp={\mathrm{\partial }}_x^{\beta +1}p$ for 0 < β < 11 with β reflecting the power law tail ( P ( | X n | > r ) r - ( β + 1 ) $ \mathbb{P}(|{X}_n|>r)\approx {r}^{-(\beta +1)}$); see Benson et al. (2000, 2004). Levy-stable-particles trace a path of dimension β + 1 namely the exponent of the x-derivative, ∂x, in t p = x β + 1 p $ {\mathrm{\partial }}_tp={\mathrm{\partial }}_x^{\beta +1}p$; in other words, the Hurst index, H, is 1/(β + 1) or 〈x 2〉 ~ t H . On the other hand, if we were to deem restrictions, impediments, intercalations, etc. exist, then we need to model subdiffusion rather than superdiffusion; that is, consider solutions of t α p = x 2 p $ {\mathrm{\partial }}_t^{\alpha }p={\mathrm{\partial }}_x^2p$ for 0 < α < 1 rather than t p = x β + 1 p $ {\mathrm{\partial }}_tp={\mathrm{\partial }}_x^{\beta +1}p$; in this case P ( T n > t ) t - α $ \mathbb{P}({T}_n>t)\approx {t}^{-\alpha }$ for 0 < α < 1, where Tn is the random waiting time at which the jump occurs. Such motions are known as a Continuous Time Random Walk (CTRW).2 Another way to arrive at this result is through waiting-time solutions as documented in Henry et al. (2010). In fact, the expression t α p = x 2 p $ {\mathrm{\partial }}_t^{\alpha }p={\mathrm{\partial }}_x^2p$ may be arrived at in a number of ways; see Montroll and Weiss (1965), Kenkre et al. (1973), Hilfer and Anton (1995). The exponent α reflects the anomalous diffusion coefficient (random walk dimension), dw (Metzler et al., 1994) for diffusion on fractal objects. The earth-science literature abounds with examples of situations where subdiffusive models apply, namely in situations wherein cracks, fissures and cervices abound (Caine et al., 1996; Evans, 1988; Mitchell and Faulkner, 2009; Savage and Brodsky, 2011; Scholz et al., 1993). Mathematical models depicting flows in porous rocks may be found in Jourde et al. (2002), Cortis and Knudby (2006) and Suzuki et al. (2016). In combined situations, we consider t α p = x β + 1 p $ {\mathrm{\partial }}_t^{\alpha }p={\mathrm{\partial }}_x^{\beta +1}p$ for 0 < α, β < 1 while recognizing that for combinations of α and β behaviors similar to classical situations may result (Chen and Raghavan, 2015; Magin et al., 2013). Excellent reviews that flesh out details may be found in Zhang et al. (2009) and Benson et al. (2013).

This work endeavors to understand the structure of the solutions of the model proposed here constrained by the premises that the geology is complex (Gale et al., 2014; Thomas et al., 2005), fluid is extracted through wellbores of complex shapes with multiple extraction points (Chen and Raghavan, 2013; Larsen and Hegre, 1991), and nonlocality by deriving a tractable, analytical solution. As noted in Ozkan and Raghavan (1991), the first step involves the development of point-source solutions that satisfy the specifications of the desired model. Should such a solution be available, as it is in many cases, prediction of the performance of a fractured well producing through a horizontal well through multiple fractures is tractable. Here, of course, we consider situations where fluids are transported in the fracture and matrix systems under anomalous diffusion, specifically, through fluid movement in the fracture system that permits for conditions instrumental in both super- and subdiffusive transport and rock properties in the matrix system conducive to subdiffusive flow.

Thus transient diffusion considered here is not construed to be diffusion in homogeneous systems in the manner envisaged by early models of naturally fractured reservoirs discussed in Barenblatt et al. (1960), Warren and Root (1963), Kazemi (1969) and de Swaan-O (1976) although we do consider responses at the identical scale. Probability density functions that result from solving fractional differential equations are different from those generated by the equation governing classical diffusion and are examined here.

2 The model

The mathematical model consists of three regions: the fracture system, Ωf, the matrix system, Ωm, and hydraulic fracture, ΩF. Figure 2 taken from Raghavan and Ohaeri (1981) is a portrayal of the model put forward in de Swaan-O (1976) which shows the colligation of the matrix and fracture elements where hf denotes the height of a fracture element and hm denotes the height of a matrix element. As we consider n elements, the total thickness, ht, is

h t = h ft + h mt , = n h f + n h m . $$ \begin{array}{c}{h}_t={h}_{{ft}}+{h}_{{mt}},\\ =n{h}_f+n{h}_m.\end{array} $$(1)

thumbnail Fig. 2

Schematic of the transient, interporosity model of Kazemi (1969) and de Swaan-O (1976); after Raghavan and Ohaeri (1981).

Although we consider linear elements, we emphasize that the results given here apply to situations where the matrix elements are in the form of spheres and cylinders as noted in Raghavan and Ohaeri (1981) except in a minor way during transitional periods. This observation is reconfirmed in Chen (1982) and Carrera et al. (1998).

Depletion of the matrix system, Ωm, occurs in the usual way as it is in perfect contact with the fracture system Ωf (Cinco-Ley and Meng, 1988; Houze et al., 1988; Yeh et al., 1986). All fluid produced through the hydraulic fracture may enter the hydraulic fracture only through the fracture system; that is, there is no communication between the matrix system and the hydraulic fracture. Perfect contact is assumed at each interface where there is hydraulic communication. The fracture tip is assumed to be sealed. As noted earlier, our representation of the fractured rock incorporates both super- and subdiffusive effects to model flow in the fracture system while flow in the matrix is presumed to be subdiffusive. Darcy flow is assumed in the hydraulic fracture. To consider a system of multiple fractures intercepting a horizontal wellbore and producing against a common wellbore pressure, we follow the outlines of Chen and Raghavan (2013) through Duhamel’s theorem.

2.1 The flux law

As noted in the Introduction we employ fractional flux laws to address heterogeneity. For subdiffusive flows, fractional-flux laws were proposed over 30 years ago; see Le Mẽhautẽ and Crepy (1983) and Nigmatullin (1984, 1986). As discussed in Raghavan (2011) the need for a fractional flux-law is implied in the development in Metzler et al. (1994) which improves upon the work of O’Shaughnessy and Procaccia (1985a) for transient diffusion in fractal structures; see also Albinali et al. (2016) and Suzuki et al. (2016). As noted earlier, ideas based on Continuous Time Random Walks (Hilfer and Anton, 1995; Kenkre et al., 1973; Montroll and Weiss, 1965; Noetinger and Estebenet, 2000) and waiting-time solutions (Henry et al., 2010) also yield such laws. For superdiffusive flows, discussions are available in Benson et al. (2000, 2004), Molz et al. (2002), Garra and Salusti (2013) and Kim et al. (2015). More recently, detailed treatments regarding experimental justification of the anomalous diffusion concept are discussed in Tao et al. (2016), Yanga et al. (2019) and Zhokh and Strizhak (2018, 2019). Many discourses on fractional Fourier laws are also found in the literature (Angulo et al., 2000; Dassas and Duby, 1995; Gurtin and Pipkin, 1968; Moodie and Tait, 1983; Norwood, 1972; Su, 2014; Su et al., 2015).

The fractional flux law discussed in Fomin et al. (2011), Magin et al. (2013), and Chen and Raghavan (2015) uses the fractional derivative, ∂ α /∂t α f(t), in the form defined by Caputo (1967)

α t α f ( t ) = 1 Γ ( 1 - α ) 0 t d t' ( t - t' ) - α t' f ( t' ) , $$ \frac{{\mathrm{\partial }}^{\alpha }}{\mathrm{\partial }{t}^{\alpha }}f(t)=\frac{1}{\mathrm{\Gamma }\left(1-\alpha \right)}{\int }_0^t \mathrm{d}{t\prime}(t-{t\prime}{)}^{-\alpha }\frac{\mathrm{\partial }}{\mathrm{\partial }{t\prime}f({t\prime}), $$(2)and assumes to govern transients in the fractured rock by an expression of the form for the velocity, v(x, t):

v ( x , t ) = - k α , β μ 1 - α t 1 - α [ β p x β ] , $$ v(x,t)=-\frac{{k}_{\alpha,\beta }}{\mu }\frac{{\mathrm{\partial }}^{1-\alpha }}{\mathrm{\partial }{t}^{1-\alpha }}\left[\frac{{\mathrm{\partial }}^{\beta }p}{\mathrm{\partial }{x}^{\beta }}\right], $$(3)with the exponents α and β both reflecting the heterogeneities in the fractured rock.3 Note that unlike most studies a one-sided form of the spatial fractional derivative is used. Denoting the exponents of the fracture system by the subscript f and that of the matrix system by the subscript m, αf and βf are presumed to be ≤1 as we consider both sub- and superdiffusive influences in the fracture system. For the matrix system, αm is presumed to be ≤1 while βm is taken to be 1 because only subdiffusive influences are considered.

2.2 The partial differential equations, associated conditions

The three partial differential equations to be solved to obtain the pressure distributions in the fracture system (Region Ωf), matrix system (Region Ωm) and the hydraulic fracture (Region ΩF) for t > 0, are respectively:

( ϕ c ) f μ k ̃ f Δ p f t = [ 1 - α f t 1 - α f ( β f Δ p f ) ] + 2 k ̃ m k ̃ f h f [ 1 - α m t 1 - α m ( Δ p m z ) ] z = h m 2 , $$ \begin{array}{ll}\frac{({\phi c}{)}_f\mu }{{\mathop{k}\limits^\tilde}_f}\frac{\mathrm{\partial \Delta }{p}_f}{\mathrm{\partial }t}& =\nabla \left[\frac{{\mathrm{\partial }}^{1-{\alpha }_f}}{\mathrm{\partial }{t}^{1-{\alpha }_f}}({\nabla }^{{\beta }_f}\Delta {p}_f)\right]\\ & +\frac{2{\mathop{k}\limits^\tilde}_m}{{\mathop{k}\limits^\tilde}_f{h}_f}{\left[\frac{{\mathrm{\partial }}^{1-{\alpha }_m}}{\mathrm{\partial }{t}^{1-{\alpha }_m}}\left(\frac{\mathrm{\partial \Delta }{p}_m}{\mathrm{\partial }z}\right)\right]}_{z=\frac{{h}_m}{2}},\end{array} $$(4)

( ϕ c ) m μ k ̃ m α m t α m Δ p m ( z , t ) = 2 Δ p m ( z , t ) z 2 , $$ \frac{({\phi c}{)}_m\mu }{{\mathop{k}\limits^\tilde}_m}\frac{{\mathrm{\partial }}^{{\alpha }_m}}{\mathrm{\partial }{t}^{{\alpha }_m}}\Delta {p}_m(z,t)=\frac{{\mathrm{\partial }}^2\Delta {p}_m(z,t)}{\mathrm{\partial }{z}^2}, $$(5)and

ϕ F c F μ k F p F ( x ; t ) t = x [ p F ( x ; t ) x ] + y [ p F ( x ; t ) y ] . $$ \frac{{\phi }_F{c}_F\mu }{{k}_F}\frac{\mathrm{\partial }{p}_F\left({x};t\right)}{\mathrm{\partial }t}=\frac{\mathrm{\partial }}{\mathrm{\partial }x}\left[\frac{\mathrm{\partial }{p}_F\left({x};t\right)}{\mathrm{\partial }x}\right]+\frac{\mathrm{\partial }}{\mathrm{\partial }y}\left[\frac{\mathrm{\partial }{p}_F\left({x};t\right)}{\mathrm{\partial }y}\right]. $$(6)

Here pi the initial pressure is identical in all parts of the system at t = 0, and Δpj = pipj(xt); j being f, m or F, k ̃ f = k α f , β f $ {\mathop{k}\limits^\tilde}_f={k}_{{\alpha }_f,{\beta }_f}$, k ̃ m = k ̃ α m $ {\mathop{k}\limits^\tilde}_m={\mathop{k}\limits^\tilde}_{{\alpha }_m}$. We require a solution to equations (4)(6) subject to the following conditions. For a quiescent system where equilibrium prevails everywhere we require

Δ p j ( x , t ) = 0 ;   for   t = 0 . $$ \Delta {p}_j(x,t)=0;\enspace \mathrm{for}\enspace t=0. $$(7)

For the matrix elements we require that the boundary at z = 0 be a no-flow boundary; thus

1 - α m t 1 - α m [ Δ p m x ] = 0 . $$ \frac{{\mathrm{\partial }}^{1-{\alpha }_m}}{\mathrm{\partial }{t}^{1-{\alpha }_m}}\left[\frac{\mathrm{\partial \Delta }{p}_m}{\mathrm{\partial }x}\right]=0. $$(8)

Further, for perfect contact we require

Δ p f ( x , t ) = Δ p m ( h m / 2 , t ) , $$ \Delta {p}_f({x},t)=\Delta {p}_m({h}_m/2,t), $$(9)as the requirement of continuity in pressure at the fracture-matrix interface, Ωf − Ωm, that is, at (z = hm/2). To satisfy the requirement that there is contact we also require that the flux at the interface Ωf − Ωm be continuous. This stipulation, continuity of flux, was incorporated in deriving equation (4); it appears as the second term on the right-hand side of the equation.

The constraints for the pressure distribution are as follows. All boundaries are assumed to be sealed. Hence we have

1 - α f t 1 - α f [ x β f ( Δ p f ) ] = 0   at   x = 0   and   at   x = x e , $$ \frac{{\mathrm{\partial }}^{1-{\alpha }_f}}{\mathrm{\partial }{t}^{1-{\alpha }_f}}\left[{\nabla }_x^{{\beta }_f}(\Delta {p}_f)\right]=0\enspace \mathrm{at}\enspace x=0\enspace \mathrm{and}\enspace \mathrm{at}\enspace x={x}_e, $$(10)and

1 - α f t 1 - α f [ y β f ( Δ p f ) ] = 0   at   y = 0   and   at   y = y e . $$ \frac{{\mathrm{\partial }}^{1-{\alpha }_f}}{\mathrm{\partial }{t}^{1-{\alpha }_f}}\left[{\nabla }_y^{{\beta }_f}(\Delta {p}_f)\right]=0\enspace \mathrm{at}\enspace y=0\enspace \mathrm{and}\enspace \mathrm{at}\enspace y={y}_e. $$(11)

The conditions to be satisfied at the fracture-reservoir interface, Ωf – ΩF, are, again, the continuity of pressure and the continuity of flux; in consequence,

p f ( x , t ) = p F ( x ; t ) ;   at   the   interface   Ω f - Ω F , $$ {p}_f\left({x},t\right)={p}_F\left({x};t\right);\enspace \mathrm{at}\enspace \mathrm{the}\enspace \mathrm{interface}\enspace {\mathrm{\Omega }}_f-{\mathrm{\Omega }}_F, $$(12)and

q ̃ f ( x , t ) = q ̃ F ( x , t ) ;   at   the   interface   Ω f - Ω F , $$ {\mathop{q}\limits^\tilde}_f\left({x},t\right)={\mathop{q}\limits^\tilde}_F\left({x},t\right);\enspace \mathrm{at}\enspace \mathrm{the}\enspace \mathrm{interface}\enspace {\mathrm{\Omega }}_f-{\mathrm{\Omega }}_F, $$(13)respectively, where the symbol q ̃ $ \mathop{q}\limits^\tilde$ represents the flux.

As is usual we follow the outline of the idealization given in Cinco-Ley and Meng (1988) to consider flow in the hydraulic fracture; that being so, the following assumptions apply: (i) the tips of the hydraulic fracture are sealed, (ii) the width of the fracture is negligibly small compared with its length; thus, one may use the pseudofunction approximation for flow the direction normal to its length, and (iii) the volume in the fracture is inconsequential compared with the reservoir volume and therefore we presume flow is steady in the fracture. The assumptions considered hitherto lead to the consideration that the fracture is a source term but of finite conductivity, kF, and finite width wF. For large values of the product kF wF we may consider a planar fracture of infinite conductivity (Prats, 1961). In addition, we assume that fluid is withdrawn from the center of the fracture.

2.2.1 Coupling of partial differential equations

We now briefly address the coupling of equations (4)(6). The simplest manner to couple these three equations is through the Laplace transformation. On combining equations (4) and (5) we have in terms of the Laplace transformation with respect to t the expression

( β f Δ p ̅ f ) - s ̂ f ( s ̃ ) Δ p ̅ f = 0 , $$ \nabla \left({\nabla }^{{\beta }_f}{\overline{\Delta p}}_f\right)-\widehat{s}f(\mathop{s}\limits^\tilde){\overline{\Delta p}}_f=0, $$(14)where

f ( s ̃ ) = 1 + λ ' 3 s ̃ g ( η ̃ f ) tanh [ 3 ω ' s ̃ g ( η ̃ f ) λ ] . $$ f(\mathop{s}\limits^\tilde)=1+\sqrt{\frac{{\lambda \prime\omega \prime}{3\mathop{s}\limits^\tildeg({\stackrel{\tilde }{\eta }}_f)}}\mathrm{tanh}\left[\sqrt{\frac{3{\omega \prime}\mathop{s}\limits^\tildeg({\stackrel{\tilde }{\eta }}_f)}{{\lambda }^\mathrm{\prime}}}\right]. $$(15)

s ̂ = s α f η ̃ f , $$ \widehat{s}=\frac{{s}^{{\alpha }_f}}{{\stackrel{\tilde }{\eta }}_f}, $$(16)

s ̃ = s α m η ̃ f , $$ \mathop{s}\limits^\tilde=\frac{{s}^{{\alpha }_m}}{{\stackrel{\tilde }{\eta }}_f}, $$(17)and

η ̃ f = k α f , β f ϕ f c f μ , $$ {\stackrel{\tilde }{\eta }}_f=\frac{{k}_{{\alpha }_f,{\beta }_f}}{{\phi }_f{c}_f\mu }, $$(18)with l $ \mathcal{l}$ being the reference length. The symbols λ′ and ω′ are the analogues of λ and ω, respectively, used in the Barenblatt et al. (1960) or Warren and Root (1963) models except that they apply to the model in Figure 2 and are expressed in terms of intrinsic properties. The expressions for λ′ and ω′ are defined, respectively, as

λ = 12 ( k ̃ m h m k ̃ f h f ) ( l ref β + 1 h m 2 ) ( η ̃ f l ref β f + 1 ) α f - α m α f , $$ \lambda \mathrm{\prime}=12\left(\frac{{\mathop{k}\limits^\tilde}_m{h}_m}{{\mathop{k}\limits^\tilde}_f{h}_f}\right)\left(\frac{{\mathcal{l}}_{{ref}}^{\beta +1}}{{h}_m^2}\right){\left(\frac{{\stackrel{\tilde }{\eta }}_f}{{\mathcal{l}}_{{ref}}^{{\beta }_f+1}}\right)}^{\frac{{\alpha }_f-{\alpha }_m}{{\alpha }_f}}, $$(19)and

ω = ϕ m c m h mt ϕ f c f h ft = ϕ m c m h m ϕ f c f h f , $$ \omega \mathrm{\prime}=\frac{{\phi }_m{c}_m{h}_{{mt}}}{{\phi }_f{c}_f{h}_{{ft}}}=\frac{{\phi }_m{c}_m{h}_m}{{\phi }_f{c}_f{h}_f}, $$(20)where

g ( η ̃ f ) = ( η ̃ f l β f + 1 ) α f - α m α f . $$ g({\stackrel{\tilde }{\eta }}_f)={\left(\frac{{\stackrel{\tilde }{\eta }}_f}{{\mathcal{l}}^{{\beta }_f+1}}\right)}^{\frac{{\alpha }_f-{\alpha }_m}{{\alpha }_f}}. $$(21)

The solution of equation (14) and its associated boundary conditions may then be coupled with the solution of equation (6) given by

Δ p ̅ w - Δ p ̅ f ( x , y w ) = 1 2 k ̃ f h ft C FD s [ | x ̃ - x w | + s 0 | x ̃ - x w | 0 x q ̃ ¯ ( x ) d x d x ] , $$ {\overline{\Delta p}}_w-{\overline{\Delta p}}_f(x,{y}_w)=\frac{1}{2{\mathop{k}\limits^\tilde}_f{h}_{{ft}}{C}_{{FD}}^{\star }s}\left[|\mathop{x}\limits^\tilde-{x}_w|+s{\int }_0^{|\mathop{x}\limits^\tilde-{x}_w|} {\int }_0^{x\mathrm{\prime} \overline{\mathop{q}\limits^\tilde}(x\mathrm{\prime}\mathrm{\prime})\mathrm{d}x\mathrm{\prime}\mathrm{\prime}\mathrm{d}x\mathrm{\prime}\right], $$(22)where | x ̃ - x w | L F $ |\mathop{x}\limits^\tilde-{x}_w|\le {L}_F$, and x ̃ $ \mathop{x}\limits^\tilde$ refers to the x coordinate along the fracture plane. The coupling noted above occurs through equation (12) and equation (13) at the interface Ωf − ΩF if we note that Δ p ̅ f ( x , y w ) $ {\overline{\Delta p}}_f(x,{y}_w)$ and q ̃ ¯ ( x ) $ \overline{\mathop{q}\limits^\tilde}(x\mathrm{\prime}\mathrm{\prime})$ are obtained from the solution of equation (14) and its associated boundary conditions. Using this procedure just outlined we arrive at the response at a well producing through a single vertical fracture for either the pressure drop, Δpw(t) if the well were produced at a constant rate, q, or the production rate, q(t), if the well were produced at a constant pressure, Δpw. The response for that of a horizontal well intercepted by multiple hydraulic fractures and producing against a common wellbore pressure is discussed below (Fig. 3).

thumbnail Fig. 3

Schematic: alignment of hydraulic fractures with wellbore.

2.2.2 The function f ( s ̃ ) $ f(\mathop{s}\limits^\tilde)$

Regardless of the approach taken (classical diffusion, fractional diffusion, CTRW, etc.) one ultimately arrives at the fact that the response of fractured reservoirs depends on a function similar to f(s); see Dontsov and Boyrchuk (1971), de Swaan-O (1976), Chen (1982), Ozkan and Raghavan (1991), Cinco-Ley and Meng (1988), Noetinger et al. (2001, 2016) and Raghavan and Chen (2017). Distinctive features of well responses are governed by the variation of tanh [ 3 ω s ̃ g ( η ̃ f ) / λ ] $ \mathrm{tanh}\left[\sqrt{3{\omega }^\mathrm{\prime}\mathop{s}\limits^\tildeg({\stackrel{\tilde }{\eta }}_f)/{\lambda }^\mathrm{\prime}}\right]$ with respect to s (time). Three options are available: For s→ ∞; early times specified as Flow Regime 1, ( f ( s ̃ ) 1 $ f(\mathop{s}\limits^\tilde)\to 1$ and tanh(x) ≈ 1; matrix has no influence), for finite values of s such that tanh(x) ≈ 1 and 3 ω s ̃ g ( η ̃ f ) / λ 1 $ \sqrt{3{\omega }^\mathrm{\prime}\mathop{s}\limits^\tildeg({\stackrel{\tilde }{\eta }}_f)/{\lambda }^\mathrm{\prime}}\gg 1$; intermediate times specified as Flow Regime 2, ( f ( s ̃ ) 3 ω s ̃ g ( η ̃ f ) / λ $ f(\mathop{s}\limits^\tilde)\approx \sqrt{3{\omega }^\mathrm{\prime}\mathop{s}\limits^\tildeg({\stackrel{\tilde }{\eta }}_f)/{\lambda }^\mathrm{\prime}}$; flow regime in the matrix system reflects transient conditions), and for s → 0; late times specified as Flow Regime 3, ( f ( s ̃ ) 1 + ω $ f(\mathop{s}\limits^\tilde)\to 1+\omega \mathrm{\prime}$ as tanh(x) ≈ 1/x; flow regime in the matrix system reflects pseudosteady conditions). These three conditions in order of appearance will exhibit slopes of m (early times), m/2 (intermediate times), and m (late times) at the well provided that boundary effects are negligibly small. The term boundary, incidentally, in the context of the transient-flow-model, includes the dashed line shown in Figure 2. Should the boundaries of the fracture system dominate the well response and the flow regime in the matrix system reflect the transient state, then Flow Regime 4, the counterpart of Flow Regime 2 first identified in Raghavan and Ohaeri (1981), will be evident. Flow Regime 4 may be readily identified through its characteristic slope as we discuss below. Finally should pseudosteady conditions prevail everywhere, then the counterpart of Flow Regime 3, namely Flow Regime 5, will be reflected (Tab. 1).

Table 1

Flow regimesa,b,c,d.

In concluding this section, we should note that additional options exist to examine solutions; for example, Chen (1982) considers f ( s ) 1 + 3 ω ' s ̃ g ( η ̃ f ) / λ ' $ f(s)\approx 1+\sqrt{3{\omega }^{\prime}\mathop{s}\limits^\tildeg({\stackrel{\tilde }{\eta }}_f)/{\lambda }^{\prime}$; this aspect is not examined here.

2.3 Source solutions

As mentioned earlier, if instantaneous source solutions were available, then the analytical option becomes an amenable proposition to obtain pressure distributions of interest. Further, we noted that no simple analytical methods are available to obtain solutions of interest for bounded systems should β < 1; that is, if superdiffusive conditions are to be modeled. Here we use options that are available and examine an instantaneous plane-source-solution in an infinite system which incorporates both sub- and superdiffusive considerations. Let γ( x ) be the instantaneous source function, located at x ′, that acts at t = 0, with γ( x ) = 0, for t < 0.

2.3.1 Plane source

The plane-source solution is particularly useful in obtaining early time responses for a fractured well. Thus, considering a linear system we may derive the plane-source solution along the lines outlined in Carslaw and Jaeger (1959). The required expression is given by

γ ̅ ( x ) = 1 2 η ̃ f 1 s 1 - α f [ 1 ( s α f f ( s ̃ ) / η ̃ f ) β f / ( β f + 1 ) E β f + 1 ( s α f η ̃ f f ( s ̃ ) x β f + 1 ) - x f β E β f + 1 , β f + 1 ( s α f η ̃ f f ( s ̃ ) x β f + 1 ) ] . $$ \begin{array}{ll}\overline{\gamma }(x)& =\frac{1}{2{\stackrel{\tilde }{\eta }}_f}\frac{1}{{s}^{1-{\alpha }_f}}\left[\frac{1}{({s}^{{\alpha }_f}f(\mathop{s}\limits^\tilde)/{\stackrel{\tilde }{\eta }}_f{)}^{{\beta }_f/({\beta }_f+1)}}{E}_{{\beta }_f+1}\left(\frac{{s}^{{\alpha }_f}}{{\stackrel{\tilde }{\eta }}_f}f(\mathop{s}\limits^\tilde){x}^{{\beta }_f+1}\right)\right.\\ & -\left.{x}_f^{\beta }{E}_{{\beta }_f+1,{\beta }_f+1}\left(\frac{{s}^{{\alpha }_f}}{{\stackrel{\tilde }{\eta }}_f}f(\mathop{s}\limits^\tilde){x}^{{\beta }_f+1}\right)\right].\end{array} $$(23)

In the above expression, Eα,β(z) is the two-parameter Mittag-Leffler function (Erdelyi et al., 1955; Mainardi, 2010; Uchaikin, 2013; and Povstenko, 2015):

E α , β ( z ) = n = 0 z n Γ ( β + α n ) ,   α ,   β > 0 ,   z     C , $$ {E}_{\alpha,\beta }(z)=\sum_{n=0}^{\mathrm{\infty }} \frac{{z}^n}{\mathrm{\Gamma }(\beta +{\alpha n})},\enspace \hspace{1em}\alpha,\enspace \beta >0,\enspace \hspace{1em}z\enspace \in \enspace C, $$(24)where C represents the complex plane, and the symbol Γ(x) is the gamma function. Commonly used functions such as the exponential, exp(z), the hyperbolic functions, sinh(z), cosh(z), and functions involving the error function complement, erfc(z), may be represented in order as: E1(z) = exp(z), E2,1(z 2) = cosh(z), E2,2(z 2) = sinh(z)/z, and exp(z 2)erfc(−z) = E1/2,1(z). Uchaikin (2013) provides convenient tabulations of expressions for obtaining inverse Laplace transforms containing Eα,β(z).

For βf = 1, the right-hand side of equation (23) reduces to an expression in terms of a stretched exponential given by

γ ̅ ( x ) = 1 2 η ̃ f 1 s 1 - α 1 λ ̃ exp ( - λ ̃ x D ) , $$ \overline{\gamma }(x)=\frac{1}{2{\stackrel{\tilde }{\eta }}_f}\frac{1}{{s}^{1-\alpha }}\frac{1}{\sqrt{\stackrel{\tilde }{\lambda }}}\mathrm{exp}\left(-\sqrt{\stackrel{\tilde }{\lambda }}{x}_D\right), $$(25)where

λ ̃ = s α f η ̃ f f ( s ̃ ) . $$ \stackrel{\tilde }{\lambda }=\frac{{s}^{{\alpha }_f}}{{\stackrel{\tilde }{\eta }}_f}f(\mathop{s}\limits^\tilde). $$(26)

The right-hand side of equation (25) may be inverted either through the Fox function or H-function, H p , q m , n [ x | ( b 1 , B 1 ) , . . . , ( b q , B q ) ( a 1 , A 1 ) , . . . , ( a p , A p ) ] $ {H}_{p,q}^{m,n}\left[x{|}_{({b}_1,{B}_1),...,({b}_q,{B}_q)}^{({a}_1,{A}_1),...,({a}_p,{A}_p)}\right]$, by using the relation

L - 1 [ s - ρ exp ( - z s ς ) ] = t ρ - 1 H 1,1 1,0 [ z t - ς | ( 0,1 ) ( ρ , ς ) ] , $$ {L}^{-1}[{s}^{-\rho }\mathrm{exp}(-z{s}^{\varsigma })]={t}^{\rho -1}{H}_{\mathrm{1,1}}^{\mathrm{1,0}}\left[z{t}^{-\varsigma }{|}_{(\mathrm{0,1})}^{(\rho,\varsigma )}\right], $$(27)where R ( s ) >   0 $ \mathbb{R}(s)>\enspace 0$, ς > 0 as noted in Saxena et al. (2006), or through the Wright function, W(α, β; t), as indicated in Povstenko (2015) by the expression

L - 1 [ s - β exp ( - λ s α ) ] = t β - 1 W ( - α , β ; - λ t - α ) ,   0 < α < 1 ,   λ > 0 . $$ {\mathrm{L}}^{-1}[{s}^{-\beta }\mathrm{exp}(-\lambda {s}^{\alpha })]={t}^{\beta -1}W(-\alpha,\beta;-\lambda {t}^{-\alpha }),\enspace \hspace{1em}0<\alpha < 1,\hspace{1em}\enspace \lambda >0. $$(28)

The Wright function is defined by

W ( α , β ; z ) = n = 0 z n n ! Γ ( β + α n ) ,   α ,   α > - 1 ,   z     C . $$ W\left(\alpha,\beta;z\right)=\sum_{n=0}^{\mathrm{\infty }} \frac{{z}^n}{n!\mathrm{\Gamma }(\beta +{\alpha n})},\enspace \hspace{1em}\alpha,\enspace \alpha >-1,\enspace \hspace{1em}z\enspace \in \enspace C. $$(29)

2.3.2 Probability density function

If f(s) = 1, the right-hand side of equation (23) is a probability density function. Denoting G ( | x | , t ) $ \mathcal{G}(|{x}|,t)$ to be the probability density function, we require that

G ( | x | , t ) > 0 for   all   x , $$ \mathcal{G}\left(\left|{x}\right|,t\right)>0\hspace{1em}\mathrm{for}\enspace \mathrm{all}\enspace x, $$(30)and

0 d x   G ( | x | , t ) = 1 2 . $$ {\int }_0^{\mathrm{\infty }} \mathrm{d}x\enspace \mathcal{G}(|{x}|,t)=\frac{1}{2}. $$(31)

To get the integral of the right-hand side of equation (23), we use

0 z d x   x β ̃ - 1 E α ̃ , β ̃ ( w x α ̃ ) = z β ̃ E α ̃ , β ̃ + 1 ( w z α ̃ ) ,   β > 0 . $$ {\int }_0^z \mathrm{d}x\enspace {x}^{\stackrel{\tilde }{\beta }-1}{E}_{\stackrel{\tilde }{\alpha },\stackrel{\tilde }{\beta }}\left(w{x}^{\stackrel{\tilde }{\alpha }}\right)={z}^{\stackrel{\tilde }{\beta }}{E}_{\stackrel{\tilde }{\alpha },\stackrel{\tilde }{\beta }+1}\left(w{z}^{\stackrel{\tilde }{\alpha }}\right),\hspace{1em}\enspace \beta >0. $$(32)and thus on integration of equation (23) the result is:

0 x d x   γ ̅ = 1 2 η ̃ 1 s 1 - α [ 1 ( s α / η ̃ ) β / ( β + 1 ) x E β + 1,2 ( s α η ̃ x β + 1 ) - x β + 1 E β + 1 , β + 2 ( s α η ̃ x β + 1 ) ] . $$ \begin{array}{ll}{\int }_0^x \mathrm{d}{x}^\mathrm{\prime}\enspace \overline{\gamma }& =\frac{1}{2\stackrel{\tilde }{\eta }}\frac{1}{{s}^{1-\alpha }}\left[\frac{1}{({s}^{\alpha }/\stackrel{\tilde }{\eta }{)}^{\beta /(\beta +1)}}x{E}_{\beta +\mathrm{1,2}}\left(\frac{{s}^{\alpha }}{\stackrel{\tilde }{\eta }}{x}^{\beta +1}\right)\right.\\ & -\left.{x}^{\beta +1}{E}_{\beta +1,\beta +2}\left(\frac{{s}^{\alpha }}{\stackrel{\tilde }{\eta }}{x}^{\beta +1}\right)\right].\end{array} $$(33)

On using the asymptotic expansion in Haubold et al. (2011),

E α ̃ , β ̃ ( z ) = 1 α ̃ z 1 - β ̃ α ̃ exp ( z 1 α ̃ ) - k = 1 z - k Γ ( β - ) , | z | ,   | arg   z | ξ , 0 < α ̃ < 2   and   π α ̃ / 2 < ξ < min   [ π , π α ̃ ] , $$ \begin{array}{c}{E}_{\stackrel{\tilde }{\alpha },\stackrel{\tilde }{\beta }}(z)=\frac{1}{\stackrel{\tilde }{\alpha }}{z}^{\frac{1-\stackrel{\tilde }{\beta }}{\stackrel{\tilde }{\alpha }}}\mathrm{exp}({z}^{\frac{1}{\stackrel{\tilde }{\alpha }}})-\sum_{k=1}^{\mathrm{\infty }} \frac{{z}^{-k}}{\mathrm{\Gamma }(\beta -{k\alpha })},\\ |z|\to \mathrm{\infty },\enspace |\mathrm{arg}\enspace z|\le \xi,0<\stackrel{\tilde }{\alpha }<2\enspace \mathrm{and}\enspace \pi \stackrel{\tilde }{\alpha }/2<\xi < \mathrm{min}\enspace [\pi,\pi \stackrel{\tilde }{\alpha }],\end{array} $$(34)for z ϵ C and ξ is a real number with β < 1, we have

0 x d x   γ ̅ = 1 2 s . $$ {\int }_0^x \mathrm{d}{x}^\mathrm{\prime}\enspace \overline{\gamma }=\frac{1}{2s}. $$(35)

Therefore, equation (31) is satisfied. Consequently, the right-hand side of equation (23) is a pdf.

If both αf and βf equal to 1, we have the Gaussian as the inversion of the right-hand side of γ ̅ ( x ) $ \overline{\gamma }(x)$ in equation (23) yields the expression for γ(x, t) to be

γ ( x , t ) = 1 4 π η ̃ t exp ( - x 2 4 η ̃ t ) . $$ \gamma (x,t)=\frac{1}{\sqrt{4\pi \stackrel{\tilde }{\eta }t}}\mathrm{exp}\left(-\frac{{x}^2}{4\stackrel{\tilde }{\eta }t}\right). $$(36)

The moments of an even order which would also be of interest are (Uchaikin, 2013). For β = 1,

R d | x | 2 n G ( | x | , t ) d x = 2 2 n Γ ( n + d 2 ) Γ ( n + 1 ) Γ ( d 2 ) Γ ( α n + 1 ) η t ,   n = 0,1 , 2 , . . . . $$ {\int }_{{\mathcal{R}}^d} |{x}{|}^{2n}\mathcal{G}(|{x}|,t)\mathrm{d}{x}={2}^{2n}\frac{\mathrm{\Gamma }\left(n+\frac{d}{2}\right)\mathrm{\Gamma }\left(n+1\right)}{\mathrm{\Gamma }\left(\frac{d}{2}\right)\mathrm{\Gamma }\left({\alpha n}+1\right)}\eta {t}^{{n\alpha }},\enspace \hspace{1em}n=\mathrm{0,1},2,.... $$(37)

2.4 Pressure distributions, analogs

Using the expressions given above, we consider the pressure distribution for a fractured well. Following Ozkan and Raghavan (1991), the pressure distribution, Δp( x t), in terms of the Laplace transformation for a continuous, plane/line source at location, M′, of strength q ̃ ( M , t ) / ϕ c $ \mathop{q}\limits^\tilde(M\mathrm{\prime},t)/{\phi c}$ is

Δ p ̅ f ( x ) = μ k α , β S q ̃ ̅ ( x ) γ ̅ f ( x - x ) d S , $$ {\overline{\Delta p}}_f({x})=\frac{\mu }{{k}_{\alpha,\beta }}{\int }_S \overline{\mathop{q}\limits^\tilde}({x}\mathrm{\prime}){\overline{\gamma }}_f({x}-{x}\mathrm{\prime})\mathrm{d}S\mathrm{\prime}, $$(38)with the symbol S in the case under study representing the surface of the source. If q ̃ ( M , t ) $ \mathop{q}\limits^\tilde(M\mathrm{\prime},t)$ were constant, then equation (38) simplifies to

Δ p ̅ f ( x ) = μ k α , β q ̃ ( x ' ) s S γ ̅ f ( x - x ) d S . $$ {\overline{\Delta p}}_f({x})=\frac{\mu }{{k}_{\alpha,\beta }}\frac{\mathop{q}\limits^\tilde({{x}}^{\prime})}{s}{\int }_S {\overline{\gamma }}_f({x}-{{x}}^\mathrm{\prime})\mathrm{d}{S}^\mathrm{\prime}. $$(39)

For the system under consideration, equations (23) and (38) yield the following expression for pressure distribution:

Δ p ̅ f ( y ̃ , s ) = Δ p ̅ f   ( 0 , s ) E β f + 1 ( λ ̃ y ̃ β f + 1 ) + [ d β f d y ̃ β f   Δ p ̅ f ( y ̃ , s ) ] y ̃ = 0 y ̃ β f E β f + 1 , β f + 1 ( λ ̃ y ̃ β f + 1 ) , $$ \begin{array}{ll}{\overline{\Delta p}}_f(\mathop{y}\limits^\tilde,s)& ={\overline{\Delta p}}_f\enspace (0,s){E}_{{\beta }_f+1}(\stackrel{\tilde }{\lambda }{\mathop{y}\limits^\tilde}^{{\beta }_f+1})\\ & +{\left[\frac{{\mathrm{d}}^{{\beta }_f}}{\mathrm{d}{\mathop{y}\limits^\tilde}^{{\beta }_f}}\enspace {\overline{\Delta p}}_f(\mathop{y}\limits^\tilde,s)\right]}_{\mathop{y}\limits^\tilde=0}{\mathop{y}\limits^\tilde}^{{\beta }_f}{E}_{{\beta }_f+1,{\beta }_f+1}(\stackrel{\tilde }{\lambda }{\mathop{y}\limits^\tilde}^{{\beta }_f+1}),\end{array} $$(40)where y ̃ 0 $ \mathop{y}\limits^\tilde\ge 0$ is the distance from the source. The flux distribution in the system corresponding to equation (40) is

[ d β f d y ̃ β f   Δ p ̅ f   ( y ̃ , s ) ] = λ ̃ Δ p ̅ f   ( 0 , s ) E β f + 1,2 ( λ ̃ y ̃ β f + 1 ) y ̃ + [ d β f d y ̃ β f   Δ p ̅ f   ( y ̃ , s ) ] y ̃ = 0 E β + 1 ( λ ̃ y ̃ β f + 1 ) , $$ \begin{array}{ll}\left[\frac{{\mathrm{d}}^{{\beta }_f}}{\mathrm{d}{\mathop{y}\limits^\tilde}^{{\beta }_f}}\enspace {\overline{\Delta p}}_f\enspace (\mathop{y}\limits^\tilde,s)\right]& =\stackrel{\tilde }{\lambda }{\overline{\Delta p}}_f\enspace (0,s){E}_{{\beta }_f+\mathrm{1,2}}\left(\stackrel{\tilde }{\lambda }{\mathop{y}\limits^\tilde}^{{\beta }_f+1}\right)\mathop{y}\limits^\tilde\\ & +{\left[\frac{{\mathrm{d}}^{{\beta }_f}}{\mathrm{d}{\mathop{y}\limits^\tilde}^{{\beta }_f}}\enspace {\overline{\Delta p}}_f\enspace (\mathop{y}\limits^\tilde,s)\right]}_{\mathop{y}\limits^\tilde=0}{E}_{\beta +1}\left(\stackrel{\tilde }{\lambda }{\mathop{y}\limits^\tilde}^{{\beta }_f+1}\right),\end{array} $$(41)with

Δ p ̅ f   ( 0 , s ) + 1 λ ̃ β f / ( β f + 1 ) [ d β f d y ̃ β f   Δ p ̅ f   ( y ̃ , s ) ] y ̃ = 0 = 0 . $$ {\overline{\Delta p}}_f\enspace (0,s)+\frac{1}{{\stackrel{\tilde }{\lambda }}^{{\beta }_f/({\beta }_f+1)}}{\left[\frac{{\mathrm{d}}^{{\beta }_f}}{\mathrm{d}{\mathop{y}\limits^\tilde}^{{\beta }_f}}\enspace {\overline{\Delta p}}_f\enspace (\mathop{y}\limits^\tilde,s)\right]}_{\mathop{y}\limits^\tilde=0}=0. $$(42)

The above expressions are also given in Chen and Raghavan (2015).

For a bounded system, the pressure distribution, Δ p f ( y ̃ , t ) $ \Delta {p}_f(\mathop{y}\limits^\tilde,t)$, in terms of the Laplace transformation may be obtained from

Δ p ̅ f ( y ̃ , s ) = Δ p ̅ f   ( 0 , s ) { E β f + 1 ( λ ̃ y ̃ β f + 1 ) - φ y ̃ β f E β f + 1 , β f + 1 ( λ ̃ y ̃ β f + 1 ) } ,   for   y ̃ y ̃ e , $$ \begin{array}{ll}{\overline{\Delta p}}_f(\mathop{y}\limits^\tilde,s)& ={\overline{\Delta p}}_f\enspace (0,s)\left\{{E}_{{\beta }_f+1}(\stackrel{\tilde }{\lambda }{\mathop{y}\limits^\tilde}^{{\beta }_f+1})\right.\\ & -\left.\phi {\mathop{y}\limits^\tilde}^{{\beta }_f}{E}_{{\beta }_f+1,{\beta }_f+1}(\stackrel{\tilde }{\lambda }{\mathop{y}\limits^\tilde}^{{\beta }_f+1})\right\},\enspace \hspace{1em}\mathrm{for}\enspace \mathop{y}\limits^\tilde\le {\mathop{y}\limits^\tilde}_e,\end{array} $$(43)where

φ = λ ̃ y ̃ e E β f + 1,2 ( λ ̃ y ̃ e β f + 1 ) E β f + 1 ( λ ̃ y ̃ e β f + 1 ) . $$ \phi =\stackrel{\tilde }{\lambda }{\mathop{y}\limits^\tilde}_e\frac{{E}_{{\beta }_f+\mathrm{1,2}}(\stackrel{\tilde }{\lambda }{\mathop{y}\limits^\tilde}_e^{{\beta }_f+1})}{{E}_{{\beta }_f+1}(\stackrel{\tilde }{\lambda }{\mathop{y}\limits^\tilde}_e^{{\beta }_f+1})}. $$(44)

We may relate the fractional gradient to the pressure drop at y ̃ = 0 $ \mathop{y}\limits^\tilde=0$ by the expression

[ d β f d y ̃ β f   Δ p ̅ f ( y ̃ , s ) ] y ̃ = 0 = - φ Δ p ̅ f ( y ̃ , s ) | y ̃ = 0 . $$ {\left[\frac{{\mathrm{d}}^{{\beta }_f}}{\mathrm{d}{\mathop{y}\limits^\tilde}^{{\beta }_f}}\enspace {\overline{\Delta p}}_f(\mathop{y}\limits^\tilde,s)\right]}_{\mathop{y}\limits^\tilde=0}=-\phi {\overline{\Delta p}}_f(\mathop{y}\limits^\tilde,s){|}_{\mathop{y}\limits^\tilde=0}. $$(45)

Equation (43) is an approximation and a general expression for the pressure distribution at any point in a composite system is available; see Raghavan and Chen (2018a). Equations (40) and (43) may be coupled with the expression for the pressure distribution within the fracture to obtain well responses; equation (22).

2.5 The multiple fracture scheme

Given the pressure response at a well draining a single fracture such as equation (43), a scheme involving the draining of a system through multiple fractures intersecting a wellbore may be addressed by Duhamel’s theorem. The scheme described below follows Spath et al. (1994) and may incorporate many kinds of variations in rock and fracture properties, such as hydraulic fracture conductivity and hydraulic fracture length. Consider a system which is produced through N fractures. On noting that the total production rate, q, at any instant in time, is given by the summation of the individual flow rates qj at each Fracture j, we may write

j = 1 N q j ( t ) = q ( t ) . $$ \sum_{j=1}^N {q}_j(t)=q(t). $$(46)

The assumption of an infinite-conductivity wellbore requires, of course, that the same pressure be maintained at every point inside the horizontal wellbore, namely, that for each fracture location, j = 1, 2, …, N we have

Δ p w ( t ) = Δ p w j ( t ) . $$ \Delta {p}_w(t)=\Delta {p}_{{w}_j}(t). $$(47)

The convolution equation to compute the pressure response at the wellbore for the multiple-fracture system is

Δ p w ( t ) = j = 1 n 0 t q j ( τ ) Δ p w jk ( t - τ ) d τ , $$ \Delta {p}_w(t)=\sum_{j=1}^n \int^t_0 {q}_j(\tau )\Delta p\mathrm{\prime}_{{w}_{{jk}}}(t-\tau )\mathrm{d}\tau, $$(48)for k = 1, 2, …, N. Here, jk represents the effect of Fracture j on the pressure response at Fracture k and Δ p wjk ( t ) = d Δ p wjk / d t $ \Delta p\mathrm{\prime}_{{wjk}}(t)=\mathrm{d}\Delta {p}_{{wjk}}/\mathrm{d}t$.

In terms of the Laplace transformation, the above three expressions, respectively, are as follows:

j = 1 N q ¯ j = q ¯ , $$ \sum_{j=1}^N {\bar{q}}_j=\bar{q}, $$(49)

Δ p ̅ w = Δ p ̅ w j = 1 , n , $$ {\overline{\Delta p}}_w={\overline{\Delta p}}_{{w}_{j=1,n}}, $$(50)and

Δ p ̅ w = j = 1 n s q ¯ j Δ p ̅ w jj . $$ {\overline{\Delta p}}_w=\sum_{j=1}^n s{\bar{q}}_j{\overline{\Delta p}}_{{w}_{{jj}}}. $$(51)

Thus we are required to solve simultaneously for Δpw and qj; j = 1, 2, …, N, namely the system

[ s Δ p ̅ w 11 - 1 s Δ p ̅ w 21 - 1 s Δ p ̅ w pq - 1 s Δ p ̅ w nn - 1 1 1 1 1 0 ] [ q ¯ 1 q ¯ 2 q ¯ n Δ p ̅ w ] = [ 0 0 0 q ¯ ] . $$ \left[\begin{array}{lllll}s{\overline{\Delta p}}_{{w}_{11}}& \cdots & \cdots & \cdots & -1\\ \cdots & s{\overline{\Delta p}}_{{w}_{21}}& \cdots & \cdots & -1\\ \cdots & \cdots & s{\overline{\Delta p}}_{{w}_{{pq}}}& \cdots & -1\\ \cdots & \cdots & \cdots & s{\overline{\Delta p}}_{{w}_{{nn}}}& -1\\ 1& 1& 1& 1& 0\end{array}\right]\left[\begin{array}{l}{\bar{q}}_1\\ {\bar{q}}_2\\ \cdots \\ {\bar{q}}_n\\ {\overline{\Delta p}}_w\end{array}\right]=\left[\begin{array}{l}0\\ 0\\ \cdots \\ 0\\ \bar{q}\\ \end{array}\right]. $$(52)

In general a formulation of this type would in practice best represent situations at early times (Raghavan et al., 1997); nonetheless, formulations of this kind are often used to represent extensions of the trilinear version of the model. The off-diagonal elements of the matrix would have reflected interactions among fractures. Each Δ p w ii $ \Delta {p}_{{w}_{{ii}}}$ element may be specified independent of other such elements; consequently, different fracture properties may be incorporated at each location.

3 Asymptotic solutions

Regardless of the method of evaluating measurements, asymptotic solutions are helpful, even essential, to meet or satisfy a number of objectives. For example, though approximate they provide the following: qualitative information as to whether initial expectations were appropriate, guidance on the appropriate physics that needs to be incorporated in evaluating the response, guidance on the type of the mathematical model to be used in quantitative evaluations, planning of future tests and the like. For example, if we consider a fractured well producing a fractured reservoir, we find that the characteristic responses expected of fractured wells, such as a slope of 1/4 on log-log coordinates during bilinear flow in the classical case, are not evident if transient effects in the matrix significantly influence the well response; see Yeh et al. (1986).

Expressions for responses at a horizontal well producing through multiple fractures that are yet to be presented in the literature are derived in the Appendix for the model considered here. Expressions are derived for all Flow Regimes considered. Table 2 is a summary of the exponents derived in the Appendix; we document results in terms of αf, βf and αm. Specifics of the exponents presented in Table 2 are discussed in the Results section. The exponents given do revert to the classical cases in the literature under identical conditions either for fissured or non-fissured idealizations. In fact, Table 2 reproduces all information available in the literature.

Table 2

Exponents of flow regimes.

4 Connection to subdiffusive solutions

To obtain a measure of the accuracy of the expression derived in equation (43), it is perhaps worthwhile to compare this solution with the solution for subdiffusive flow as it is possible to derive solutions from first principles through standard operational techniques such as through the method of images, Fourier series, etc. On recasting equation (43) and equation (44), we obtain

Δ p ̅ f ( y ̃ , s ) = Δ p ̅ f   ( 0 , s ) E β f + 1 ( λ ̃ y ̃ e β f + 1 ) [ E β f + 1 ( λ ̃ y ̃ β f + 1 ) E β f + 1 ( λ ̃ y ̃ e β f + 1 ) - λ ̃ y ̃ β f y ̃ e E β f + 1,2 ( λ ̃ y ̃ e β f + 1 ) E β f + 1 , β f + 1 ( λ ̃ y ̃ β f + 1 ) ] . $$ \begin{array}{ll}{\overline{\Delta p}}_f(\mathop{y}\limits^\tilde,s)& =\frac{{\overline{\Delta p}}_f\enspace (0,s)}{{E}_{{\beta }_f+1}(\stackrel{\tilde }{\lambda }{\mathop{y}\limits^\tilde}_e^{{\beta }_f+1})}\left[{E}_{{\beta }_f+1}(\stackrel{\tilde }{\lambda }{\mathop{y}\limits^\tilde}^{{\beta }_f+1}){E}_{{\beta }_f+1}(\stackrel{\tilde }{\lambda }{\mathop{y}\limits^\tilde}_e^{{\beta }_f+1})\right.\\ & -\left.\stackrel{\tilde }{\lambda }{\mathop{y}\limits^\tilde}^{{\beta }_f}{\mathop{y}\limits^\tilde}_e{E}_{{\beta }_f+\mathrm{1,2}}(\stackrel{\tilde }{\lambda }{\mathop{y}\limits^\tilde}_e^{{\beta }_f+1}){E}_{{\beta }_f+1,{\beta }_f+1}(\stackrel{\tilde }{\lambda }{\mathop{y}\limits^\tilde}^{{\beta }_f+1})\right].\end{array} $$(53)

For subdiffusive flow; that is, for βf = 1, equation (53) yields the expression

Δ p ̅ f ( y ̃ , s ) = Δ p ̅ f ( 0 , s ) cosh [ λ ̃ ( y ̃ e - y ̃ ) ] cosh ( λ ̃ y ̃ e ) , $$ {\overline{\Delta p}}_f(\mathop{y}\limits^\tilde,s)={\overline{\Delta p}}_f(0,s)\frac{\mathrm{cosh}[\sqrt{\stackrel{\tilde }{\lambda }}({\mathop{y}\limits^\tilde}_e-\mathop{y}\limits^\tilde)]}{\mathrm{cosh}(\sqrt{\stackrel{\tilde }{\lambda }}{\mathop{y}\limits^\tilde}_e)}, $$(54)which is identical to what we would obtain through standard operational techniques for subdiffusive flow.

It is noteworthy and compelling that equation (53) yields an expression identical to the one in equation (54) for βf = 1; we find it to be so because equation (53) is an approximation obtained through the solution for a composite reservoir which in itself was derived because options such methods as the method of images are not applicable for superdiffusive flow (Metzler et al., 2007). For us, the correspondence shown above is of more than academic interest as it provides a glimpse of the structure of the solution for superdiffusion, particularly in 2D for those who wish to pursue analytical options; it is clear that deep relations that are yet to be explored, expressions of the form E β f + 1,2 ( z 1 ) / E β f + 1 , β f + 1 ( z 2 ) $ {E}_{{\beta }_f+\mathrm{1,2}}\left({z}_1\right)/{E}_{{\beta }_f+1,{\beta }_f+1}\left({z}_2\right)$ would exist for situations wherein the fracture length and the width of the rectangle are not identical as indicated in Figure 4 (xe/xf > 1); see Appendix.

thumbnail Fig. 4

A fractured well in a rectangular drainage region. The fracture with identical wing lengths with its center at (xw, yw) is located anywhere in the drainage region and is parallel to one of the boundaries. The length, width and height of the fracture are xF, wF, and h, respectively.

In deriving an asymptotic solution for Flow Regimes 4 and 5, we show in the Appendix that two-term approximations for Eα,β(z) as z → 0 provide the needed expression for pressure distributions rather than the indirect approach used in Raghavan and Chen (2018a). We also illustrate that this approximation reduces to that obtained from the analog of equation (54) as s → 0 through the results in (1.445,2) and (1.445,2) of Gradshteyn and Ryzhik (1965).

5 Results

Thus far we have summarized the outline needed to conduct simulations. In this section we document actual results based on our simulations. We first discussed the results regarding exponents summarized in Table 2 and explain factors that govern the magnitude of the exponents in terms of its constituents, αf, βf and αm. We then document results of simulations for a few cases and show that the exponents predicted in Table 2 are, indeed, reproduced in the simulations. Third we expand upon an actual analysis, determine ranges for the three components of the exponents, point out the consequences of these estimates by determining time-distance relationships for the transient front and show they are much different from expectations for classical diffusion. We also discuss the purport of the ranges of the exponents in practical terms such as well spacing. All results are computed through the Stehfest Algorithm (1970a, b) and the Mittag-Leffler functions are calculated in the manner indicated in Gorenflo et al. (2002).

5.1 Exponents of flow regimes

The following comments on the exponent during the transient period, though obvious, are important. Although dependent on the specific values of αf, βf and αm some of the observations for classical diffusion do hold; for example, on considering Flow Regimes 1 and 3 the observation that parallel straight lines should be evident is valid for anomalous diffusion, and the ratio of the slopes (factor of 2 or doubling of slopes) between bilinear and linear Flow Regimes is preserved as in the case of classical diffusion. If we consider Flow Regime 3 our expectations of the roles of αf and βf are met. A reduction in βf (more communication) results in a reduction in the value of the exponent; however, as expected the opposite happens if there were a reduction in αf (greater restriction). Inspection of the role of the exponents for Flow Regime 2, however, suggests a more complicated picture: For fixed values of αf and βf, the slope increases as αm decreases. But changes in the slope are more complicated for fixed values of αf as βf changes, and depends on the magnitude of αm. If αm were <1, slopes changes depend on the values of αm and αf for fixed values of βf; then slope changes depend on the values of αm, αf for fixed values of βf; however, if αm = 1, then for fixed values of αf the slope decreases as βf decreases. The product αm βf in the numerator causes the differences just noted.

If β were < 1, then all other things being equal the exponent (slope of the line on log-log coordinates) will be even much smaller. Chu (2018), who framed his comments in terms of the Chow (1952) Pressure Group, Δpc, has noted that we need to consider situations when slopes are very, very small implying situations beyond subdiffusive effects. And super-connected networks resulting from βf < 1 may result in slopes (exponents of t) of lines on log-log plots being much less than the slopes based on classical expectations given in Yeh et al. (1986); that is, being < 0.125.

While on this subject, we should note that if the values of αf and βf were such that β f + 1 - 2 α f = 0 $ {\beta }_f+1-2{\alpha }_f=0$, then the exponents for the linear flow period, Flow Regimes 1 or 3, noted in Table 2 suggest that the well response is identical to the classical case; that is, the slope would be 1/2. Similarly, if β f + 1 - 2 α f - α m β f $ {\beta }_f+1-2{\alpha }_f-{\alpha }_m{\beta }_f$ were equal to 0, then the models of anomalous diffusion for Flow Regime 2 in Table 2 suggest that the trend in behaviors for anomalous diffusion would be similar to those predicted by classical models. Essentially, many points on the (αf − βf) or (αf − βf – αm) phase planes yield slopes identical to those of classical diffusion because the points representing classical diffusion, namely, αf of 1, βf of 1, and αm of 1, are simply coordinates on the phase plane in the context of anomalous diffusion.

Finally, we should note that the exponents given here for Flow Regimes 4 and 5 also hold for 2D systems.

5.2 Computational results

Figures 5 and 6 present sample calculations of responses at a well produced at a constant rate through a horizontal well intercepted by multiple fractures. We presume the well drains the system through 300 fractures. For the discussion of these two figures, we posit that fracture properties are identical. Dimensionless pressure, pwD(tD), dimensionless time, tD, and dimensionless fracture conductivity, CFD, are defined, respectively, as

p wD ( t D ) = 2 π k α f , β f h t η ̃ f 1 - α f α f L F β f + 1 - 2 α f α f Δ p w ( t ) , $$ {p}_{{wD}}\left({t}_D\right)=\frac{2\pi {k}_{{\alpha }_f,{\beta }_f}{h}_t}{{q\mu }}\frac{{\stackrel{\tilde }{\eta }}_f^{\frac{1-{\alpha }_f}{{\alpha }_f}}}{{L}_F^{\frac{{\beta }_f+1-2{\alpha }_f}{{\alpha }_f}}}\Delta {p}_w(t), $$(55)

t D = η ̃ L F β f + 1 t α f , $$ {t}_D=\frac{\stackrel{\tilde }{\eta }}{{L}_F^{{\beta }_f+1}}{t}^{{\alpha }_f}, $$(56)and

C FD = k F w F k α f , β f L F L F β f + 1 - 2 α f α f η ̃ f 1 - α f α f . $$ {C}_{{FD}}=\frac{{k}_F{w}_F}{{k}_{{\alpha }_f,{\beta }_f}{L}_F}\frac{{L}_F^{\frac{{\beta }_f+1-2{\alpha }_f}{{\alpha }_f}}}{{\stackrel{\tilde }{\eta }}_f^{\frac{1-{\alpha }_f}{{\alpha }_f}}}. $$(57)

thumbnail Fig. 5

Production response at a horizontal well produced through multiple fractures at a constant rate to highlight Flow Regimes 1 and 3; the fracture conductivity is assumed to be infinitely large, CFD ≥ 500. The unbroken lines are the pressure, pwD, and derivative, p wD $ {p}_{{wD}}^\mathrm{\prime}$, responses respectively. Asymptotic solutions corresponding to short and long times, Flow Regime 1 and Flow Regime 3, respectively, are shown as dashed lines. The circles are solutions for a single fracture presented in the manner indicated in equation (58).

thumbnail Fig. 6

Production response at a horizontal well produced through multiple fractures at a constant rate to highlight Flow Regime 2; two values of fracture conductivity are considered, CFD of 1 and 5. The unbroken lines are the pressure, pwD, and derivative, p wD $ {p}_{{wD}}^\mathrm{\prime}$, responses respectively. The dashed lines reflect the asymptotic solution, Flow Regime 2. Values of λ′ and ω′ are so chosen to facilitate fluid transfer and promote the duration of Flow Regime 2. The circles are solutions for a single fracture presented in the manner indicated by equation (58).

These definitions follow from the more general definitions discussed in the Appendix.

Figure 5 highlights responses during Flow Regimes 1 and 3; specific properties used to calculate dimensionless pressure, pwD(tD), and its derivative p' wD ( t D ) $ {{p\prime}_{{wD}}({t}_D)$, (t|dpw/dt|) are noted on Figure 5. The dashed lines reflect the asymptotic solutions derived in the Appendix for Flow Regimes 1 and 3; the exponents noted in Table 2 are reproduced. The circles highlighted by the abbreviation sf represent calculated responses for a well produced through a single fracture and are correlated in the manner shown because the expression

p wD ( t D ) | sf = N p wD ( t D ) , $$ {{p}_{{wD}}\left({t}_D\right)|}_{{sf}}=N{p}_{{wD}}({t}_D), $$(58)where N is the number of fractures, holds during times when there is no interaction between the fractures (Raghavan et al., 1997). This kind of correlation is possible even if fracture properties are not identical.

The intention of Figure 6 is to demonstrate similar information for Flow Regime 2 and two values of CFD, namely 1 and 5 are considered. The chosen values of λ′ and ω′ (in comparison to those used in Fig. 5) enhance rapid fluid transfer from the matrix system and also extend the duration of the transient period in the matrix. The dashed lines represent responses during the asymptotic solution derived in the Appendix, Flow Regime 2 and the exponent given in Table 2. Finally, the computations of the single-fracture solution identified by the circles in Figure 6 again confirm that the single and multiple fracture solutions correlate in the manner indicated by equation (58).

5.3 Pressure tests: α and β

We noted in the Introduction that our intent is to consider patterns similar to those in Doe et al. (2014) shown in Figure 1. Figure 7 taken from Raghavan and Chen (2019) is one such case, and offers a plot of flowing and buildup responses in the standard way; that is, pressure and the corresponding derivative at an oil well with the filled-in circles being buildup responses (Δp vs. Δt) and the squares representing flowing responses IRNP(te) vs. te, in terms of the nomenclature used in Rate-Transient-Analysis, RTA. The symbol IRNP(te) is the integral of rate-normalized pressure defined in Palacio and Blasingame (1993) as

IRNP ( t e ) = 1 t e 0 t e p i - p wf ( t ) q ( t )   d t , $$ \mathrm{IRNP}\left({t}_e\right)=\frac{1}{{t}_e}{\int }_0^{{t}_e} \frac{{p}_i-{p}_{{wf}}\left({t}^\mathrm{\prime}\right)}{q\left({t}^\mathrm{\prime}\right)}\enspace \mathrm{d}{t}^\mathrm{\prime}, $$(59)and the symbol te represents material-balance time defined in Camacho-Velázquez and Raghavan (1989) as

t e = 1 q ( t ) 0 t q ( t ) d t . $$ {t}_e=\frac{1}{q(t)}{\int }_0^t q\left({t}^\mathrm{\prime}\right)\mathrm{d}{t}^\mathrm{\prime}. $$(60)

thumbnail Fig. 7

Buildup (circles) and flowing well (squares) responses for a shale, oil well.

(In this example, no other information other than that shown here is available to us as we were asked to show a proof of the concept.) Here we consider issues beyond those identified in Raghavan and Chen (2019) who showed that the flowing and buildup responses correspond to Flow Regime 4 and Flow Regime 2, respectively. The alignment of the flowing and buildup responses is striking and establishes the need to buildup data which is not often measured although early-time flowing-well responses are often quite suspect. The basis for combining buildup and responses is discussed in Raghavan (1980). The slope at late times corresponding to Flow Regime 4 is much larger than what we would expect for classical diffusion and based on the exponent in Table 2 suggests that the value of αm is 0.6. By the same token the slope of the first straight line in Figure 7, namely, 0.1 reflective of Flow Regime 2 is much smaller than expectations based on classical diffusion, 1/4 or 1/8. The estimates provide a well-founded basis for a model of the kind proffered here. The value of αm along with the estimate of 0.1 for the early-time slope may be translated to a range of values for αf and βf. Scouring a map of possible ranges of αf and βf, we arrive at 0.8 ≤ αf ≤ 1 and 0.1 ≤ βf ≤ 0.4. This analysis suggests that long-range interactions are important to the performance of this well. Consequently, it behooves the operator to estimate values of αf and βf at other wells draining the same formation. If that were the case, it is incumbent to evaluate whether the origins lie in the mechanical properties of the porous rock or in the geological environment. Once, this is determined, it seems sensible to obtain geological concurrence and advisable to make modifications such as fracture spacing, well spacing, etc. Such an exercise may have significant economic benefit.

Figure 8 presents distributions of the Green’s function for three situations based on the interpretation in Figure 7. The top curve representing purely subdiffusive flow reflects Flow Regime 4 and the two bottom curves include the combined effects of sub- and superdiffusion subject to the constraint, namely, that the exponent, m, has the value of 0.1, where

m = 2 ( β f + 1 ) - 2 α f - α m β f 4 ( β f + 1 ) , $$ m=\frac{2({\beta }_f+1)-2{\alpha }_f-{\alpha }_m{\beta }_f}{4({\beta }_f+1)}, $$(61)with αm = 0.6. Differences between the bottom two curves are to only be expected; all that is indicated is the well-known point that many combinations of properties yield the same pressure response.

thumbnail Fig. 8

Green’s function responses corresponding to the test discussed in Figure 7; subdiffusive responses correspond to Flow Regime 4 and combined responses to Flow Regime 2.

All curves shown here, and others, of course, including the Gaussian may be expressed through the similarity variable, ς,

ς = ( t D x D β + 1 ) 1 β + 1 , $$ \varsigma ={\left(\frac{{t}_D}{{x}_D^{\beta +1}}\right)}^{\frac{1}{\beta +1}}, $$(62)by the expression

x D G ( x D , t D ) = f ( ς ) , $$ {x}_D\mathcal{G}({x}_D,{t}_D)=f(\varsigma ), $$(63)where ς is deduced by considering the arguments of the Mittag-Leffler, Fox and Wright functions on the right-hand sides of equations (24), (27), and (29), respectively, and the solutions for Flow Regime 1 in the Appendix; see equation (A.36). The ability to correlate solutions in terms of ς is shown in Figure 9 where a typical plot of x D G ( x D , t D ) $ {x}_D\mathcal{G}({x}_D,{t}_D)$ as function of ς we have constructed is shown for several values of xD and tD. The circles, squares, diamonds, etc. represent specific vales of xD as noted in Figure 9. Curves for distinct values of locations, xD, may be correlated in this manner for specific combinations of α and β.

thumbnail Fig. 9

Correlation of x D G ( x D , t D ) $ {x}_D\mathcal{G}({x}_D,{t}_D)$ as a function of the similarity variable, ς, of plane source solutions. Values of xD and tD in the range 0.2 ≤ xD ≤ 5 and 3 × 10−2 ≤ tD ≤ 500 are considered.

For α ≤ 1 and β = 1, f(ς) in equation (63) may be expressed analytically through the Wright function, W(α, β; t), as was noted in equation (28); thus

x D G ( x D , t D ) = 1 2 ς W ( - α / 2,1 - α / 2 ; - 1 / ς ) , $$ {x}_D\mathcal{G}({x}_D,{t}_D)=\frac{1}{2\varsigma }W(-\alpha /\mathrm{2,1}-\alpha /2;-1/\varsigma ), $$(64)noting that

W ( - 1 / 2,1 / 2 ; - z ) = 1 π exp ( - z 2 4 ) . $$ W(-1/\mathrm{2,1}/2;-z)=\frac{1}{\sqrt{\pi }}\mathrm{exp}\left(-\frac{{z}^2}{4}\right). $$(65)

6 Discussion and concluding remarks

This work attempts to identify and quantify patterns in the spatiotemporal evolution in transient behavior wherein the geology is complex, where the interconnected structure is composed of multiple links and microscale discontinuities at the Theis (1935) scale. A motivation for using fractional flux laws is presented. As such patterns are repeatable at many levels both in individual wells and across the resource, it is imperative that new methods of analysis be available. Patterns of derivative responses shown in Figure 1 are a typical example that may not be evaluated through classical methods; in shale systems, in most cases, accepted relationships based on classical norms demonstrated in the literature, even when officially sanctioned, are no more than speculation. Exponents displayed in the curves in Figure 1 convey important information, particularly regarding scaling behaviors and no longer need to be ignored. Results reported here may be used for further analysis and/or modelling at any scale (equivalent porous media models, DFN schemes) through the constants α and β; for regardless of the approach used the following factors are common for predicting the performance of fractured rocks: High and low conductivity paths, discontinuous properties that encompass several orders of magnitude, convoluted flow-paths that serve as conduits and barriers and complex geology (see Graphical Abstract). None of these may be addressed in the classical context. One particular disadvantage of DFN models is that, to our knowledge, at this time no particular rules or guidelines of a general kind that are normally available in evaluating transient responses exist. Our work provides an option to overcome this limitation. A recent contribution (Artus, 2020) that came to our attention as this work was being completed gives much credence to a model that includes sub- and superdiffusive options of the kind discussed here.

A.1 Asymptotic solutions, development

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. To obtain characteristic responses at early times, namely Flow Regimes 1, 2 and 3, we proceed under the assumption of 1D flow in an infinite reservoir. Modelling, for convenience, one fourth of the drainage region shown in Figure 4, namely x  ≥  x w  = 0 and assuming all distances are measured from the well center, equations (6), (12) and (13) subject to the assumptions stated below equation (13) lead to

2 Δ p ̅ F x 2 + 2 k α f , β f k F w s 1 - α f β f Δ p ̅ f y ̃ β f | y ̃ = 0 = s η F Δ p ̅ F , $$ \frac{{\mathrm{\partial }}^2{\overline{\Delta p}}_F}{\mathrm{\partial }{x}^2}+2\frac{{k}_{{\alpha }_f,{\beta }_f}}{{k}_Fw}{s}^{1-{\alpha }_f}\frac{{\mathrm{\partial }}^{{\beta }_f}{\overline{\Delta p}}_f}{\mathrm{\partial }{\mathop{y}\limits^\tilde}^{{\beta }_f}}{|}_{\mathop{y}\limits^\tilde=0}=\frac{s}{{\eta }_F}{\overline{\Delta p}}_F, $$(A.1)with

η F = k F ϕ F c F μ . $$ {\eta }_F=\frac{{k}_F}{{\phi }_F{c}_F\mu }. $$(A.2)

For a system produced at a constant rate, q, noting that the well tips are assumed to be sealed, equation (A.1) yields the pressure distribution in the fracture, Δ p ̅ F ( y D , s ) $ {\overline{\Delta p}}_F({y}_D,s)$, to be

Δ p ̅ F ( x , s ) = 2 s q ̃ F k F w F h t cosh [ q ̃ F ( L F - x ) ] sinh q ̃ F L F , $$ {\overline{\Delta p}}_F(x,s)=\frac{{q\mu }}{2s{\mathop{q}\limits^\tilde}_F{k}_F{w}_F{h}_t}\frac{\mathrm{cosh}[{\mathop{q}\limits^\tilde}_F({L}_F-x)]}{\mathrm{sinh}{\mathop{q}\limits^\tilde}_F{L}_F}, $$(A.3)with q ̃ F $ {\mathop{q}\limits^\tilde}_F$ given as

q ̃ F 2 = 2 k α f , β f k F w F s 1 - α f λ ̃ β ̃ f + s η F , $$ {\mathop{q}\limits^\tilde}_F^2=2\frac{{k}_{{\alpha }_f,{\beta }_f}}{{k}_F{w}_F}{s}^{1-{\alpha }_f}{\stackrel{\tilde }{\lambda }}^{{\stackrel{\tilde }{\beta }}_f}+\frac{s}{{\eta }_F}, $$(A.4)where

β ̃ f = β f β f + 1 , $$ {\stackrel{\tilde }{\beta }}_f=\frac{{\beta }_f}{{\beta }_f+1}, $$(A.5)by virtue of equation (42). We thus find that the well pressure depends on a second tanh(x) function, as the well pressure for equation (A.3) reduces to

Δ p ̅ w ( s ) = 2 s h t [ q ̃ F k F w F tanh ( q ̃ F L F ) ] - 1 . $$ {\overline{\Delta p}}_w(s)=\frac{{q\mu }}{2s{h}_t}{\left[{\mathop{q}\limits^\tilde}_F{k}_F{w}_F\mathrm{tanh}({\mathop{q}\limits^\tilde}_F{L}_F)\right]}^{-1}. $$(A.6)

Thus, on defining

p D ( x , t ) = 2 π k α f , β f h t Δ p ( x , t ) , $$ {p}_D^{\star }(x,t)=\frac{2\pi {k}_{{\alpha }_f,{\beta }_f}{h}_t}{{q\mu }}\Delta p(x,t), $$(A.7)and

C FD = k F w F k α f , β f , $$ {C}_{{FD}}^{\star }=\frac{{k}_F{w}_F}{{k}_{{\alpha }_f,{\beta }_f}}, $$(A.8)we obtain the well pressure in terms of the Laplace transformation to be

p ̅ wD ( s ) = π s [ q ̃ F C FD tanh ( q ̃ F L F ) ] - 1 . $$ {\bar{p}}_{{wD}}^{\star }(s)=\frac{\pi }{s}{\left[{\mathop{q}\limits^\tilde}_F{C}_{{FD}}^{\star }\mathrm{tanh}({\mathop{q}\limits^\tilde}_F{L}_F)\right]}^{-1}. $$(A.9)

The above expression may be readily extended to the situation where a horizontal well is produced through multiple fractures because at early times each fracture acts individually; that is, each fracture acts as if the others do not exist, therefore it is possible to amalgamate responses. Converting equation (A.9) to the situation where the fracture produces at a constant pressure, noting that the total rate for constant-pressure-production may thus be obtained by summing each individual response as all fractures produce against a common well pressure and reconverting the resulting sum to the situation where the horizontal well produces at a constant rate, q, we have

p ¯ wD ( s ) = π s [ l = 1 N q ̃ F l C FD l tanh ( q ̃ F l L F l ) ] - 1 , $$ {\bar{p}}_{{wD}}^{\star }(s)=\frac{\pi }{s}{\left[\sum_{\mathcal{l}=1}^N {\mathop{q}\limits^\tilde}_{F\mathcal{l}}{C}_{{FD}\mathcal{l}}^{\star }\mathrm{tanh}({\mathop{q}\limits^\tilde}_{F\mathcal{l}}{L}_{F\mathcal{l}})\right]}^{-1}, $$(A.10)where the symbol q ̃ F l $ {\mathop{q}\limits^\tilde}_{F\mathcal{l}}$ is now the right-hand side of equation (A.4) for each of the N fractures. Basically, in arriving at equation (A.10) we have applied Duhamel’s theorem repeatedly. The above developments permit us to examine several asymptotic solutions. For illustration we first consider Flow Regimes 2 and 3.

A.1.1 Flow Regime 2

As noted earlier, Flow Regime 2 holds when the argument of tanh(x) in equation (15) is such that

tanh [ 3 ω s ̃ g ( η ̃ f ) λ ] 1 ; $$ \mathrm{tanh}\left[\sqrt{\frac{3{\omega }^\mathrm{\prime}\mathop{s}\limits^\tildeg({\stackrel{\tilde }{\eta }}_f)}{{\lambda }^\mathrm{\prime}}}\right]\approx 1; $$(A.11)and f ( s ̃ ) $ f(\mathop{s}\limits^\tilde)$ is given by

f ( s ̃ ) = λ ω 3 s ̃ g ( η ̃ f ) . $$ f(\mathop{s}\limits^\tilde)=\sqrt{\frac{{\lambda }^\mathrm{\prime}{\omega }^\mathrm{\prime}}{3\mathop{s}\limits^\tildeg({\stackrel{\tilde }{\eta }}_f)}}. $$(A.12)

Using this result we consider the variation of the second tanh(x) of the solution developed above, namely tanh q ̃ F $ \mathrm{tanh}{\mathop{q}\limits^\tilde}_F$ in equation (A.6) to illustrate responses for the Bilinear and Linear flow regimes.

A.1.2 Bilinear flow

For time ranges when tanh q ̃ F L F $ \mathrm{tanh}{\mathop{q}\limits^\tilde}_F{L}_F$ is ≈1 (Eq. (A.9) or Eq. (A.10)), the expression for the pressure response simplifies considerably. For example, in the case of equation (A.10) we have

p ¯ wD ( s ) = π s [ l = 1 N q ̃ F l C FD l ] - 1 . $$ {\bar{p}}_{{wD}}^{\star }(s)=\frac{\pi }{s}{\left[\sum_{\mathcal{l}=1}^N {\mathop{q}\limits^\tilde}_{F\mathcal{l}}{C}_{{FD}\mathcal{l}}^{\star }\right]}^{-1}. $$(A.13)

On substituting the right-hand side of equation (A.4) for q ̃ F $ {\mathop{q}\limits^\tilde}_F$ after noting the dependence on l $ \mathcal{l}$, equation (A.13) is

p ¯ wD ( s ) = π s 1 l = 1 N q ̃ F l C FD l , $$ {\bar{p}}_{{wD}}^{\star }(s)=\frac{\pi }{s}\frac{1}{\sum_{\mathcal{l}=1}^N {\mathop{q}\limits^\tilde}_{F\mathcal{l}}{C}_{{FD}\mathcal{l}}^{\star }}, $$(A.14)with q ̃ F l $ {\mathop{q}\limits^\tilde}_{F\mathcal{l}}$ given as

q ̃ F l 2 = ( 2 k α f , β f k F l w F l s 1 - α f λ ̃ β ̃ f + s η F ) . $$ {\mathop{q}\limits^\tilde}_{F\mathcal{l}}^2=\left(2\frac{{k}_{{\alpha }_f,{\beta }_f}}{{k}_{F\mathcal{l}}{w}_{F\mathcal{l}}}{s}^{1-{\alpha }_f}{\stackrel{\tilde }{\lambda }}^{{\stackrel{\tilde }{\beta }}_f}+\frac{s}{{\eta }_F}\right). $$(A.15)

On extracting the first term from the parentheses on the right-hand side of equation (A.15), expanding the resulting expression by the binomial theorem and dropping all but the leading terms we replace q ̃ F l $ {\mathop{q}\limits^\tilde}_{F\mathcal{l}}$ in equation (A.10) by the final result to obtain

p ¯ wD ( s ) = π s 1 2 l = 1 N C FD l s 1 [ s 1 - α f λ ̃ β f β f + 1 ] 1 2 . $$ {\bar{p}}_{{wD}}^{\star }(s)=\frac{\pi }{s}\frac{1}{\sqrt{2\sum_{\mathcal{l}=1}^N {C}_{{FD}{\mathcal{l}}_s}^{\star }}}\frac{1}{{\left[{s}^{1-{\alpha }_f}{\stackrel{\tilde }{\lambda }}^{\frac{{\beta }_f}{{\beta }_f+1}}\right]}^{\frac{1}{2}}}. $$(A.16)

Now on using the expression for f ( s ̃ ) $ f(\mathop{s}\limits^\tilde)$ in equation (A.12), we obtain the dimensionless pressure response, pwD(tD), in terms of dimensionless time, tD, and dimensionless fracture conductivity, C FD l $ {C}_{{FD}\mathcal{l}}$, to be

p wD ( t D ) = π 2 l = 1 N C FD l ( 3 λ ω ) β t D 2 ( β f + 1 ) - 2 α f - α m β f 4 α f ( β f + 1 ) Γ [ 1 + 2 ( β f + 1 ) - 2 α f - α m β f 4 ( β f + 1 ) ] , $$ {p}_{{wD}}({t}_D)=\frac{\pi }{\sqrt{2\sum_{\mathcal{l}=1}^N {C}_{{FD}\mathcal{l}}}}{\left(\frac{3}{{\lambda }^\mathrm{\prime}{\omega }^\mathrm{\prime}}\right)}^{{\beta }^{\star }}\frac{{t}_D^{\frac{2({\beta }_f+1)-2{\alpha }_f-{\alpha }_m{\beta }_f}{4{\alpha }_f({\beta }_f+1)}}}{\mathrm{\Gamma }\left[1+\frac{2({\beta }_f+1)-2{\alpha }_f-{\alpha }_m{\beta }_f}{4({\beta }_f+1)}\right]}, $$(A.17)with

β = β f 4 ( β f + 1 ) . $$ {\beta }^{\star }=\frac{{\beta }_f}{4({\beta }_f+1)}. $$(A.18)

The symbols pwD(tD), tD, and C FD l $ {C}_{{FD}\mathcal{l}}$, are defined, respectively, by the following expressions

p wD ( t D ) = 2 π k α f , β f h t η ̃ f 1 - α f α f l ref β f + 1 - 2 α f α f Δ p w ( t ) , $$ {p}_{{wD}}({t}_D)=\frac{2\pi {k}_{{\alpha }_f,{\beta }_f}{h}_t}{{q\mu }}\frac{{\stackrel{\tilde }{\eta }}_f^{\frac{1-{\alpha }_f}{{\alpha }_f}}}{{\mathcal{l}}_{{ref}}^{\frac{{\beta }_f+1-2{\alpha }_f}{{\alpha }_f}}}\Delta {p}_w(t), $$(A.19)

t D = η ̃ l ref β f + 1 t α f , $$ {t}_D=\frac{\stackrel{\tilde }{\eta }}{{\mathcal{l}}_{{ref}}^{{\beta }_f+1}}{t}^{{\alpha }_f}, $$(A.20)and

C FD l = k F l w F l k α f , β f l ref l ref β f + 1 - 2 α f α f η ̃ f 1 - α f α f . $$ {C}_{{FD}\mathcal{l}}=\frac{{k}_{F\mathcal{l}}{w}_{F\mathcal{l}}}{{k}_{{\alpha }_f,{\beta }_f}{\mathcal{l}}_{{ref}}}\frac{{\mathcal{l}}_{{ref}}^{\frac{{\beta }_f+1-2{\alpha }_f}{{\alpha }_f}}}{{\stackrel{\tilde }{\eta }}_f^{\frac{1-{\alpha }_f}{{\alpha }_f}}}. $$(A.21)

If all three exponents (αf, βf and βm) are unity; that is, classical diffusion prevails, then the exponent of t in equation (A.17) is 1/8 reflecting the influence of the matrix (Yeh et al., 1986).

A.1.3 Linear flow

This flow regime reflects situations in which tanh q ̃ F q ̃ F $ \mathrm{tanh}{\mathop{q}\limits^\tilde}_F\approx {\mathop{q}\limits^\tilde}_F$ or tanh(x) ≈ x − x 3/3 for tanh q ̃ F $ \mathrm{tanh}{\mathop{q}\limits^\tilde}_F$ (Camacho-Velázquez, 1984; Camacho-Velázquez et al., 1987). We follow the latter option, and using the right-hand side of equation (A.4) for q ̃ F 2 $ {\mathop{q}\limits^\tilde}_F^2$ we may write

l = 1 N q ̃ F l k F l w F l tanh ( q ̃ F l L F l ) [ l = 1 N k F l w F l Ω 1 l L F l ( 1 + Ω 2 l ) ] [ 1 - l = 1 N k F l w F l Ω 3 l ( 1 + Ω 2 l ) 2 ] , $$ \begin{array}{ll}\sum_{\mathcal{l}=1}^N {\mathop{q}\limits^\tilde}_{F\mathcal{l}}{k}_{F\mathcal{l}}{w}_{F\mathcal{l}}\mathrm{tanh}({\mathop{q}\limits^\tilde}_{F\mathcal{l}}{L}_{F\mathcal{l}})& \approx \left[\sum_{\mathcal{l}=1}^N {k}_{F\mathcal{l}}{w}_{F\mathcal{l}}{\mathrm{\Omega }}_{1\mathcal{l}}{L}_{F\mathcal{l}}(1+{\mathrm{\Omega }}_{2\mathcal{l}})\right]\\ & \\ & \left[1-\sum_{\mathcal{l}=1}^N {k}_{F\mathcal{l}}{w}_{F\mathcal{l}}{\mathrm{\Omega }}_{3\mathcal{l}}(1+{\mathrm{\Omega }}_{2\mathcal{l}}{)}^2\right],\end{array} $$(A.22)where

Ω 1 l = 2 k α f , β f k F l w F l s a λ ̃ β ̃ f , $$ {\mathrm{\Omega }}_{1\mathcal{l}}=2\frac{{k}_{{\alpha }_f,{\beta }_f}}{{k}_{F\mathcal{l}}{w}_{F\mathcal{l}}}{s}^a{\stackrel{\tilde }{\lambda }}^{{\stackrel{\tilde }{\beta }}_f}, $$(A.23)

Ω 2 l = s 1 - a 2 η F Ω 1 l , $$ {\mathrm{\Omega }}_{2\mathcal{l}}=\frac{{s}^{1-a}}{2{\eta }_F{\mathrm{\Omega }}_{1\mathcal{l}}}, $$(A.24)

Ω 3 l = 2 s 1 - a 3 Ω 1 l 2 L F 3 k = 1 N k Fk w Fk Ω Fk L Fk $$ {\mathrm{\Omega }}_{3\mathcal{l}}=\frac{2{s}^{1-a}}{3}\frac{{\mathrm{\Omega }}_{1\mathcal{l}}^2{L}_F^3}{\sum_{k=1}^N {k}_{{Fk}}{w}_{{Fk}}{\mathrm{\Omega }}_{{Fk}}{L}_{{Fk}}} $$(A.25)with β ̃ f $ {\stackrel{\tilde }{\beta }}_f$ and a, respectively, given as

β ̃ f = β f β f + 1 , $$ {\stackrel{\tilde }{\beta }}_f=\frac{{\beta }_f}{{\beta }_f+1}, $$(A.26)and

a = 1 - α f ( 1 - β ̃ f ) . $$ a=1-{\alpha }_f(1-{\stackrel{\tilde }{\beta }}_f). $$(A.27)

With the result noted in equation (A.22) we may obtain well responses for linear flow, again, on substituting appropriate expressions for λ ̃ β ̃ f $ {\stackrel{\tilde }{\lambda }}^{{\stackrel{\tilde }{\beta }}_f}$ and we do so on the assumption that Ω 2 l = 1 $ {\mathrm{\Omega }}_{2\mathcal{l}}=1$. On using the right-hand side of equation (A.22) in equation (A.10) for l = 1 N q ̃ F l k F l w F l tanh ( q ̃ F l L F l ) $ {\sum }_{\mathcal{l}=1}^N {\mathop{q}\limits^\tilde}_{F\mathcal{l}}{k}_{F\mathcal{l}}{w}_{F\mathcal{l}}\mathrm{tanh}({\mathop{q}\limits^\tilde}_{F\mathcal{l}}{L}_{F\mathcal{l}})$, along with the expression in equation (A.12) for f ( s ̃ ) $ f(\mathop{s}\limits^\tilde)$, simplifying and inverting the result, the dimensionless pressure, pwD(tD), during linear flow under the constraint that Flow Regime 2 $ 2$ prevails is given by

p wD ( t D ) = π 2 ( 3 λ ω ) β 1 j = 1 n L Dj t D 2 ( β f + 1 ) - 2 α f - α m β f 2 α f ( β f + 1 ) Γ [ 1 + 2 ( β f + 1 ) - 2 α f - α m β f 2 ( β f + 1 ) ]   + π 3 j = 1 n ( L Dj 3 / C fDj ) ( j = 1 n L Dj ) 2 , $$ \begin{array}{ll}{p}_{{wD}}({t}_D)& =\frac{\pi }{2}{\left(\frac{3}{{\lambda }^\mathrm{\prime}{\omega }^\mathrm{\prime}}\right)}^{{\beta }^{\star \star }}\frac{1}{\sum_{j=1}^n {L}_{{Dj}}}\frac{{t}_D^{\frac{2({\beta }_f+1)-2{\alpha }_f-{\alpha }_m{\beta }_f}{2{\alpha }_f({\beta }_f+1)}}}{\mathrm{\Gamma }\left[1+\frac{2({\beta }_f+1)-2{\alpha }_f-{\alpha }_m{\beta }_f}{2({\beta }_f+1)}\right]}\\ & \enspace +\frac{\pi }{3}\frac{\sum_{j=1}^n \left({L}_{{Dj}}^3/{C}_{{fDj}}\right)}{{\left(\sum_{j=1}^n {L}_{{Dj}}\right)}^2},\end{array} $$(A.28)with

β = β f 2 ( β f + 1 ) , $$ {\beta }^{\star \star }=\frac{{\beta }_f}{2({\beta }_f+1)}, $$(A.29)and

L Dj = L fj L f 1 , $$ {L}_{{Dj}}=\frac{{L}_{{fj}}}{{L}_{f1}}, $$(A.30)where LDj reflects the fact that the lengths of the fractures may be variable. Equation (A.28) indicates that the slopes of the straight line reflecting linear flow during Flow Regime 2 will, again, be much smaller than the slope of 1/4 that would result if the classical diffusion were to prevail. By the same token, for high conductivity or planar fractures we have

p wD ( t D ) = π 2 ( 3 λ ω ) β 1 j = 1 n L Dj t D 2 ( β f + 1 ) - 2 α f - α m β f 2 α f ( β f + 1 ) Γ [ 1 + 2 ( β f + 1 ) - 2 α f - α m β f 2 ( β f + 1 ) ] . $$ {p}_{{wD}}({t}_D)=\frac{\pi }{2}{\left(\frac{3}{{\lambda }^\mathrm{\prime}{\omega }^\mathrm{\prime}}\right)}^{{\beta }^{\star \star }}\frac{1}{\sum_{j=1}^n {L}_{{Dj}}}\frac{{t}_D^{\frac{2({\beta }_f+1)-2{\alpha }_f-{\alpha }_m{\beta }_f}{2{\alpha }_f({\beta }_f+1)}}}{\mathrm{\Gamma }\left[1+\frac{2({\beta }_f+1)-2{\alpha }_f-{\alpha }_m{\beta }_f}{2({\beta }_f+1)}\right]}. $$(A.31)

A.1.4 Flow Regime 3

In this section we address our expectations of well responses for the linear-flow regime. As noted earlier, Flow Regime 3 (and also Flow Regime 5, see below) applies for time ranges when the argument of tanh(x) in equation (15) is such that

tanh [ 3 ω s ̃ g ( η ̃ f ) λ ] [ 3 ω s ̃ g ( η ̃ f ) λ ] ; $$ \mathrm{tanh}\left[\sqrt{\frac{3{\omega }^\mathrm{\prime}\mathop{s}\limits^\tildeg({\stackrel{\tilde }{\eta }}_f)}{{\lambda }^\mathrm{\prime}}}\right]\approx \left[\sqrt{\frac{3{\omega }^\mathrm{\prime}\mathop{s}\limits^\tildeg({\stackrel{\tilde }{\eta }}_f)}{{\lambda }^\mathrm{\prime}}}\right]; $$(A.32)consequently, f ( s ̃ ) $ f(\mathop{s}\limits^\tilde)$ may be written as

f ( s ̃ ) = 1 + ω . $$ f\left(\mathop{s}\limits^\tilde\right)=1+{\omega }^\mathrm{\prime}. $$(A.33)

With the expression in equation (A.22), we may write

p ̅ wD ( s ) = π η ̃ β ̃ f 2 s 1 + a l = 1 N k F l w F l L F l Ω F l ( 1 + l = 1 N 2 3 k F l w F l Ω F l 2 L F l 3 s a l = 1 N k F l w F l L F l Ω F l ) , $$ {\bar{p}}_{{wD}}^{\star }(s)=\frac{\pi {\stackrel{\tilde }{\eta }}^{{\stackrel{\tilde }{\beta }}_f}}{2{s}^{1+a}\sum_{\mathcal{l}=1}^N {k}_{F\mathcal{l}}{w}_{F\mathcal{l}}{L}_{F\mathcal{l}}{\mathrm{\Omega }}_{F\mathcal{l}}}\left(1+\frac{\sum_{\mathcal{l}=1}^N \frac{2}{3}{k}_{F\mathcal{l}}{w}_{F\mathcal{l}}{\mathrm{\Omega }}_{F\mathcal{l}}^2{L}_{F\mathcal{l}}^3{s}^a}{\sum_{\mathcal{l}=1}^N {k}_{F\mathcal{l}}{w}_{F\mathcal{l}}{L}_{F\mathcal{l}}{\mathrm{\Omega }}_{F\mathcal{l}}}\right), $$(A.34)where

Ω F l = ( 1 + ω ) β ̃ f k F l w F l . $$ {\mathrm{\Omega }}_{F\mathcal{l}}=\frac{(1+\omega \mathrm{\prime}{)}^{{\stackrel{\tilde }{\beta }}_f}}{{k}_{F\mathcal{l}}{w}_{F\mathcal{l}}}. $$(A.35)

Consequently, the well response, pwD(tD), is given by

p wD ( t D ) = π 2 ( 1 + ω ) β ̃ f 1 j = 1 n L Dj t D ( β f + 1 ) - α f α f ( β f + 1 ) Γ [ 1 + ( β f + 1 ) - α f ( β f + 1 ) ] + π 3 j = 1 n L Dj 3 / C fDj ( j = 1 n L Dj ) 2 . $$ {p}_{{wD}}({t}_D)=\frac{\pi }{2(1+{\omega }^\mathrm{\prime}{)}^{{\stackrel{\tilde }{\beta }}_f}}\frac{1}{\sum_{j=1}^n {L}_{{Dj}}}\frac{{t}_D^{\frac{({\beta }_f+1)-{\alpha }_f}{{\alpha }_f({\beta }_f+1)}}}{\mathrm{\Gamma }\left[1+\frac{({\beta }_f+1)-{\alpha }_f}{({\beta }_f+1)}\right]}+\frac{\pi }{3}\frac{\sum_{j=1}^n {L}_{{Dj}}^3/{C}_{{fDj}}}{{\left(\sum_{j=1}^n {L}_{{Dj}}\right)}^2}. $$(A.36)

During Flow Regime 1 the well performs as if the matrix system were nonexistent; consequently, the response during that period may be obtained from equation (A.36) by setting ω′ = 0. In other words, much like the Barenblatt et al. (1960) formulation, responses during Flow Regime 3 parallel that of Flow Regime 1 with the displacement being proportional to ( 1 + ω ) β ̃ f $ (1+{\omega }^\mathrm{\prime}{)}^{{\stackrel{\tilde }{\beta }}_f}$.

A.1.5 Flow Regime 4 and Flow Regime 5

We consider responses to a single fracture for times that are long enough; also, if times are long enough the influence of the conductivity of the fracture may be ignored. The needed expressions are derived directly from the development given here rather than the indirect approach used in Raghavan and Chen (2018a). We then compare the derived expression with that obtained for subdiffusion. As subdiffusive solutions obtained are exact, again as in the text, parallels are drawn to the expression for the pressure distribution in 2D.

As the solutions involve the ratio of the Mittag-Leffler functions, for small values of s, the ratio of Mittag-Leffler Functions, E α 1 , β 1 ( z ) $ {E}_{{\alpha }_1,{\beta }_1}(z)$ and E α 2 , β 2 ( z ) $ {E}_{{\alpha }_2,{\beta }_2}(z)$, may be written as

E α 1 , β 1 ( z ) E α 2 , β 2 ( z ) = 1 [ Γ ( β 1 ) ] 2 1 [ 1 Γ ( β 1 ) - z Γ ( α 1 + β 1 ) ] [ 1 Γ ( β 2 ) + z Γ ( α 2 + β 2 ) ] , $$ \frac{{E}_{{\alpha }_1,{\beta }_1}(z)}{{E}_{{\alpha }_2,{\beta }_2}(z)}=\frac{1}{[\mathrm{\Gamma }({\beta }_1){]}^2}\frac{1}{\left[\frac{1}{\mathrm{\Gamma }({\beta }_1)}-\frac{z}{\mathrm{\Gamma }({\alpha }_1+{\beta }_1)}\right]\left[\frac{1}{\mathrm{\Gamma }({\beta }_2)}+\frac{z}{\mathrm{\Gamma }({\alpha }_2+{\beta }_2)}\right]}, $$(A.37)or

E α 1 , β 1 ( z ) E α 2 , β 2 ( z ) = 1 [ Γ ( β 1 ) ] 2 Γ ( β 1 ) Γ ( β 2 ) a ( z + 1 a ) , $$ \frac{{E}_{{\alpha }_1,{\beta }_1}(z)}{{E}_{{\alpha }_2,{\beta }_2}(z)}=\frac{1}{[\mathrm{\Gamma }({\beta }_1){]}^2}\frac{\mathrm{\Gamma }({\beta }_1)\mathrm{\Gamma }({\beta }_2)}{a\left(z+\frac{1}{a}\right)}, $$(A.38)where

a = Γ ( β 2 ) Γ ( α 1 + β 1 ) - Γ ( β 1 ) Γ ( α 2 + β 2 ) Γ ( α 1 + β 1 ) Γ ( α 2 + β 2 ) . $$ a=\frac{\mathrm{\Gamma }({\beta }_2)\mathrm{\Gamma }({\alpha }_1+{\beta }_1)-\mathrm{\Gamma }({\beta }_1)\mathrm{\Gamma }({\alpha }_2+{\beta }_2)}{\mathrm{\Gamma }({\alpha }_1+{\beta }_1)\mathrm{\Gamma }({\alpha }_2+{\beta }_2)}. $$(A.39)

Thus the approximation for φ in equation (44) is then given by

φ = λ ̃ y ̃ e a ̃ ( λ ̃ y ̃ e + 1 a ̃ ) , $$ \phi =\frac{\stackrel{\tilde }{\lambda }{\mathop{y}\limits^\tilde}_e}{\mathop{a}\limits^\tilde(\stackrel{\tilde }{\lambda }{\mathop{y}\limits^\tilde}_e+\frac{1}{\mathop{a}\limits^\tilde})}, $$(A.40)where

a ̃ = 1 ( β f + 2 ) Γ ( β f + 1 ) . $$ \mathop{a}\limits^\tilde=\frac{1}{({\beta }_f+2)\mathrm{\Gamma }({\beta }_f+1)}. $$(A.41)

On again using a two-term approximation for the Eα,β(z) expressions, we may write equation (43) with the aid of the expression in equation (A.40) to be

Δ p ̅ f ( y ̃ , s ) = - [ d β f d y ̃ β f   Δ p ̅ f ( y ̃ , s ) ] y ̃ = 0 y ̃ e { η ̃ f s α f f ( s ) + [ y ̃ β f + 1 Γ ( β f + 2 ) - y ̃ e y ̃ β f Γ ( β f + 1 ) + y ̃ e β f + 1 ( β f + 2 ) Γ ( β f + 1 ) ] } . $$ \begin{array}{ll}{\overline{\Delta p}}_f(\mathop{y}\limits^\tilde,s)& =-\frac{{\left[\frac{{d}^{{\beta }_f}}{d{\mathop{y}\limits^\tilde}^{{\beta }_f}}\enspace {\overline{\Delta p}}_f(\mathop{y}\limits^\tilde,s)\right]}_{\mathop{y}\limits^\tilde=0}}{{\mathop{y}\limits^\tilde}_e}\\ & \left\{\frac{{\stackrel{\tilde }{\eta }}_f}{{s}^{{\alpha }_f}f(s)}+\left[\frac{{\mathop{y}\limits^\tilde}^{{\beta }_f+1}}{\mathrm{\Gamma }({\beta }_f+2)}-\frac{{\mathop{y}\limits^\tilde}_e{\mathop{y}\limits^\tilde}^{{\beta }_f}}{\mathrm{\Gamma }({\beta }_f+1)}+\frac{{\mathop{y}\limits^\tilde}_e^{{\beta }_f+1}}{({\beta }_f+2)\mathrm{\Gamma }({\beta }_f+1)}\right]\right\}.\end{array} $$(A.42)

In the above expression, the flow rate, q ̅ , is [ d β f   Δ p ̅ f ( y ̃ , s ) ] / d y ̃ β f | y ̃ = 0 $ \bar{q},{is}\propto \left[{d}^{{\beta }_f}\enspace {\overline{\Delta p}}_f(\mathop{y}\limits^\tilde,s)\right]/d{\mathop{y}\limits^\tilde}^{{\beta }_f}{|}_{\mathop{y}\limits^\tilde=0}$. In order to relate the solutions to subdiffusion, we now modify the solution in equation (A.42), replace y ̃ $ \mathop{y}\limits^\tilde$ by y with the presumption that |y| ≥ 0 represents the distance from the fracture plane located at y = 0 and consider a system of drainage area, A,

A = L F y e . $$ A={L}_F{y}_e. $$(A.43)

The expression for dimensionless pressure is thus given by

p ̅ D ( 0 x D x eD , y D , s 0 ) = π 2 { ( η ̃ f A ̃ ) 1 α f 1 s 2 f ( s ) + [ | y | D β f + 1 Γ ( β f + 2 ) - y eD | y | D β f Γ ( β f + 1 ) + y eD β f + 1 ( β f + 2 ) Γ ( β f + 1 ) ] ( η ̃ f A ̃ ) 1 - α f α f 1 s 2 - α f } . $$ \begin{array}{ll}{\bar{p}}_D(0\le {x}_D\le {x}_{{eD}},{y}_D,s\to 0)& =\frac{\pi }{2}\left\{{\left(\frac{{\stackrel{\tilde }{\eta }}_f}{\mathop{A}\limits^\tilde}\right)}^{\frac{1}{{\alpha }_f}}\frac{1}{{s}^2f(s)}+\left[\frac{|y{|}_D^{{\beta }_f+1}}{\mathrm{\Gamma }({\beta }_f+2)}-\frac{{y}_{{eD}}|y{|}_D^{{\beta }_f}}{\mathrm{\Gamma }({\beta }_f+1)}\right.\right.\\ & \left.\left.+\frac{{{y}_{{eD}}}^{{\beta }_f+1}}{({\beta }_f+2)\mathrm{\Gamma }({\beta }_f+1)}\right]{\left(\frac{{\stackrel{\tilde }{\eta }}_f}{\mathop{A}\limits^\tilde}\right)}^{\frac{1-{\alpha }_f}{{\alpha }_f}}\frac{1}{{s}^{2-{\alpha }_f}}\right\}.\end{array} $$(A.44)where

A ̃ = A ( β f + 1 ) / 2 . $$ \mathop{A}\limits^\tilde={A}^{({\beta }_f+1)/2}. $$(A.45)

The analog of equation (54) for constant-rate production is

Δ p ̅ f ( y ̃ , s ) = 4 k α f β f h L F 1 s 2 - α f cosh [ λ ̃ ( y ̃ e - y ̃ ) ] λ ̃ sinh ( λ ̃ y ̃ e ) . $$ {\overline{\Delta p}}_f(\mathop{y}\limits^\tilde,s)=\frac{{q\mu }}{4{k}_{{\alpha }_f{\beta }_f}h{L}_F}\frac{1}{{s}^{2-{\alpha }_f}}\frac{\mathrm{cosh}[\sqrt{\stackrel{\tilde }{\lambda }}({\mathop{y}\limits^\tilde}_e-\mathop{y}\limits^\tilde)]}{\stackrel{\tilde }{\lambda }\mathrm{sinh}(\sqrt{\stackrel{\tilde }{\lambda }}{\mathop{y}\limits^\tilde}_e)}. $$(A.46)

Using the result given on page 47 (1.445,2) of Gradshteyn and Ryzhik (1965), namely,

k = 1 cos kx k 2 + a 2 = π 2 a cosh [ a ( π - x ) ] sinh ( ) - 1 2 a 2 ,   [ 0 x 2 π ] , $$ \sum_{k=1}^{\mathrm{\infty }} \frac{\mathrm{cos}{kx}}{{k}^2+{a}^2}=\frac{\pi }{2a}\frac{\mathrm{cosh}[a(\pi -x)]}{\mathrm{sinh}({a\pi })}-\frac{1}{2{a}^2},\enspace \hspace{1em}[0\le x\le 2\pi ], $$(A.47)we recast equation (A.46) as

p ̅ D ( 0 x D x eD , y D , s 0 ) = π 2 { ( η ̃ f A ̃ ) 1 α f 1 s 2 f ( s ) + 2 y ̃ e 2 π 2 A ̃ k = 1 cos y / y ̃ e k 2 + λ ̃ y e / π 2 ( η ̃ A ̃ ) 1 - α f α f 1 s 2 - α f } . $$ {\bar{p}}_D(0\le {x}_D\le {x}_{{eD}},{y}_D,s\to 0)=\frac{\pi }{2}\left\{{\left(\frac{{\stackrel{\tilde }{\eta }}_f}{\mathop{A}\limits^\tilde}\right)}^{\frac{1}{{\alpha }_f}}\frac{1}{{s}^2f(s)}+\frac{2{\mathop{y}\limits^\tilde}_e^2}{{\pi }^2\mathop{A}\limits^\tilde}\sum_{k=1}^{\mathrm{\infty }} \frac{\mathrm{cos}{k\pi y}/{\mathop{y}\limits^\tilde}_e}{{k}^2+\stackrel{\tilde }{\lambda }{y}_e/{\pi }^2}{\left(\frac{\stackrel{\tilde }{\eta }}{\mathop{A}\limits^\tilde}\right)}^{\frac{1-{\alpha }_f}{{\alpha }_f}}\frac{1}{{s}^{2-{\alpha }_f}}\right\}. $$(A.48)

If we consider the solution given in equation (A.48) as λ ̃ 0 $ \stackrel{\tilde }{\lambda }\to 0$, and use the result of Gradshteyn and Ryzhik (1965) given on page 46 (1.445,2) namely that

k = 1 cos kx k 2 = π 2 6 - π x 2 + x 2 4 ,   [ 0 x 2 π ] , $$ \sum_{k=1}^{\mathrm{\infty }} \frac{\mathrm{cos}{kx}}{{k}^2}=\frac{{\pi }^2}{6}-\frac{{\pi x}}{2}+\frac{{x}^2}{4},\enspace \hspace{1em}[0\le x\le 2\pi ], $$(A.49)we will arrive at the result in equation (A.44) for βf = 1.

The belaboring which has occurred to establish the correspondences noted here becomes more evident if we consider subdiffusive solutions in 2D. The pressure distribution in a rectangle wherein the fracture length is not equal to its width (xe/xf ≠ 1) as shown in Figure 4 is given as

Δ p ̅ f ( x , y ) = π x e s 1 - α f x w - x F 1 x w + x F 2 dx q ̃ ( x ) [ ch λ ̃ y e 1 + ch λ ̃ y e 2 λ ̃   sh λ ̃ y e + 2 k = 1 cos   x x e cos   x x e ch ε k y e 1 + ch ε k y e 2 ε k sh ε k y e ] , $$ \begin{array}{ll}{\overline{\Delta p}}_f(x,y)& =\frac{\pi }{{x}_e{s}^{1-{\alpha }_f}}{\int }_{{x}_w-{x}_{F1}}^{{x}_w+{x}_{F2}} {dx}\mathrm{\prime}\mathop{q}\limits^\tilde(x\mathrm{\prime})\\ & \left[\frac{{ch}\sqrt{\stackrel{\tilde }{\lambda }}{y}_{e1}+{ch}\sqrt{\stackrel{\tilde }{\lambda }}{y}_{e2}}{\sqrt{\stackrel{\tilde }{\lambda }}\enspace {sh}\sqrt{\stackrel{\tilde }{\lambda }}{y}_e}\right.\\ & +2\sum_{k=1}^{\mathrm{\infty }} \mathrm{cos}\enspace {k\pi }\frac{x\mathrm{\prime}{{x}_e}\mathrm{cos}\enspace {k\pi }\frac{x}{{x}_e}\\ & \left.\frac{{ch}{\sqrt{\epsilon }}_k{y}_{e1}+{ch}{\sqrt{\epsilon }}_k{y}_{e2}}{{\sqrt{\epsilon }}_k{sh}{\sqrt{\epsilon }}_k{y}_e}\right],\end{array} $$(A.50)where q ̃ ( x ) $ \mathop{q}\limits^\tilde(x)$ is the flux distribution along the fracture face, ε ̃ k = λ ̃ + k 2 π 2 / x e 2 $ {\stackrel{\tilde }{\epsilon }}_k=\stackrel{\tilde }{\lambda }+{k}^2{\pi }^2/{x}_e^2$, y ̃ e 1 = | y - y w | $ {\mathop{y}\limits^\tilde}_{e1}=|y-{y}_w|$, and y ̃ e 2 = y + y w $ {\mathop{y}\limits^\tilde}_{e2}=y+{y}_w$, and xF = xF1 + xF2.

In equation (A.50), at long enough times, the well response for xe/xf of 1, the situation corresponding to where the fracture length equals the width of the rectangle, is given by the first term and will be identical to that given in equation (A.42) for βf = 1. This implies that all hyperbolic terms should, perhaps, be replaced by Mittag-Leffler functions. A glimpse of the structure of the solution for superdiffusion then becomes available. All this, of course, would have to be verified.

A.1.5.1 Pressure distributions, constant rate

Here we use the expressions in equation (A.12) and in equation (A.33) for f ( s ̃ ) $ f(\mathop{s}\limits^\tilde)$ in equation (A.44). Inversion of the resulting expressions yields the pressure distributions during Flow Regimes 4 and 5 which are, respectively,

p D ( x D , t DA 0 ) = 2 π { 3 λ ω t DA 2 - α m 2 α f Γ ( 2 - α m / 2 ) + [ y D β f + 1 Γ ( β f + 2 ) - y LD y D β f Γ ( β f + 1 ) + y eD β f + 1 ( β f + 2 ) Γ ( β f + 1 ) ] t D 1 - α f α f Γ ( 2 - α f ) } , $$ \begin{array}{ll}{p}_D({x}_D,{t}_{{DA}}\to 0)& =2\pi \left\{\sqrt{\frac{3}{{\lambda }^\mathrm{\prime}{\omega }^\mathrm{\prime}}}\frac{{t}_{{DA}}^{\frac{2-{\alpha }_m}{2{\alpha }_f}}}{\mathrm{\Gamma }(2-{\alpha }_m/2)}+\right.\\ & \left.\left[\frac{{y}_D^{{\beta }_f+1}}{\mathrm{\Gamma }({\beta }_f+2)}-\frac{{y}_{{LD}}{y}_D^{{\beta }_f}}{\mathrm{\Gamma }({\beta }_f+1)}+\frac{{y}_{{eD}}^{{\beta }_f+1}}{({\beta }_f+2)\mathrm{\Gamma }({\beta }_f+1)}\right]\frac{{t}_D^{\frac{1-{\alpha }_f}{{\alpha }_f}}}{\mathrm{\Gamma }(2-{\alpha }_f)}\right\},\end{array} $$(A.51)and

p D ( x D , t DA ) = 2 π { 1 ( 1 + ω ) t DA 1 α f + [ y D β f + 1 Γ ( β f + 2 ) - y eD x D β f Γ ( β f + 1 ) + y eD β f + 1 ( β f + 2 ) Γ ( β f + 1 ) ] t DA 1 - α f α f Γ ( 2 - α f ) } . $$ \begin{array}{ll}{p}_D({x}_D,{t}_{{DA}}\to \mathrm{\infty })& =2\pi \left\{\frac{1}{(1+\omega \mathrm{\prime})}{t}_{{DA}}^{\frac{1}{{\alpha }_f}}+\right.\\ & \left.\left[\frac{{y}_D^{{\beta }_f+1}}{\mathrm{\Gamma }({\beta }_f+2)}-\frac{{y}_{{eD}}{x}_D^{{\beta }_f}}{\mathrm{\Gamma }({\beta }_f+1)}+\frac{{y}_{{eD}}^{{\beta }_f+1}}{({\beta }_f+2)\mathrm{\Gamma }({\beta }_f+1)}\right]\frac{{t}_{{DA}}^{\frac{1-{\alpha }_f}{{\alpha }_f}}}{\mathrm{\Gamma }(2-{\alpha }_f)}\right\}.\end{array} $$(A.52)

Turning our attention to production rate, q(t), for situations wherein the well produces at a constant pressure; that is, when Δ p ̅ f ( y ̃ , s ) = Δ p wf / s $ {\overline{\Delta p}}_f(\mathop{y}\limits^\tilde,s)=\Delta {p}_{{wf}}/s$ at y ̃ = 0 $ \mathop{y}\limits^\tilde=0$ in equation (A.42) as given in van Everdingen and Hurst (1949), we may readily show that for Flow Regime 4, q ( t ) 1 / t ( 2 - α m ) / 2 E α f - α m / 2 , α f ( - a t ( α f - α m / 2 ) ) $ q(t)\sim 1/{t}^{(2-{\alpha }_m)/2}{E}_{{\alpha }_f-{\alpha }_m/2,{\alpha }_f}\left(-a{t}^{({\alpha }_f-{\alpha }_m/2)}\right)$, where a is a constant or q ( t ) 1 / t ( 1 - α m / 2 ) $ q(t)\sim 1/{t}^{(1-{\alpha }_m/2)}$ (q(t) ~ 1/t 1/2) if αm = 1) at terminal conditions because at long enough times Eα,β(−x) ~ 1/x for αβ. Such a power-law trend is reflective of disconnectedness, disjointedness, ruptures, and dead ends as a result of intercalations and the like. Similarly, for Flow Regime 5, q ( t ) 1 / t 1 - α f E α f , α f ( - at ) $ q(t)\sim 1/{t}^{1-{\alpha }_f}{E}_{{\alpha }_f,{\alpha }_f}\left(-{at}\right)$ or a power-law trend of the form q ( t ) 1 / t 1 + α f $ q(t)\sim 1/{t}^{1+{\alpha }_f}$ as Eα,α(−x) ~ 1/x 2 at long enough times provided that α ≠ 1. Should αf = 1; we, of course, as noted, earlier, have q(t) ~ exp(−at).

A.1.5.2 Addressing multiple fractures

Following Spath et al. (1994) we may address the performance of a well producing through multiple fractures through the expression

p ̅ wD ( s ) = 1 s 2 k = 1 N q ̅ wDk , $$ {\bar{p}}_{{wD}}(s)=\frac{1}{{s}^2\sum_{k=1}^N {\bar{q}}_{{wDk}}}, $$(A.53)where qwDk(t) is the flow rate for a well produced at a constant pressure, pwf, through a fracture with the properties of Fracture k. If we restrict our attention to the boundary-dominated period, then equation (A.42) yields the following expression for determining the production rate, qwp(t), at Fracture k:

q ̅ kwp ( s ) = 4 k ̃ α f , β f y e L Fk h t ( p i - p wf ) μ 1 { η ̃ f f ( s ) + y e β f + 1 ( β f + 2 ) Γ ( β f + 1 ) s α f } . $$ {\bar{q}}_{{kwp}}(s)=\frac{4{\mathop{k}\limits^\tilde}_{{\alpha }_f,{\beta }_f}{y}_e{L}_{{Fk}}{h}_t}{({p}_i-{p}_{{wf}})\mu }\frac{1}{\{\frac{{\stackrel{\tilde }{\eta }}_f}{f(s)}+\frac{{y}_e^{{\beta }_f+1}}{({\beta }_f+2)\mathrm{\Gamma }({\beta }_f+1)}{s}^{{\alpha }_f}\}}. $$(A.54)


1

For superdiffusion, many studies solve the pde t p = x β p $ {\mathrm{\partial }}_tp={\mathrm{\partial }}_x^{\beta }p$ for 0 < β < 2. We simply note, however, that the range of β may have to be circumscribed to solve for certain auxiliary conditions. We shall return to this point later.

2

For a discussion of addressing superdiffusion through the exponent α see Metzler et al. (2014). We do not consider this perspective. It is an open problem for situations considered here.

3

The fractional flux law for space fractional diffusion may be stated in other forms. Often a pde in the form of t p = x β p $ {\mathrm{\partial }}_tp={\mathrm{\partial }}_x^{\beta }p$ with 0 < β < 2 is considered to model superdiffusion. In such cases a fractional flux law of the form v ( β ) = ( k 1 , β / μ ) x β - 1 p   $ v(\beta )=({k}_{1,\beta }/\mu ){\mathrm{\partial }}_x^{\beta -1}p\enspace $ needs to be used; see Benson et al. (2004). For boundary conditions of the kind considered in this work the applicable range for β would be much less and happens to be identical to the one used in this work.

References

  • Albinali A., Holy R., Sarak H., Ozkan E. (2016) Modeling of 1D anomalous diffusion in fractured nanoporous media, presented at the low permeability media and nanoporous materials from characterization to modeling, Can we do better? Rueil-Malmaison, France, Oil Gas Sci. Technol. – Rev. IFP Energies nouvelles 71, 4, 56. https://doi.org/10.2516/ogst/2016008. [CrossRef] [Google Scholar]
  • Angulo J.M., Ruiz-Medina M.D., Anh V.V., Grecksch W. (2000) Fractional diffusion and fractional heat equation, Adv. Appl. Probab. 32, 4, 1077–1099. [Google Scholar]
  • Artus V. (2020) Numerical upscaling of discrete fracture networks for transient analysis URTEC-2020-3087-MS, in: Presentation at the Unconventional Resources Technology Conference held in Austin, Texas, USA, pp. 20–22. [Google Scholar]
  • Barenblatt G.I., Zheltov Iu.P., Kochina I.N. (1960) Basic concepts in the theory of seepage of homogeneous liquids in fissured rocks, J. Appl. Math. Mech. 24, 5, 1286–1303. [CrossRef] [Google Scholar]
  • Barker J.A. (1988) A Generalized Radial Flow Model for Hydraulic Tests in Fractured Rock, Water Resour. Res. 24, 10, 1796–1804. [Google Scholar]
  • Beier R.A. (1994) Pressure-transient model for a vertically fractured well in a fractal reservoir, SPE Form. Eval. 9, 2, 122–128. https://doi.org/10.2118/20582-PA. [CrossRef] [Google Scholar]
  • Belayneh M., Masihi M., Matthäi S.K., King P.R. (2006) Prediction of vein connectivity using the percolation approach: model test with field data, J. Geophys. Eng. 33, 219–229. https://doi.org/10.1088/1742-2132/3/3/003. [CrossRef] [Google Scholar]
  • Benson D.A., Wheatcraft S.W., Meerschaert M.M. (2000) Application of a fractional advection-dispersion equation, Water Resour. Res. 36, 6, 1403–1412. [Google Scholar]
  • Benson D.A., Tadjeran C., Meerschaert M.M., Farnham I., Pohll G. (2004) Radial fractional-order dispersion through fractured rock, Water Resour. Res. 40, W12416. https://doi.org/10.1029/2004WR003314. [Google Scholar]
  • Benson D.A., Meerschaert M.M., Revielle J. (2013) Fractional calculus in hydrologic modeling: A numerical perspective, Adv. Water Resour. 51, 479–497. [CrossRef] [PubMed] [Google Scholar]
  • Bisdom K., Bertotti G., Nick H.M. (2016) The impact of different aperture distribution models and critical stress criteria on equivalent permeability in fractured rocks, J. Geophys. Res. Solid Earth 121, 5, 2169–9356. https://doi.org/10.1002/2015JB012657. [Google Scholar]
  • Cacas M.C., Ledoux E., de Marsily G., Tillie B., Barbreau A., Durand E., Feuga B., Peaudecerf P. (1990a) Modeling fracture flow with a stochastic discrete fracture network: calibration and validation: 1. The flow model, Water Resour. Res. 26, 3, 479–489. https://doi.org/10.1029/WR026i003p00479. [Google Scholar]
  • Cacas M.C., Ledoux E., de Marsily G., Barbreau A., Calmels P., Gaillard B., Margritta R. (1990b) Modeling fracture flow with a stochastic discrete fracture network: Calibration and validation: 2. The transport model, Water Resour. Res. 26, 3, 491–500. https://doi.org/10.1029/WR026i003p00491. [Google Scholar]
  • Cacas M.C., Daniel J.M., Letouzey J. (2001) Nested geological modelling of naturally fractured reservoirs, Petrol. Geosci. 7, S43–S52. https://doi.org/10.1144/petgeo.7.S.S43. [CrossRef] [Google Scholar]
  • Caine J.S., Evans J.P., Forster C.B. (1996) Fault zone architecture and permeability structure, Geology 24, 11, 1025–1028. [Google Scholar]
  • Camacho-Velázquez R.G. (1984) Response of wells producing commingled reservoirs: Unequal fracture length, Master’s thesis, University of Tulsa, Tulsa, OK. [Google Scholar]
  • Camacho-Velázquez R.G., Raghavan R., Reynolds A.C. (1987) Response of wells producing layered reservoirs: unequal fracture length, SPE Form. Eval. 2, 1, 9–28. https://doi.org/10.2118/12844-PA. [CrossRef] [Google Scholar]
  • Camacho-Velázquez R.G., Raghavan R. (1989) Boundary-dominated flow in solutions-gas-drive reservoirs, SPE Reserv. Eng. 4, 4, 503–512. https://doi.org/10.2118/18562-PA. [CrossRef] [Google Scholar]
  • Caputo M. (1967) Linear models of dissipation whose Q is almost Frequency Independent-II, Geophys. J. Roy. Astron. Soc. 13, 5, 529–539. [NASA ADS] [CrossRef] [Google Scholar]
  • Carrera J., Sánchez-Vila X., Benet I., Medina A., Galarza G., Guimerá J. (1998) On matrix diffusion: formulations, solution methods and qualitative effects, Hydrogeol. J. 6, 1, 178–190. [Google Scholar]
  • Carslaw H.S., Jaeger J.C. (1959) Conduction of heat in solids, 2nd edn., Clarendon Press, Oxford, p. 510. [Google Scholar]
  • Chang J., Yortsos Y.C. (1990) Pressure-transient analysis of Fractal Reservoirs, SPE Form. Eval. 5, 1, 31–39. [CrossRef] [Google Scholar]
  • Chen C. (1982) A study of naturally fractured reservoirs, MS Thesis, University of Tulsa, Tulsa, OK. [Google Scholar]
  • Chen C., Raghavan R. (2013) On some characteristic features of fractured-horizontal wells and conclusions drawn thereof, SPE Reserv. Eval. Eng. 16, 1, 19–28. https://doi.org/10.2118/163104-PA. [CrossRef] [Google Scholar]
  • Chen C., Raghavan R. (2015) Transient flow in a linear reservoir for space-time fractional diffusion, J. Pet. Sci. Eng. 128, 194–202. [Google Scholar]
  • Chow V.T. (1952) On the determination of transmissibility and storage coefficients from pumping test data, Trans. Am. Geophys. Un. 33, 397–404. [CrossRef] [Google Scholar]
  • Chu W., Pandya N., Flumerfelt R.W., Chen C. (2019) Rate-transient analysis based on power-law behavior for Permian wells, SPE Res. Eval. Eng. 22, 4, 1360–1370 https://doi.org/10.2118/187180-PA. [Google Scholar]
  • Chu W., Scott K., Flumerfelt R.W., Chen C. (2020) A new technique for quantifying pressure interference in fractured horizontal shale wells. SPE Res. Eval. Eng. 23, 1, 143–157 https://doi.org/10.2118/191407-PA. [CrossRef] [Google Scholar]
  • Chu W. (2018) Personal Communication. [Google Scholar]
  • Cortis A., Knudby C. (2006) A continuous time random walk approach to transient flow in heterogeneous porous media, LBNL-59885, Water Resour. Res. 42, W10201. [Google Scholar]
  • Cinco-Ley H., Meng H.-Z. (1988) Pressure transient analysis of wells with finite conductivity vertical fractures in double porosity reservoirs, in: Presented at the SPE Annual Technical Conference and Exhibition, 2–5 October, Houston, Texas. https://doi.org/10.2118/18172-MS. [Google Scholar]
  • Dassas Y., Duby Y. (1995) Diffusion toward fractal interfaces, potentiostatic, galvanostatic, and linear sweep voltammetric techniques, J. Electrochem. Soc. 142, 12, 4175–4180. [Google Scholar]
  • de Swaan-O A. (1976) Analytic solutions for determining naturally fractured reservoir properties by well testing, Soc. Pet. Eng. J. 16, 3, 117–122. https://doi.org/10.2118/5346-PA. [CrossRef] [Google Scholar]
  • Defterli O., D’Elia M., Du Q., Gunzburger M., Lehoucq R., Meerschaert M.M. (2015) Fractional Diffusion on Bounded Domains, Fract. Calc. Appl. Anal. 18, 2, 342–360. [Google Scholar]
  • Dershowitz W., Klise K., Fox A., Takeuchi S., Uchida M. (2002) Channel network and discrete fracture network modeling of hydraulic interference and tracer experiments and the prediction of phase C experiments, SKB Report IPR-02-29, SKB, Stockholm. [Google Scholar]
  • Doe T.W. (1991) Fractional dimension analysis of constant-pressure well tests, in: Paper SPE-22702-MS, Presented at the SPE Annual Technical Conference and Exhibition, 6–9 October, Dallas, Texas. https://doi.org/10.2118/22702-MS. [Google Scholar]
  • Doe T., Shi C., Knitter C., Enachescu C. (2014) Discrete fracture network simulation of production data from unconventional wells, in: Paper URTeC 1923802, Proceedings of The Unconventional Resources Technology Conference, Denver, CO. [Google Scholar]
  • Dong Y., Fu Y., Yeh T.-C.J., Wang Y.-L., Zha Y., Wang L., Hao Y. (2019) Equivalence of discrete fracture network and porous media models by hydraulic tomography, Water Resour. Res. 55, 4, 3234–3247. https://doi.org/10.1029/2018WR024290. [Google Scholar]
  • Dontsov K.M., Boyrchuk B.T. (1971) Effect of characteristics of fractured media on pressure buildup behavior, Izvestia VUZ, Oil and Gas N1, 42–46 (in Russian). [Google Scholar]
  • Doyle P.G., Snell J.L. (1984) Random walks and electric networks, Carus Mathematical Monographs, Vol. 22, Mathematical Association of America, 174 pp. https://www.jstor.org/stable/10.4169/j.ctt5hh804. [Google Scholar]
  • Erdelyi A., Magnus W.F., Oberhettinger F., Tricomi F.G. (1955) Higher Transcendental Functions, Chapter 18: Miscellaneous Functions, Vol. 3, McGraw-Hill, New York, pp. 206–227. [Google Scholar]
  • Evans J.P. (1988) Deformation mechanisms in granitic rocks at shallow crustal levels, J. Struct. Geol. 10, 5, 437–443. [Google Scholar]
  • Fomin S., Chugunov V., Hashida T. (2011) Mathematical modeling of anomalous diffusion in porous media, Fract. Differ. Calc. 1, 1–28. [Google Scholar]
  • Gale J.F.W., Laubach S., Olson J.E., Eichhuble P., Fall A. (2014) Natural Fractures in shale: A review and new observations, AAPG Bull. 98, 11, 2165–2216. [CrossRef] [Google Scholar]
  • Garra R., Salusti E. (2013) Application of the nonlocal Darcy law to the propagation of nonlinear thermoelastic waves in fluid saturated porous media, Phys. D Nonlinear Phenom. 250, 52–57. https://doi.org/10.1016/j.physd.2013.01.014. [CrossRef] [Google Scholar]
  • Gorenflo R., Loutchko J., Luchko Yu. (2002) Computation of the Mittag-Leffler function and its derivatives, Fract. Calc. Appl. Anal. 5, 491–518. [Google Scholar]
  • Gradshteyn S., Ryzhik I.M. (1965) Table of integrals, series, and products, 5th edn, in: Jeffrey A. (eds.), Academic Press, New York. [Google Scholar]
  • Gurtin M.E., Pipkin A.C. (1968) A general theory of heat conduction with finite wave speeds, Arch. Rational Mech. Anal. 31, 2, 113–126. https://doi.org/10.1007/BF00281373. [CrossRef] [MathSciNet] [Google Scholar]
  • Haggerty R., Gorelick S.M. (1995) Multiple-rate mass transfer for modeling diffusion and surface reactions in media with pore-scale heterogeneity, Water Resour. Res. 31, 10, 2383–2400. https://doi.org/10.1029/95WR10583. [Google Scholar]
  • Haubold H.J., Mathai A.M., Saxena R.K. (2011) Mittag-Leffler functions and their applications, J. Appl. Math. 2011, 298628. https://doi.org/10.1155/2011/298628. [Google Scholar]
  • Henry B.I., Langlands T.A.M., Straka P. (2010) An Introduction to Fractional Diffusion, in: Dewar R.L., Detering F. (eds), Complex Physical, Biophysical and Econophysical Systems, World Scientific, Hackensack, NJ, p. 400. [Google Scholar]
  • Hilfer R., Anton L. (1995) Fractional master equations and fractal time random walks, Phys. Rev. E 51, 2, R848–R851. [Google Scholar]
  • Holy R.W., Ozkan E. (2016) A Practical and Rigorous Approach for Production Data Analysis in Unconventional Wells, in: Paper 181662-MS presented at the SPE Low Perm Symposium, 5–6 May, Denver, Colorado, USA. [Google Scholar]
  • Houze O.P., Horne R.N., Ramey H.J. (1988) Pressure-transient response of an infinite-conductivity vertical fracture in a reservoir with double-porosity behavior, SPE Form. Eval. 3, 3, 510–518. https://doi.org/10.2118/12778-PA. [CrossRef] [Google Scholar]
  • Jourde H., Pistrea S., Perrochet P., Droguea C. (2002) Origin of fractional flow dimension to a partially penetrating well in stratified fractured reservoirs. New results based on the study of synthetic fracture networks, Adv. Water Resour. 25, 4, 371–387. [Google Scholar]
  • Kang P.K., Le Borgne T., Dentz M., Bour O., Juanes R. (2015) Impact of velocity correlation and distribution on transport in fractured media: Field evidence and theoretical model, Water Resour. Res. 51, 2, 940–959. https://doi.org/10.1002/2014WR015799. [Google Scholar]
  • Karimi-Fard M., Durlofsky L.J., Aziz K. (2004) An efficient discrete-fracture model applicable for general-purpose reservoir simulators, SPE J. 9, 2, 227–236. https://doi.org/10.2118/88812-PA. [CrossRef] [Google Scholar]
  • Kazemi H. (1969) Pressure transient analysis of naturally fractured reservoirs, Trans. AIME 256, 451–461. [Google Scholar]
  • Kenkre V.M., Montroll E.W., Shlesinger M.F. (1973) Generalized master equations for continuous-time random walks, J. Stat. Phys. 9, 1, 45–50. [Google Scholar]
  • Kim S., Kavvas M.L., Ercan A. (2015) Fractional ensemble average governing equations of transport by time-space nonstationary stochastic fractional advective velocity and fractional dispersion, II: Numerical investigation, J. Hydrol. Eng. 20, 2, 04014040. [Google Scholar]
  • Klafter J., Sokolov I.M. (2011) First steps in random walks, Oxford University Press, p. 152. [Google Scholar]
  • Larsen L., Hegre T.M. (1991) Pressure-transient behavior of horizontal wells with finite-conductivity vertical fractures, in: Paper SPE 22076, Presented at the International Arctic Technology Conference, May 29–31, Anchorage, Alaska, USA. [Google Scholar]
  • Le Mẽhautẽ A., Crepy G. (1983) Introduction to transfer and motion in fractal media: The geometry of kinetics, Solid State Ion. 1, 9–10, 17–30. [Google Scholar]
  • Magin R.L., Ingo C., Colon-Perez L., Triplett W., Mareci T.H. (2013) Characterization of anomalous diffusion in porous biological tissues using fractional order derivatives and entropy, Micropor. Mesopor. Mater. 178, 15, 39–43. https://doi.org/10.1016/j.micromeso.2013.02.054. [CrossRef] [PubMed] [Google Scholar]
  • Mainardi F. (2010) Fractional calculus and waves in linear viscoelasticity, Imperial College Press, London, p. 344. [Google Scholar]
  • Metzler R., Glockle W.G., Nonnenmacher T.F. (1994) Fractional model equation for anomalous diffusion, Phys. A 211, 1, 13–24. [CrossRef] [Google Scholar]
  • Metzler R., Chechkin A.V., Goncharb V.Yu., Klafter J. (2007) Some fundamental aspects of Lévy flights, Chaos Soliton. Fract. 34, 1, 129–142. [CrossRef] [Google Scholar]
  • Metzler R., Jeon J.-H., Cherstvy A.G., Barkai E. (2014) Anomalous diffusion models and their properties: non-stationarity, non-ergodicity, and ageing at the centenary of single particle tracking, Phys. Chem. Chem. Phys. 16, 24128–24164. [CrossRef] [PubMed] [Google Scholar]
  • Mitchell T.M., Faulkner D.R. (2009) The nature and origin of off-fault damage surrounding strike-slip fault zones with a wide range of displacements: A field study from the Atacama fault system, northern Chile, J. Struct. Geol. 31, 8, 802–816. https://doi.org/10.1016/j.jsg.2009.05.002. [Google Scholar]
  • Molz F.J. III, Fix G.J. III, Lu S.S. (2002) A physical interpretation for the fractional derivative in Levy diffusion, Appl. Math. Lett. 15, 7, 907–911. [Google Scholar]
  • Montroll E.W., Weiss G.H. (1965) Random walks on lattices II, J. Math. Phys. 6, 167–181. [Google Scholar]
  • Moodie T.B., Tait R. (1983) On thermal transients with finite wave speeds, J. Acta Mech. 50, 1–2, 97–104. https://doi.org/10.1007/BF01170443. [CrossRef] [Google Scholar]
  • Nigmatullin R.R. (1984) To the theoretical explanation of the universal response, Phys. Status Solidi B Basic Res. 123, 2, 739–745. [CrossRef] [Google Scholar]
  • Nigmatullin R.R. (1986) The realization of the generalized transfer equation in a medium with fractal geometry, Phys. Status Solidi B Basic Res. 133, 1, 425–430. [CrossRef] [Google Scholar]
  • Noetinger B. (2015) A quasi steady state method for solving transient Darcy flow in complex 3D fractured networks accounting for matrix to fracture flow, J. Comput. Phys. 283, 205–223. [Google Scholar]
  • Noetinger B., Estebenet T. (2000) Up-scaling of double porosity fractured media using continuous-time random walks methods, Transp. Porous Med. 39, 3, 315–337. [CrossRef] [Google Scholar]
  • Noetinger B., Estebenet T., Landereau P. (2001) A direct determination of the transient exchange term of fractured media using a continuous time random walk method, Transp. Porous Med. 44, 3, 539–557. https://doi.org/10.1023/A:1010647108341. [CrossRef] [Google Scholar]
  • Noetinger B., Roubinet D., Russian A., Le Borgne T., Delay F., Dentz M., Gouze P. (2016) Random walk methods for modeling hydrodynamic transport in porous and fractured media from pore to reservoir scale, Transp. Porous Med. 115, 2, 345–385. https://doi.org/10.1007/s11242-016-0693-z. [CrossRef] [Google Scholar]
  • Norwood F.R. (1972) Transient thermal waves in the general theory of heat conduction with finite wave speeds, ASME. J. Appl. Mech. 39, 3, 673–676. https://doi.org/10.1115/1.3422771. [CrossRef] [Google Scholar]
  • O’Shaughnessy B., Procaccia I. (1985a) Analytical solutions for diffusion on fractal objects, Phys. Rev. Lett. 54, 5, 455–458. [Google Scholar]
  • O’Shaughnessy B., Procaccia I. (1985b) Diffusion on fractals, Phys. Rev. A 32, 5, 3073–3083. [Google Scholar]
  • Ozkan E., Raghavan R. (1991) Some new solutions to solve problems in well test analysis: I-analytical considerations, SPE Form. Eval. 3, 359–368. https://doi.org/10.2118/18615-PA. [CrossRef] [Google Scholar]
  • Palacio J.C., Blasingame T.A. (1993) Decline-curve analysis with type curves – analysis of gas well production data, Presented at the SPE Low Permeability Reservoirs Symposium, 26–28 April, Denver. https://doi.org/10.2118/25909-MS. [Google Scholar]
  • Povstenko Y. (2015) Linear fractional diffusion-wave equation for scientists and engineers, Birkhäuser, p. 460. [Google Scholar]
  • Prats M. (1961) Effect of vertical fractures on reservoir behavior – incompressible fluid case, Soc. Pet. Eng. J. 1, 2, 105–118. https://doi.org/10.2118/1575-G. [CrossRef] [Google Scholar]
  • Pruess K., Narasimhan T. (1985) A practical method for modeling fluid and heat flow in fractured porous media, SPE J. 25(1), 14–26. https://doi.org/10.2118/10509-PA. [Google Scholar]
  • Raghavan R. (1980) The effect of producing time on type curve analysis, J. Petrol. Technol. 32, 6, 1053–1064. https://doi.org/10.2118/6997-PA. [CrossRef] [Google Scholar]
  • Raghavan R. (2004) A review of applications to constrain pumping test responses to improve on geological description and uncertainty, Rev. Geophys. 42, RG4001. https://doi.org/10.1029/2003RG000142. [Google Scholar]
  • Raghavan R. (2008) A note on the drawdown, diffusive behavior of fractured rocks, Water Resour. Res. 45, 2, W02502. [Google Scholar]
  • Raghavan R. (2011) Fractional derivatives: Application to transient flow, J. Pet. Sci. Eng. 801, 7–13. https://doi.org/10.1016/j.petrol.2011.10.003. [Google Scholar]
  • Raghavan R., Ohaeri C.U. (1981) Unsteady Flow to a Well Produced at Constant Pressure in a Fractured Reservoir, in: Paper SPE-9902MS, Presented at the SPE California Regional Meeting, 25–27 March, Bakersfield, California. https://doi.org/10.2118/9902-MS. [Google Scholar]
  • Raghavan R., Chen C. (2017) Addressing the influence of a heterogeneous matrix on well performance in fractured rocks, Transp. Porous Med. 117, 1, 69–102. https://doi.org/10.1007/s11242-017-0820-5. [CrossRef] [Google Scholar]
  • Raghavan R., Chen C. (2018a) Time and space fractional diffusion in finite systems, Transp. Porous Med. 1231, 173–193. https://doi.org/10.1007/s11242-018-1031-4. [CrossRef] [Google Scholar]
  • Raghavan R., Chen C. (2018b) A conceptual structure to evaluate wells producing fractured rocks of the Permian basin, in: Paper SPE-191484-MS, Presented at the Annual Technical Conference and Exhibition, 24–28 September, Dallas, TX, USA. [Google Scholar]
  • Raghavan R., Chen C. (2019) Evaluation of Fractured Rocks through Anomalous Diffusion, in: Paper SPE-195305-MS, Presented at the SPE Western Regional Meeting, 23–26 April, San Jose, California, USA. [Google Scholar]
  • Raghavan R., Chen C., Agarwal B. (1997) An analysis of horizontal wells intercepted by multiple fractures, SPE J. 2, 3, 235–245. https://doi.org/10.2118/27652-PA. [CrossRef] [Google Scholar]
  • Raghavan R., Dixon T.N., Phan V.Q., Robinson S.W. (2001) Integration of geology, geophysics, and numerical simulation in the interpretation of a well test in a fluvial reservoir, SPE Reserv. Eval. Eng. Soc. Petrol. Eng. 4, 3, 201–208. https://doi.org/10.2118/72097-PA. [CrossRef] [Google Scholar]
  • Savage H.M., Brodsky E.E. (2011) Collateral damage: Evolution with displacement of fracture distribution and secondary fault strands in fault damage zones, J. Geophys. Res. 116, B03405. https://doi.org/10.1029/2010JB007665. [Google Scholar]
  • Saxena R.K., Mathai A.M., Haubold H.J. (2006) Fractional reaction-diffusion equations, Astrophys. Space Sci. 305, 3, 289–296. [Google Scholar]
  • Scholz C.H., Dawers N.H., Yu J.Z., Anders M.H., Cowie P.A. (1993) Fault growth and fault scaling laws: Preliminary results, J. Geophys. Res. 98, B12, 21951–21961. [Google Scholar]
  • Scott K.D., Chu W.-C., Flumerfelt R.W. (2015) Application of real-time bottom-hole pressure to improve field development strategies in the Midland Basin Wolfcamp Shale, in: Paper URTEC-2154675, Proceedings of Unconventional Resources Technology Conference, San Antonio, Texas. https://doi.org/10.15530/URTEC-2015-2154675. [Google Scholar]
  • Spath J.B., Ozkan E., Raghavan R. (1994) An efficient algorithm for computation of well responses in commingled reservoirs, SPE Form. Eval. 9, 2, 115–121. https://doi.org/10.2118/21550-PA. [CrossRef] [Google Scholar]
  • Stehfest H. (1970a) Algorithm 368: Numerical inversion of Laplace transforms [D5], Commun. ACM 13, 1, 47–49. [Google Scholar]
  • Stehfest H. (1970b) Remark on algorithm 368: Numerical inversion of Laplace transforms, Commun. ACM 13, 10, 624. [Google Scholar]
  • Su N. (2014) Mass-time and space-time fractional partial differential equations of water movement in soils: Theoretical framework and application to infiltration, J. Hydrol. 519, B, 1792–1803. [CrossRef] [Google Scholar]
  • Su N., Nelson P.N., Connor S. (2015) The distributed-order fractional diffusion-wave equation of groundwater flow: Theory and application to pumping and slug tests, J. Hydrol. 529, 1262–1273. https://doi.org/10.1016/j.jhydrol.2015.09.033. [CrossRef] [Google Scholar]
  • Suzuki A., Hashida T., Li K., Horne R.N. (2016) Experimental tests of truncated diffusion in fault damage zones, Water Resour. Res. 52, 8578–8589. https://doi.org/10.1002/2016WR019017. [Google Scholar]
  • Tao S., Gao X., Li C., Zeng J., Zhang X., Yang C., Zhang J., Gong Y. (2016) The experimental modeling of gas percolation mechanisms in a coal-measure tight sandstone reservoir: A case study on the coal-measure tight sandstone gas in the Upper Triassic Xujiahe Formation, Sichuan Basin, China, J. Nat. Gas Geosci. 1, 6, 445–455. https://doi.org/10.1016/j.jnggs.2016.11.009. [CrossRef] [Google Scholar]
  • Theis C.V. (1935) The relation between the lowering of the piezometric surface and the rate and duration of discharge of a well using ground-water storage. American Geophysical Union Transactions 16, 519–524. [CrossRef] [Google Scholar]
  • Thomas O.O., Raghavan R., Dixon T.N. (2005) Effect of scaleup and aggregation on the use of well tests to identify geological properties, SPE Res. Eval. Eng. 8, 3, 248–254. https://doi.org/10.2118/77452-PA. [CrossRef] [Google Scholar]
  • Uchaikin V.V. (2013) Fractional derivatives for physicists and engineers, Volume I: Background and theory, Springer, New York, p. 384. [Google Scholar]
  • van Everdingen A.F., Hurst W. (1949) The application of the LaPlace transformation to flow problems in reservoirs, Trans. AIME 186, 305–324. [Google Scholar]
  • Warren J.E., Root P.J. (1963) The Behavior of Naturally Fractured Reservoirs, Soc. Pet. Eng. J. 3, 3, 245–255. https://doi.org/10.2118/426-PA. [CrossRef] [Google Scholar]
  • Yanga S., Zhoua H.W., Zhang S.Q., Ren W.G. (2019) A fractional derivative perspective on transient pulse test for determining the permeability of rocks, Int. J. Rock Mech. Min. Sci. 113, 92–98. [CrossRef] [Google Scholar]
  • Yeh N.S., Davison M.J., Raghavan R. (1986) Fractured well responses in heterogeneous systems-application to Devonian Shale and Austin Chalk reservoirs, ASME. J. Energy Resour. Technol. 108, 2, 120–130. https://doi.org/10.1115/1.3231251. [CrossRef] [Google Scholar]
  • Yeh T.-C.J., Mao D., Zha Y., Wen J., Wan L., Hsu K., Lee C. (2015) Uniqueness, scale, and resolution issues in groundwater model parameter identification, Water Sci. Eng. 8, 3, 175–194. https://doi.org/10.1016/j.wse.2015.08.002. [CrossRef] [Google Scholar]
  • Zhang Y., Benson D.A., Reeves D.M. (2009) Time and space nonlocalities underlying fractional-derivative models: Distinction and literature review of field applications, Adv. Water Res. 32, 561–581. [CrossRef] [Google Scholar]
  • Zhokh A., Strizhak P. (2018) Non-Fickian transport in porous media: Always temporally anomalous? Transp. Porous Med. 124, 2, 309–323. https://doi.org/10.1007/s11242-018-1066-6. [CrossRef] [Google Scholar]
  • Zhokh A., Strizhak P. (2019) Investigation of the anomalous diffusion in the porous media: a spatiotemporal scaling, Heat Mass Trans. 55, 1–10. https://doi.org/10.1007/s00231-019-02602-4. [CrossRef] [Google Scholar]

All Tables

Table 1

Flow regimesa,b,c,d.

Table 2

Exponents of flow regimes.

All Figures

thumbnail Fig. 1

Rate-normalized responses at five wells producing the same lease in the Eagle Ford Group; after Doe et al. (2014).

In the text
thumbnail Fig. 2

Schematic of the transient, interporosity model of Kazemi (1969) and de Swaan-O (1976); after Raghavan and Ohaeri (1981).

In the text
thumbnail Fig. 3

Schematic: alignment of hydraulic fractures with wellbore.

In the text
thumbnail Fig. 4

A fractured well in a rectangular drainage region. The fracture with identical wing lengths with its center at (xw, yw) is located anywhere in the drainage region and is parallel to one of the boundaries. The length, width and height of the fracture are xF, wF, and h, respectively.

In the text
thumbnail Fig. 5

Production response at a horizontal well produced through multiple fractures at a constant rate to highlight Flow Regimes 1 and 3; the fracture conductivity is assumed to be infinitely large, CFD ≥ 500. The unbroken lines are the pressure, pwD, and derivative, p wD $ {p}_{{wD}}^\mathrm{\prime}$, responses respectively. Asymptotic solutions corresponding to short and long times, Flow Regime 1 and Flow Regime 3, respectively, are shown as dashed lines. The circles are solutions for a single fracture presented in the manner indicated in equation (58).

In the text
thumbnail Fig. 6

Production response at a horizontal well produced through multiple fractures at a constant rate to highlight Flow Regime 2; two values of fracture conductivity are considered, CFD of 1 and 5. The unbroken lines are the pressure, pwD, and derivative, p wD $ {p}_{{wD}}^\mathrm{\prime}$, responses respectively. The dashed lines reflect the asymptotic solution, Flow Regime 2. Values of λ′ and ω′ are so chosen to facilitate fluid transfer and promote the duration of Flow Regime 2. The circles are solutions for a single fracture presented in the manner indicated by equation (58).

In the text
thumbnail Fig. 7

Buildup (circles) and flowing well (squares) responses for a shale, oil well.

In the text
thumbnail Fig. 8

Green’s function responses corresponding to the test discussed in Figure 7; subdiffusive responses correspond to Flow Regime 4 and combined responses to Flow Regime 2.

In the text
thumbnail Fig. 9

Correlation of x D G ( x D , t D ) $ {x}_D\mathcal{G}({x}_D,{t}_D)$ as a function of the similarity variable, ς, of plane source solutions. Values of xD and tD in the range 0.2 ≤ xD ≤ 5 and 3 × 10−2 ≤ tD ≤ 500 are considered.

In the text

Les statistiques affichées correspondent au cumul d'une part des vues des résumés de l'article et d'autre part des vues et téléchargements de l'article plein-texte (PDF, Full-HTML, ePub... selon les formats disponibles) sur la platefome Vision4Press.

Les statistiques sont disponibles avec un délai de 48 à 96 heures et sont mises à jour quotidiennement en semaine.

Le chargement des statistiques peut être long.