Regular Article
A study in fractional diffusion: Fractured rocks produced through horizontal wells with multiple, hydraulic fractures
^{1}
Raghavan Group, PO Box 52756, Tulsa, OK 74152, USA
^{2}
Kappa Engineering, 11767 Katy Freeway, #500, Houston, TX 77079, USA
^{*} Corresponding author: raghavan.raj@gmail.com
Received:
11
June
2020
Accepted:
28
July
2020
The spatiotemporal evolution of transients in fractured rocks often displays unusual characteristics and is traced to multifaceted origins such as microdiscontinuity in rock properties, rock fragmentation, longrange connectivity and complex flow paths. A physicsbased model that incorporates transient propagation wherein the mean square displacement of the diffusion front follows a nonlinear scaling with time is proposed. This model is based on fractional diffusion. The motivation for fractional flux laws follows from the existence of longrange connectivity that results in the mean square displacement of fronts moving faster than predicted by classical models; correspondingly, obstructions and discontinuities, local flow reversals, intercalations, etc. produce the opposite effect with fronts moving at a slower rate than classical predictions. The interplay of these two competing behaviors is quantified. We simulate transient production in a porous rock at the Theis scale as a result of production through a horizontal well consisting of multiple hydraulic fractures. Asymptotic solutions are derived and computations verified. The practical potential of this model is described through an example and the movement of fronts under the constraints of this model is demonstrated through the new expressions developed in this work. We demonstrate that this model offers a potential avenue to explain other behaviors noted in the literature. Though this work is developed in the context of applications to the earth sciences (production of hydrocarbons, extraction of geothermal resources, sequestration of radioactive waste and other fluids, groundwater flow), a minimal change in the Nomenclature permits application to other contexts. Ideas proposed here are particularly useful in the context of superdiffusion in bounded systems which until now, in many ways, has been considered to be an open problem.
© R. Raghavan & C.C. Chen, published by IFP Energies nouvelles, 2020
This 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)
: See equation (A.41)
C_{FD} : Dimensionless fracture conductivity
B : Formation volume factor [L^{3}/L^{3}]
E_{α,β}(z): MittagLeffler function; see equation (24)
: See equation (15)
: See equation (21)
IRNP(t_{e}): Integral of ratenormalized pressure
k_{α,β} : See equation (3)
L : Width of linear system [L]
M_{ν}(z): Mainardi or MWright function
RNP(t_{e}): Ratenormalized pressure
t_{e} : Material balance time [T]
x_{L} : Distance to interface [L]
β_{f} : See equation (A.5)
β ^{⋆} : See equation (A.18)
β ^{⋆⋆} : See equation (A.29)
: Equation (18) []
λ′: See equation (19)
ς : Similarity variable; see equation (64)
ω′: See equation (20)
Subscripts
Superscript
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 welldefined, 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; KarimiFard 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 equivalentporousmedia 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 timedistance 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 welltests). 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: with the principal conclusion suggesting that at long enough times transients display characteristic powerlaw features, namely, Δp ~ t ^{ ϑ }, where ϑ is a constant, contrary to the widely known result Δp ~ lnt. The finding that Δp ~ Γ(a, x) 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 powerlaw 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 superdiffusion 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). Earlytime responses, as is usual in most wells producing unconventional reservoirs, reflect the lack of appropriate measurements particularly bottomhole 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.
Fig. 1 Ratenormalized 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 reservoirengineering 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, X_{n}; the first moment (the conservation of mass principle) and second moment (the timedistance 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 crosssectional 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 for 0 < β < 1^{1} with β reflecting the power law tail (); see Benson et al. (2000, 2004). Levystableparticles trace a path of dimension β + 1 namely the exponent of the xderivative, ∂_{x}, in ; 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 for 0 < α < 1 rather than ; in this case for 0 < α < 1, where T_{n} 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 waitingtime solutions as documented in Henry et al. (2010). In fact, the expression 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), d_{w} (Metzler et al., 1994) for diffusion on fractal objects. The earthscience 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 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 pointsource 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 SwaanO (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 SwaanO (1976) which shows the colligation of the matrix and fracture elements where h_{f} denotes the height of a fracture element and h_{m} denotes the height of a matrix element. As we consider n elements, the total thickness, h_{t}, is
Fig. 2 Schematic of the transient, interporosity model of Kazemi (1969) and de SwaanO (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} (CincoLey 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, fractionalflux 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 fluxlaw 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 waitingtime 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)
(2)and assumes to govern transients in the fractured rock by an expression of the form for the velocity, v(x, t):
(3)with the exponents α and β both reflecting the heterogeneities in the fractured rock.^{3} Note that unlike most studies a onesided 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:
Here p_{i} the initial pressure is identical in all parts of the system at t = 0, and Δp_{j} = p_{i} − p_{j}(x, t); j being f, m or F, , . We require a solution to equations (4)–(6) subject to the following conditions. For a quiescent system where equilibrium prevails everywhere we require
For the matrix elements we require that the boundary at z = 0 be a noflow boundary; thus
Further, for perfect contact we require
(9)as the requirement of continuity in pressure at the fracturematrix interface, Ω_{f} − Ω_{m}, that is, at (z = h_{m}/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 righthand side of the equation.
The constraints for the pressure distribution are as follows. All boundaries are assumed to be sealed. Hence we have
The conditions to be satisfied at the fracturereservoir interface, Ω_{f} – Ω_{F}, are, again, the continuity of pressure and the continuity of flux; in consequence,
(13)respectively, where the symbol represents the flux.
As is usual we follow the outline of the idealization given in CincoLey 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, k_{F}, and finite width w_{F}. For large values of the product k_{F} w_{F} 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
(18)with 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
The solution of equation (14) and its associated boundary conditions may then be coupled with the solution of equation (6) given by
(22)where , and 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 and 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, Δp_{w}(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, Δp_{w}. 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).
Fig. 3 Schematic: alignment of hydraulic fractures with wellbore. 
2.2.2 The function
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 SwaanO (1976), Chen (1982), Ozkan and Raghavan (1991), CincoLey and Meng (1988), Noetinger et al. (2001, 2016) and Raghavan and Chen (2017). Distinctive features of well responses are governed by the variation of with respect to s (time). Three options are available: For s→ ∞; early times specified as Flow Regime 1, ( and tanh(x) ≈ 1; matrix has no influence), for finite values of s such that tanh(x) ≈ 1 and ; intermediate times specified as Flow Regime 2, (; flow regime in the matrix system reflects transient conditions), and for s → 0; late times specified as Flow Regime 3, ( 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 transientflowmodel, 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).
Flow regimes^{a}^{,}^{b}^{,}^{c}^{,}^{d}.
In concluding this section, we should note that additional options exist to examine solutions; for example, Chen (1982) considers ; 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 planesourcesolution 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 planesource solution is particularly useful in obtaining early time responses for a fractured well. Thus, considering a linear system we may derive the planesource solution along the lines outlined in Carslaw and Jaeger (1959). The required expression is given by
In the above expression, E_{α,β}(z) is the twoparameter MittagLeffler function (Erdelyi et al., 1955; Mainardi, 2010; Uchaikin, 2013; and Povstenko, 2015):
(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: E_{1}(z) = exp(z), E_{2,1}(z ^{2}) = cosh(z), E_{2,2}(z ^{2}) = sinh(z)/z, and exp(z ^{2})erfc(−z) = E_{1/2,1}(z). Uchaikin (2013) provides convenient tabulations of expressions for obtaining inverse Laplace transforms containing E_{α,β}(z).
For β_{f} = 1, the righthand side of equation (23) reduces to an expression in terms of a stretched exponential given by
The righthand side of equation (25) may be inverted either through the Fox function or Hfunction, , by using the relation
(27)where , ς > 0 as noted in Saxena et al. (2006), or through the Wright function, W(α, β; t), as indicated in Povstenko (2015) by the expression
The Wright function is defined by
2.3.2 Probability density function
If f(s) = 1, the righthand side of equation (23) is a probability density function. Denoting to be the probability density function, we require that
To get the integral of the righthand side of equation (23), we use
(32)and thus on integration of equation (23) the result is:
On using the asymptotic expansion in Haubold et al. (2011),
(34)for z ϵ C and ξ is a real number with β < 1, we have
Therefore, equation (31) is satisfied. Consequently, the righthand side of equation (23) is a pdf.
If both α_{f} and β_{f} equal to 1, we have the Gaussian as the inversion of the righthand side of in equation (23) yields the expression for γ(x, t) to be
The moments of an even order which would also be of interest are (Uchaikin, 2013). For β = 1,
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 is
(38)with the symbol S in the case under study representing the surface of the source. If were constant, then equation (38) simplifies to
For the system under consideration, equations (23) and (38) yield the following expression for pressure distribution:
(40)where is the distance from the source. The flux distribution in the system corresponding to equation (40) is
The above expressions are also given in Chen and Raghavan (2015).
For a bounded system, the pressure distribution, , in terms of the Laplace transformation may be obtained from
We may relate the fractional gradient to the pressure drop at by the expression
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 q_{j} at each Fracture j, we may write
The assumption of an infiniteconductivity 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
The convolution equation to compute the pressure response at the wellbore for the multiplefracture system is
(48)for k = 1, 2, …, N. Here, jk represents the effect of Fracture j on the pressure response at Fracture k and .
In terms of the Laplace transformation, the above three expressions, respectively, are as follows:
Thus we are required to solve simultaneously for Δp_{w} and q_{j}; j = 1, 2, …, N, namely the system
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 offdiagonal elements of the matrix would have reflected interactions among fractures. Each 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 loglog 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 nonfissured idealizations. In fact, Table 2 reproduces all information available in the literature.
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
For subdiffusive flow; that is, for β_{f} = 1, equation (53) yields the expression
(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 would exist for situations wherein the fracture length and the width of the rectangle are not identical as indicated in Figure 4 (x_{e}/x_{f} > 1); see Appendix.
Fig. 4 A fractured well in a rectangular drainage region. The fracture with identical wing lengths with its center at (x_{w}, y_{w}) is located anywhere in the drainage region and is parallel to one of the boundaries. The length, width and height of the fracture are x_{F}, w_{F}, and h, respectively. 
In deriving an asymptotic solution for Flow Regimes 4 and 5, we show in the Appendix that twoterm 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 timedistance 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 MittagLeffler 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 loglog coordinates) will be even much smaller. Chu (2018), who framed his comments in terms of the Chow (1952) Pressure Group, Δp_{c}, has noted that we need to consider situations when slopes are very, very small implying situations beyond subdiffusive effects. And superconnected networks resulting from β_{f} < 1 may result in slopes (exponents of t) of lines on loglog 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 , 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 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, p_{wD}(t_{D}), dimensionless time, t_{D}, and dimensionless fracture conductivity, C_{FD}, are defined, respectively, as
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, C_{FD} ≥ 500. The unbroken lines are the pressure, p_{wD}, and derivative, , 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). 
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, C_{FD} of 1 and 5. The unbroken lines are the pressure, p_{wD}, and derivative, , 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, p_{wD}(t_{D}), and its derivative , (tdp_{w/}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
(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 C_{FD}, 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 singlefracture 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 filledin circles being buildup responses (Δp vs. Δt) and the squares representing flowing responses IRNP(t_{e}) vs. t_{e}, in terms of the nomenclature used in RateTransientAnalysis, RTA. The symbol IRNP(t_{e}) is the integral of ratenormalized pressure defined in Palacio and Blasingame (1993) as
(59)and the symbol t_{e} represents materialbalance time defined in CamachoVelázquez and Raghavan (1989) as
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 earlytime flowingwell 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 wellfounded basis for a model of the kind proffered here. The value of α_{m} along with the estimate of 0.1 for the earlytime 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 longrange 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
(61)with α_{m} = 0.6. Differences between the bottom two curves are to only be expected; all that is indicated is the wellknown point that many combinations of properties yield the same pressure response.
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, ς,
(63)where ς is deduced by considering the arguments of the MittagLeffler, Fox and Wright functions on the righthand 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 as function of ς we have constructed is shown for several values of x_{D} and t_{D}. The circles, squares, diamonds, etc. represent specific vales of x_{D} as noted in Figure 9. Curves for distinct values of locations, x_{D}, may be correlated in this manner for specific combinations of α and β.
Fig. 9 Correlation of as a function of the similarity variable, ς, of plane source solutions. Values of x_{D} and t_{D} in the range 0.2 ≤ x_{D} ≤ 5 and 3 × 10^{−2} ≤ t_{D} ≤ 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
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 flowpaths 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 finiteconductivity model for a single fracture in a rectangular drainage region and then sketch out the details of the multiplefracture 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
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, , to be
(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
Thus, on defining
(A.8)we obtain the well pressure in terms of the Laplace transformation to be
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 constantpressureproduction 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
(A.10)where the symbol is now the righthand 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
Using this result we consider the variation of the second tanh(x) of the solution developed above, namely in equation (A.6) to illustrate responses for the Bilinear and Linear flow regimes.
A.1.2 Bilinear flow
For time ranges when 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
On substituting the righthand side of equation (A.4) for after noting the dependence on , equation (A.13) is
On extracting the first term from the parentheses on the righthand side of equation (A.15), expanding the resulting expression by the binomial theorem and dropping all but the leading terms we replace in equation (A.10) by the final result to obtain
Now on using the expression for in equation (A.12), we obtain the dimensionless pressure response, p_{wD}(t_{D}), in terms of dimensionless time, t_{D}, and dimensionless fracture conductivity, , to be
The symbols p_{wD}(t_{D}), t_{D}, and , are defined, respectively, by the following expressions
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 or tanh(x) ≈ x − x ^{3}/3 for (CamachoVelázquez, 1984; CamachoVelázquez et al., 1987). We follow the latter option, and using the righthand side of equation (A.4) for we may write
(A.25)with and a, respectively, given as
With the result noted in equation (A.22) we may obtain well responses for linear flow, again, on substituting appropriate expressions for and we do so on the assumption that . On using the righthand side of equation (A.22) in equation (A.10) for , along with the expression in equation (A.12) for , simplifying and inverting the result, the dimensionless pressure, p_{wD}(t_{D}), during linear flow under the constraint that Flow Regime prevails is given by
(A.30)where L_{Dj} 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
A.1.4 Flow Regime 3
In this section we address our expectations of well responses for the linearflow 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
(A.32)consequently, may be written as
With the expression in equation (A.22), we may write
Consequently, the well response, p_{wD}(t_{D}), is given by
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 .
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 MittagLeffler functions, for small values of s, the ratio of MittagLeffler Functions, and , may be written as
Thus the approximation for φ in equation (44) is then given by
On again using a twoterm approximation for the E_{α,β}(z) expressions, we may write equation (43) with the aid of the expression in equation (A.40) to be
In the above expression, the flow rate, . In order to relate the solutions to subdiffusion, we now modify the solution in equation (A.42), replace 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,
The expression for dimensionless pressure is thus given by
The analog of equation (54) for constantrate production is
Using the result given on page 47 (1.445,2) of Gradshteyn and Ryzhik (1965), namely,
(A.47)we recast equation (A.46) as
If we consider the solution given in equation (A.48) as , and use the result of Gradshteyn and Ryzhik (1965) given on page 46 (1.445,2) namely that
(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 (x_{e}/x_{f} ≠ 1) as shown in Figure 4 is given as
(A.50)where is the flux distribution along the fracture face, , , and , and x_{F} = x_{F1} + x_{F2}.
In equation (A.50), at long enough times, the well response for x_{e}/x_{f} 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 MittagLeffler 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 in equation (A.44). Inversion of the resulting expressions yields the pressure distributions during Flow Regimes 4 and 5 which are, respectively,
Turning our attention to production rate, q(t), for situations wherein the well produces at a constant pressure; that is, when at in equation (A.42) as given in van Everdingen and Hurst (1949), we may readily show that for Flow Regime 4, , where a is a constant or (q(t) ~ 1/t ^{1/2}) if α_{m} = 1) at terminal conditions because at long enough times E_{α,β}(−x) ~ 1/x for α ≠ β. Such a powerlaw trend is reflective of disconnectedness, disjointedness, ruptures, and dead ends as a result of intercalations and the like. Similarly, for Flow Regime 5, or a powerlaw trend of the form 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
(A.53)where q_{wDk}(t) is the flow rate for a well produced at a constant pressure, p_{wf}, through a fracture with the properties of Fracture k. If we restrict our attention to the boundarydominated period, then equation (A.42) yields the following expression for determining the production rate, q_{wp}(t), at Fracture k:
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.
The fractional flux law for space fractional diffusion may be stated in other forms. Often a pde in the form of with 0 < β < 2 is considered to model superdiffusion. In such cases a fractional flux law of the form 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? RueilMalmaison, 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., RuizMedina 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 URTEC20203087MS, 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) Pressuretransient model for a vertically fractured well in a fractal reservoir, SPE Form. Eval. 9, 2, 122–128. https://doi.org/10.2118/20582PA. [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/17422132/3/3/003. [CrossRef] [Google Scholar]
 Benson D.A., Wheatcraft S.W., Meerschaert M.M. (2000) Application of a fractional advectiondispersion equation, Water Resour. Res. 36, 6, 1403–1412. [Google Scholar]
 Benson D.A., Tadjeran C., Meerschaert M.M., Farnham I., Pohll G. (2004) Radial fractionalorder 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]
 CamachoVelázquez R.G. (1984) Response of wells producing commingled reservoirs: Unequal fracture length, Master’s thesis, University of Tulsa, Tulsa, OK. [Google Scholar]
 CamachoVelá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/12844PA. [CrossRef] [Google Scholar]
 CamachoVelázquez R.G., Raghavan R. (1989) Boundarydominated flow in solutionsgasdrive reservoirs, SPE Reserv. Eng. 4, 4, 503–512. https://doi.org/10.2118/18562PA. [CrossRef] [Google Scholar]
 Caputo M. (1967) Linear models of dissipation whose Q is almost Frequency IndependentII, Geophys. J. Roy. Astron. Soc. 13, 5, 529–539. [NASA ADS] [CrossRef] [Google Scholar]
 Carrera J., SánchezVila 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) Pressuretransient 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 fracturedhorizontal wells and conclusions drawn thereof, SPE Reserv. Eval. Eng. 16, 1, 19–28. https://doi.org/10.2118/163104PA. [CrossRef] [Google Scholar]
 Chen C., Raghavan R. (2015) Transient flow in a linear reservoir for spacetime 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) Ratetransient analysis based on powerlaw behavior for Permian wells, SPE Res. Eval. Eng. 22, 4, 1360–1370 https://doi.org/10.2118/187180PA. [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/191407PA. [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, LBNL59885, Water Resour. Res. 42, W10201. [Google Scholar]
 CincoLey 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/18172MS. [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 SwaanO 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/5346PA. [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 IPR0229, SKB, Stockholm. [Google Scholar]
 Doe T.W. (1991) Fractional dimension analysis of constantpressure well tests, in: Paper SPE22702MS, Presented at the SPE Annual Technical Conference and Exhibition, 6–9 October, Dallas, Texas. https://doi.org/10.2118/22702MS. [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, McGrawHill, 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 MittagLeffler 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) Multiplerate mass transfer for modeling diffusion and surface reactions in media with porescale 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) MittagLeffler 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 181662MS 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) Pressuretransient response of an infiniteconductivity vertical fracture in a reservoir with doubleporosity behavior, SPE Form. Eval. 3, 3, 510–518. https://doi.org/10.2118/12778PA. [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]
 KarimiFard M., Durlofsky L.J., Aziz K. (2004) An efficient discretefracture model applicable for generalpurpose reservoir simulators, SPE J. 9, 2, 227–236. https://doi.org/10.2118/88812PA. [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 continuoustime 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 timespace 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) Pressuretransient behavior of horizontal wells with finiteconductivity 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., ColonPerez 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: nonstationarity, nonergodicity, 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 offfault damage surrounding strikeslip 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) Upscaling of double porosity fractured media using continuoustime 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/s112420160693z. [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: Ianalytical considerations, SPE Form. Eval. 3, 359–368. https://doi.org/10.2118/18615PA. [CrossRef] [Google Scholar]
 Palacio J.C., Blasingame T.A. (1993) Declinecurve 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/25909MS. [Google Scholar]
 Povstenko Y. (2015) Linear fractional diffusionwave 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/1575G. [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/10509PA. [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/6997PA. [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 SPE9902MS, Presented at the SPE California Regional Meeting, 25–27 March, Bakersfield, California. https://doi.org/10.2118/9902MS. [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/s1124201708205. [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/s1124201810314. [CrossRef] [Google Scholar]
 Raghavan R., Chen C. (2018b) A conceptual structure to evaluate wells producing fractured rocks of the Permian basin, in: Paper SPE191484MS, 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 SPE195305MS, 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/27652PA. [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/72097PA. [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 reactiondiffusion 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 realtime bottomhole pressure to improve field development strategies in the Midland Basin Wolfcamp Shale, in: Paper URTEC2154675, Proceedings of Unconventional Resources Technology Conference, San Antonio, Texas. https://doi.org/10.15530/URTEC20152154675. [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/21550PA. [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) Masstime and spacetime 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 distributedorder fractional diffusionwave 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 coalmeasure tight sandstone reservoir: A case study on the coalmeasure 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 groundwater 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/77452PA. [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/426PA. [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 systemsapplication 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 fractionalderivative models: Distinction and literature review of field applications, Adv. Water Res. 32, 561–581. [CrossRef] [Google Scholar]
 Zhokh A., Strizhak P. (2018) NonFickian transport in porous media: Always temporally anomalous? Transp. Porous Med. 124, 2, 309–323. https://doi.org/10.1007/s1124201810666. [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/s00231019026024. [CrossRef] [Google Scholar]
All Tables
All Figures
Fig. 1 Ratenormalized responses at five wells producing the same lease in the Eagle Ford Group; after Doe et al. (2014). 

In the text 
Fig. 2 Schematic of the transient, interporosity model of Kazemi (1969) and de SwaanO (1976); after Raghavan and Ohaeri (1981). 

In the text 
Fig. 3 Schematic: alignment of hydraulic fractures with wellbore. 

In the text 
Fig. 4 A fractured well in a rectangular drainage region. The fracture with identical wing lengths with its center at (x_{w}, y_{w}) is located anywhere in the drainage region and is parallel to one of the boundaries. The length, width and height of the fracture are x_{F}, w_{F}, and h, respectively. 

In the text 
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, C_{FD} ≥ 500. The unbroken lines are the pressure, p_{wD}, and derivative, , 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 
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, C_{FD} of 1 and 5. The unbroken lines are the pressure, p_{wD}, and derivative, , 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 
Fig. 7 Buildup (circles) and flowing well (squares) responses for a shale, oil well. 

In the text 
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 
Fig. 9 Correlation of as a function of the similarity variable, ς, of plane source solutions. Values of x_{D} and t_{D} in the range 0.2 ≤ x_{D} ≤ 5 and 3 × 10^{−2} ≤ t_{D} ≤ 500 are considered. 

In the text 