Regular Article
Spacetime fractional diffusion: transient flow to a line source
^{1}
The Raghavan Group, PO Box 52756, Tulsa, OK 74152, USA
^{2}
Kappa Engineering, 11767 Katy Freeway #500, Houston, TX 77079, USA
^{*} Corresponding author: chen@kappaeng.com
Received:
3
August
2021
Accepted:
27
October
2021
Nonlocal diffusion to a line source well is addressed by spacetime fractional diffusion to model transients governed by both longrange connectivity and distorted flow paths that result in interruptions in the geological medium as a consequence of intercalations, dead ends, etc. The former, superdiffusion, results in longdistance runs and the latter, subdiffusion, in pauses. Both phenomena are quantified through fractional constitutive laws, and two exponents α and β are used to model subdiffusion and superdiffusion, respectively. Consequently, we employ both time and space fractional derivatives. The spatiotemporal evolution of transients in 2D is evaluated numerically and insights on the structure of solutions described through asymptotic solutions are confirmed numerically. Pressure distributions may be classified through two situations (i) wherein 2α = β + 1 in which case solutions may be grouped on the basis of the classical Theis solution, and (ii) wherein 2α ≠ β + 1 in which case conventional expectations do not hold; regardless, at long enough times for the combined case, powerlaw responses are similar to those for pure subdiffusive flows. Pure superdiffusion on the other hand, although we consider a system that is infinite in its areal extent, interestingly, results in behaviors similar to steadystate flow. To our knowledge, documented behaviors are yet to be reported.
© R. Raghavan & C. Chen, published by IFP Energies nouvelles, 2021
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
c : Compressibility [L T^{2}/M]
E_{α,β} (z): Mittag–Leffler function; see (11)
k_{α,β } : Fractional permeability; see (2) [L^{β+1}T^{1−α}]
p′: Logarithmic derivative [M/L/T^{2}]
r_{w} : Radius of well bore [L]
φ_{α,β} (x,t): Green’s function
: Fractional diffusivity; see (6) [L^{β+1}T^{1−α}]
Subscripts
Superscript
1 Introduction
The ubiquity of heterogeneity in geologic media hardly needs an introduction – every occasional gardener knows firsthand. Technical ideas of the concept of heterogeneity preceded the introduction of Darcy’s fluxlaw with the introduction of the definition of facies, we are told, in 1838. Most working models of well performance, in practice, often degenerate to diffusion in homogeneous systems where fluid movement is presumed to follow a velocityjump random walk following the working hypothesis that the mean square displacement or mean square deviation, 〈(Δx)^{2}〉, is proportional to time, t; that is, 〈(Δx)^{2}〉 ∝ t. Over time this idea has been expanded by resorting either to numerical means as summarized in Mattax and Dalton (1990) or analytically by considering Continuous Time Random Walks, CTRW, (Cortis and Knudby, 2006; Noetinger and Estebenet, 2000; Noetinger et al., 2016), fractal objects (Ali et al., 2014; Chang and Yortsos, 1990; FlamencoLopez and CamachoVelazquez, 2001; Metzler et al., 1994; Park et al., 2000; Wang, 2013), waitingtime solutions (Henry et al., 2010) and Levystable random walks (Benson et al., 2000, 2004; Cloot and Botha, 2006; Gehlhausen, 2009; Lenormand, 1992). Many field evaluations of well responses in recent times (Bernard et al., 2006; Chu et al., 2019a, b; Cloot and Botha, 2006; Doe et al., 2014; Meerschaert et al., 2008; Raghavan and Chen, 2019; Scott et al., 2015, Sinkov et al., 2021; Thomas et al., 2005) note that the proper assessment of responses in systems producing with complex geology requires going beyond the classical transient diffusion; that is, beyond the premises of the Theis (1935), solution. They have proposed incorporation of the internal topology of the complex porous network along with actual paths where distant locations may be connected rather well whereas nearby locations are isolated because of localized barriers and dead ends which lead to mean square displacements of the form 〈(Δx)^{2}〉 ∝ t^{a} with . Such an approach also requires reconsideration of the flux law – a new flux law in terms of fractional derivatives; see, for example, Fomin et al. (2011) and Chen and Raghavan (2015).
In this study we explore for the first time, to our knowledge, the combined effects of super and subdiffusion by considering pressure distributions in 2D around a linesource well by examining solutions of the fractional differential equation, for 0 < α, β < 1. The properties of statistical distributions that govern the influence of obstructions, restrictions, impediments, intercalations, faultdamage zones, etc. are reflected through the exponent α, and the connectivity between distant points is reflected through the exponent β. Specifically, α reflects the waiting time in a CTRW (Henry et al., 2010; Noetinger and Estebenet, 2000), for 0 < α < 1, where T_{n} is the random waiting time at which the jump occurs, and β reflects the path taken by Levystableparticles through the power law tail () (Benson et al., 2004). Considered individually, the exponent, α, represents the reciprocal of the anomalous diffusion coefficient (random walk dimension), d_{w}, (Metzler et al., 1994) for diffusion on fractal objects, , and the exponent, β, represents the Hurst index, H, for Levystable diffusion; H is 1/(β + 1) or <r^{2}> ~ t^{H}. Elaborations of these ideas with respect to transient, anomalous diffusion in porous solids are given in Raghavan and Chen (2020). For α = β = 1, the fractional equation reverts to the classical form with alternating sequences of fixedvelocity runs and reorientations. The analytical expressions of the Greens functions in Luchko (2016), Bazhlekova (2019) and Bazhlekova and Bazhlekov (2019) form the cornerstone of this work. The analysis of numerical experiments of expressions for the pressure distribution to a linesource well sheds light on behaviors yet to be reported; unusual features of pure superdiffusion in 2D which one may not expect are reported and differences with respect to 1D systems are documented. The results of the numerical experiments also establish the viability of the solutions for practical use and form the centrepiece of more complex well completions to be reported elsewhere. Recall that we consider both super and subdiffusion behaviors.
2 The differential equation, auxillary conditions
This work, as is usual for singlephase flow, concentrates on the flow of a slightly compressible liquid to a cylindrical well which is treated as a line source; the well is considered to be located in a system that is infinite in extent. The structure considered is described by a set of scaling exponents, α and β, and we work in terms of pressure, p(r, t), permeability, k_{α,β}, porosity, ϕ, thickness, h, fluid compressibility, c, and viscosity, μ in the cylindrical coordinate system assuming the axis of the well passes through the origin. Initially, at time t = 0, the pressure distribution is distributed uniformly at the value, p_{i}; that is, equilibrium prevails at t = 0. Wellbore storage and skin effects are presumed to be zero.
The conservation equation of mass principle to a control volume in cylindrical coordinates for a slightly compressible, constant viscosity fluid is
where v(r, t) is the velocity. The fractional flux law we use in this work is of the form
where λ_{α,β} = k_{α,β}/μ, with k_{α,β} being the effective “fractional permeability” with dimensions of . Both the exponents α and β are <1 with the temporal fractional derivative, ∂^{α}f(t)/∂t^{α}, defined in the Caputo (1967) sense, namely
and the spatial fractional derivative, ∂^{α} f(r)/∂r^{α}, in the Riemann–Liouville form, namely
where .
As is well known fractional derivatives are convolutions; that is, nonlocal operators; for example, in [Oldham and Spanier (1974), Podlubny (1998)] equation (4) the value of the fractional derivative at the point r is a function of all points to the left of r. Expressions of the form given in equations (3) and (4) have been used previously (Fomin et al., 2011; Chen and Raghavan, 2015) but the manner in which fractional derivatives are used here is slightly different. For superdiffusive flows; that is, α = 1, Benson et al. (2004) use an expression similar to that given here. Many discussions of fractional flux laws are available; see Noetinger and Estebenet (2000), Suzuki et al. (2016), Henry et al. (2010), Su (2014), and Su et al. (2015) for subdiffusive flow and Molz et al. (2002), Garra and Salusti (2013), Kim et al. (2015) and Sapora et al. (2017) for superdiffusive flow. A number of experimental works also exist: Tao et al. (2016), Yanga et al. (2019), Zhokh and Strizhak (2018, 2019) and Chang et al. (2019). Fractional flux laws are also needed to address flow in fractal systems; see Le Mehaute and Crepy (1983), Nigmatullin (1984) and Beier (1990). 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). Raghavan and Chen (2019) summarize estimates of α based on a number of field studies for rocks of different types; they report the estimates in the range: 0.56 ≤ α ≤ 0.91 for β = 1. Expressions for the flux distributions are not always addressed as the preference is to use pseudodifferential operators (Balakrishnan, 1960).
On combining the conservation of mass principle and the flux law after noting that we consider flow to a well located at the origin we, as usual, arrive at the differential equation given by:
where , the “fractional diffusivity” of the medium, is given as
We seek a solution to equation (5) subject to the following initial and boundary conditions:
and
where q is the rate. Equation (9) specifies the condition for a linesource well; that is, lim r_{w} → 0.
The analog of an expression corresponding to the lefthand side of equation (5) in 1D and 2D Cartesian coordinates given by the Riesz fractional derivative, , has been studied extensively as in Mainardi et al. (2001) and Huang and Liu (2005) for 1D, and as a fractional Laplacian, −(−Δ)^{α}p(x, t); 0 < α < 1, as in Luchko (2019), Bazhlekova (2019) and Bazhlekova and Bazhlekov (2019) for 2D. The correspondence with classical diffusion is usually established through Fourier transforms (Gorenflo and Mainardi, 1998); it is given in a pointwise sense in Cai and Li (2019) and the fractional flux law becomes Darcy’s law in the limit. Kwaśnicki (2017) notes that 10 equivalent definitions of the fractional Laplacian have been used.
Incidentally,
see EstradaRodriguez et al. (2018).
3 The Mittag–Leffler functions, E_{α} (z) and E_{α,β}(z)
As is inevitable in most analytical studies of anomalous diffusion, this work requires the use of the Mittag–Leffler functions, E_{α} (z) and E_{α,β} (z), which are briefly reviewed. Extensive discussions may be found in Erdelyi et al. (1955), Mainardi (2010), Uchaikin (2012) and Gorenflo et al. (2014) to name a few sources.
The function is defined by the representation
where denotes the complex plane. Implicit in all of the following is that .
Of particular interest to us is the Laplace transformation of . From the references cited, we have
where denotes the Laplace transformation, s is the Laplace transform variable, and both α and β are real. Thus, for β = 1 we have
The following asymptotic expressions for and are often needed: For 0 < α < 2 and with we have
and
We compute the Mittag–Leffler functions by the algorithm in Gorenflo et al. (2002) and Podlubny (2005).
4 Subordination, solution
The principle of subordination is a standard procedure for solving problems of the type considered in this work (Bazhlekova and Bazhlekov, 2019; Benson et al., 2000; Luchko, 2019; Raghavan and Chen, 2015; Stanislavsky et al., 2014). The procedure recasts the classical Green’s (or source) function in the classical diffusion equation, often called the parent process as in Benson et al. (2000), by expressing it in terms of an operational time, τ; that is, the classical solution is now expressed as ( rather than . Operational time, on the other hand, is defined through a subordinator, , a nondecreasing and right continuous function which ensures that its inverse subordinator, , exists. In essence, as outlined in Raghavan and Chen (2015) if the pressure distribution reflecting the process under consideration were , then and are related by the integral transform
where is the probability density of the inverse subordinator, ; it contains all the requirements for the process under consideration (subdiffusion, superdiffusion or a combination). For example, Raghavan and Chen (2015) show that for subdiffusive processes where is
Fortunately, for the combined case, the needed expressions, and are available (Bazhlekova, 2019; Bazhlekova and Bazhlekov, 2019; Luchko, 2019). Analogous to equation (16) for the combined process, we may show that the Green’s function, (x, t) is
where (Bazhlekova, 2019; Bazhlekova and Bazhlekov, 2019)
with being the classical expression given by
As noted in Raghavan and Ozkan (1994) and Raghavan and Chen (2018), Green’s (or source) functions in transient diffusion in porous media reflect the instantaneous change in pressure, d(Δp)/dt, and on using the results of Bazhlekova and Bazhlekov (2019) given by
we arrive at the expression for the logarithmicpressurederivative or pressure derivative, p′,
The symbols , t_{D}, and x_{D}, reflect the dimensionless pressure, dimensionless time and dimensionless distance, respectively, defined by the following expressions
and
In the following, as all measurements are expressed with respect to the distance from a linesource well, henceforth, we replace x_{D} by r_{D}, the distance from the center of the well; we take l_{ref} to be the well radius, r_{w}.
The analog of equation (22) for the homogeneous case α = 1 and β = 1 is
To our knowledge, the logarithmic pressure derivative was introduced by Chow (1952). But it was not adopted presumably because of the lack of precision in pressure measuring devices that measured pressures suitable for computing pressure derivatives. Procedures such as the maximum slope method of GarciaRivera and Raghavan (1979) are an attempt to circumvent this limitation. The introduction of high precision gauges in the late 1970s was a prime factor in the widespread adoption of the pressure derivative technique; see Bourdet et al. (1983). For completeness, we note that other methods to arrive at equation (22) are available in Luchko (2016, 2017) as well as Boyadjiev and Luchko (2017).
For the special situation 2α = β + 1, we have the closed form expression (Luchko, 2016)
where
The righthand side of equation (28) is the grouping we would expect to use based on our knowledge of the Theis (1935) solution and the subdiffusive solutions in Raghavan and Chen (2019); in fact, Gehlhausen (2009) presumes this to be the case a priori in all situations which do not appear to be the case, as we shall see. This expression is important in many ways. First, if α and β were both 1 it reduces to the classsical expression, the expression that follows from the Exponential Integral solution; see equation (26). Second, it suggests that ideas which we intuitively expect from classical diffusion and even subdiffusion should hold in the combined situation for this condition. Third, it provides a method for correlating and understanding transient responses in a general way; basically one may show all results on a single chart if α were used as a parameter of interest much like the Theis (1935) solution.
Most important, however, it is not clear whether solutions may be correlated in terms of in general. This matter would have to be determined solely through computations, although we do get a glimpse of what is to be expected through considering the work of Deng and Schilling (2019) who have presented the longtime approximation corresponding to the above expression in equation (22) given by
where n is number of dimensions. Again, interestingly, many observations are evident: this expression suggests that the time dependence of the above expression is essentially similar to that of Raghavan and Chen (2019) for purely subdiffusive flow where (with the influence of the logarithmic term being small) as is given by
and where ψ (·) is the Digamma function. Most interesting and surprising, the superdiffusion component plays no role in the exponent of time, t, which is not the case in 1D flow as discussed in Chen and Raghavan (2015); furthermore, the expression makes clear that in the general case unless 2α were equal to β + 1. This issue could become important when we consider more complex geometries for wellbores such as fractured or horizontal wells which are not considered here. Other points of note in this regard are explored below. Note that for α = 1, equation (30) does yield the semilogarithmic approximation of the Theis (1935) solution, namely the Cooper and Jacob (1946) approximation (also, see Eq. (32)).
In concluding this section, we note that the results in EstradaRodriguez et al. (2018) may also be useful and they have derived Green’s functions.
5 Computational results
The principal goal of the computations considered below is to establish the viability of the proposals made to ensure that they are suitable for applications. Second, we establish the validity of the observations made concerning the structure of the solutions; these are particularly revealing in that they illustrate the nature of downhole transients to arrive at reservoir descriptions. Third, we demonstrate some unusual features of superdiffusion which are surprising at least to us. All solutions discussed here were computed through equation (22). Asymptotic behaviors are demonstrated through equation (29) or through other expressions derived in this work. For convenience and also consistency, computations are displayed in terms of . Pressure responses are computed by straight forward integration of
Figures 1 and 2 present results wherein 2α = β + 1; as indicated by equation (22) we expect solutions to be correlatable in terms of which is the case. Figure 1 considers responses in terms of the logarithmic derivative, , for four values of r_{D} [1 (unbroken line), 25 (circles), 50 (squares) and 100 (triangles)]. The results in Figure 1 are for α = 0.9 and β = 0.8. In actuality, we have compared and correlated solutions for r_{D} as large as 250. The dashed line reflects the expression on the righthand side of equation (29); at long enough times . For completeness the derivative response for the Theis (1935) solution is also shown as a dashed line; see equation (26). Not unexpectedly, responses for the classical case are quite distinct from the heterogeneous solutions. Based on this result, we expect the dimensionless pressure, p_{D}(r_{D},t_{D}), to also correlate as a function of ; see Figure 2. Here, we present both p_{D}(r_{D},t_{D}) and and two combinations of α and β (α = 0.9, β = 0.8 and α = 0.8, β = 0.6). The correlation in terms of the circles, triangles, etc. is compatible with the results shown in Figure 1 for α = 0.9, β = 0.8. The bottom two curves (shown as dashes) are the corresponding responses; at long enough times they, as already noted, will follow the trend predicted by equation (29); the circles denote the start of the straight line, equation (29). Pressure responses at long enough times, , do also follow the trend provided that ; the longtime expression derived through the Laplace transformation corresponds to:
and, incidentally, for purely subdiffusive flows, β = 1, we have
Fig. 1 Derivative responses reflecting anomalous diffusion (super and subdiffusion) under the constraint 2α = β + 1, equation (22). Three values of r_{D} besides the wellbore shown as an unbroken line are considered to demonstrate that pressure responses are a function of . The dashed line reflects equation (29). For each value of r_{D}, in agreement with equation (29). 
Fig. 2 Pressure and derivative distributions reflecting anomalous diffusion (super and subdiffusion) under the constraint 2α = β + 1, equation (22). Results shown apply for all values of r_{D}. Two sets of α and β values are considered; for one of the cases we specifically show that results may be correlated in terms of . For each value of r_{D}, both p_{D} (r_{D}, t_{D} → ∞) and . The topmost dashed line (α = 0.8, β = 0.6) reflects the trend noted by equation (31). 
The straight line commensurate to equation (31), however, unlike the derivative responses which begin much earlier, is only evident at times much longer than those considered in Figure 2; and this is often the case. The topmost dashed line shows the pressure trend reflected by equation (31) for α = 0.8, β = 0.6; it is clear that the line will merge with the pressure response.
Figure 3 considers responses where pure superdiffusion is the governing mechanism; that is α = 1. Responses considered are at the wellbore and are quite distinct from those considered so far or those shown in the literature for subdiffusive flow. The curve β = 1 is the linesource solution and as expected. For the other cases although equation (31) does not strictly apply if α = 1, by considering α = 1 − ε we expect p_{D} (r_{D}, t_{D} → ∞) ≈ a constant for β < 1; that is, we expect responses to be akin to a constant pressure boundary. This is the case as both curves in Figure 3 display a downward trend. As β decreases further the trend to a steadystatetype response will be even more rapid. Recall that all this is occurring in a system that is infinite in its areal extent and where steady flow should not be the norm. In summary, for the situation under discussion, pure superdiffusion in 2D leads to, as alluded earlier a steady type behavior at long enough times. It does not appear to be possible to make reference to any specific characteristic that would by itself permit us to readily identify superdiffusion in the pressure trace which would lead to the determination of β in a manner similar through equation (29). Note that in equation (32) the logarithmic nature of the pressure response is preserved for α = 1 whereas in the case of equation (31) no conclusion may be drawn. Presumably because of the lack of an expression similar to equation (31) neither Cloot and Botha (2006), nor Gehlhausen (2009) recognize the possibility of steady flow. The above observations also apply to r_{D} > 1. The slopes of the downward sections of the derivative curves are independent of r_{D} and depend only on β. As no analytical options are available at this time, one option is to determine β empirically through the exponents after determining the slope of the downward section; see Figure 4 which may be represented by
Fig. 3 Pressure and derivative responses at the wellbore, r_{D} = 1, reflecting the influence of pure anomalous, superdiffusion; that is, responses obtained under the presumption that α = 1. The downward trend in the derivative curves is to be expected as the pressure response at the well approaches a steady response as the predicted by equations (22) and (29). 
Fig. 4 Exponents of the downward sections of the derivative curves for superdiffusion to determine β. 
Simply put, in essence, pure superdiffusion leads to negative slopes at long times.
Because of the unusual aspects seen in Figure 3, Figure 5 explores the characteristic behaviors near the classicalsubdiffusive limit by considering the role of subdiffusion in the range and we consider for four values of α at the wellbore in Figure 5 (0:9000; 0:9250; 0:9900 and 0:9999). Pressure responses are shown for only α of 0.9000. Straight lines corresponding to the asymptotic solution, namely equation (31), are evident for . For larger values no straight line is evident for the time ranges considered here.
Fig. 5 Pressure and derivative responses at the wellbore, r_{D} = 1, reflecting influence of subdifussion near the classical limit of α of 1 under the presumption that both super and subdiffusion effects govern transients. Downward trend in the derivative curves for α close to the limit of classical diffusion suggests that the well pressure approaches a steady value. 
Returning to the results in Figure 3, we call attention to the observation not explicitly detailed earlier that in linear systems the pressure does depend on the longrange exponent, β, for α = 1, for Chen and Raghavan (2015) show that
Furthermore, one important insight that should be garnered on comparing equation (34) with equation (27) is that for 1D flow the response for the neutral case, 2α = β + 1, equation (34) yields the homogeneous solution; see Magin et al. (2013) and Chen and Raghavan (2015) and this is not the case in 2D. Incidentally, the works by Chu and coworkers tellingly illustrate advantages that accrue through the use of equation (34).
Figure 6 considers the general case of combined diffusion for the situation ; nevertheless pressure distributions are plotted in terms of for ease of comparison. Responses are typical of the myriad solutions considered by us. It is clear that in general solutions depend on both α and β; the general shapes are no different from those considered earlier. Indeed at long enough times the exponent of t is independent of β, and equation (29) does remain available. The lines are not shown here as the onset of the straight lines happens to be much later than the times considered here. Longtime responses, however, are correlatable in terms of ; see equation (27).
Fig. 6 Pressure derivative distributions reflecting anomalous diffusion (super and subdiffusion) for the general case. The constraint, 2α = β + 1, does not apply. Results are presented in terms of principally for convenience. 
5 Discussion and concluding remarks
The principal goal of this work is to incorporate the influences of macroand microscale heterogeneities at the Theis scale around a linesource well. Anomalous processes of fractional calculus are used to evaluate pressure responses governing spacetime fractional diffusion. Two dimensional super and subdiffusion solutions are, to our knowledge, yet to be explored. One immediate advantage of this work is that it provides an understanding of long term behaviors of the flow of fluids produced through complex wellbores (fractured wells, horizontal wells), for pseudoradial flow will ultimately prevail in all situations. For example, if a well produces through a finiteconductivity fracture, then exponents of and 1α corresponding to the bilinear, linear and pseudoradial flowperiods, respectively, should be evident on the derivative trace whenever all three flow regimes are evident; see Raghavan and Chen (2020). Responses to a linesource well are presented both analytically and graphically. Salient features of the results obtained are as follows: (i) characteristic features of super and subdiffusion and pure subdiffusion at long enough times are identical, (ii) derivative curves for pure superdiffusion yield negative slopes, and (iii) pressure distributions fall under two classifications, those that meet the restriction 2α = β + 1 in which case , the grouping expected from the Theis solution, and the general case wherein ; regardless, at long enough times provided that . Computational results presented here may be used to address many applications in the earth sciences: movement of groundwater, production of hydrocarbon and geothermal resources, sequestration etc. Studies of this type provide a richer understanding of transient diffusion in heterogeneous porous solids and may be extended to other similar applications.
Acknowledgments
We thank Dr. Wei Chun Chu for reviewing the typescript and his comments.
References
 Ali I., Malik N., Chanane B. (2014) Fractional diffusion model for transport through porous media, in 5th International Conference on Porous Media and Their Applications in Science, Engineering and Industry, K. Vafai, A. Bejan, A. Nakayama, O. Manca (eds). ECI Symposium Series. [Google Scholar]
 Angulo J.M., RuizMedina M.D., Anh V.V., Grecksch W. (2000) Fractional diffusion and fractional heat equation, Adv. Appl. Prob. 32, 4, 1077–1099. [Google Scholar]
 Balakrishnan A.V. (1960) Fractional powers of closed operators and the semigroups generated by them, Pacific J. Math. 10, 2, 419–437. https://doi.org/10.2140/pjm.1960.10.419. [Google Scholar]
 Bazhlekova E. (2019) Subordination principle for spacetime fractional evolution equations and some applications, Integr. Transform. Spec. Funct., 431–452 https://doi.org/10.1080/10652469.2019.1581186. [Google Scholar]
 Bazhlekova E., Bazhlekov I. (2019) Subordination approach to spacetime fractional diffusion, Mathematics 2019, 7, 415. https://doi.org/10.3390/math7050415. [Google Scholar]
 Beier R.A. (1990) Pressure Transient Field Data Showing Fractal Reservoir Structure, in: Paper 90–04 presented at the Annual Technical Meeting, Calgary, Alberta Petroleum Society of Canada. https://doi.org/10.2118/9004. [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]
 Bernard S., Delay F., Porel G. (2006) A new method of data inversion for the identification of fractal characteristics and homogenization scale from hydraulic pumping tests in fractured aquifers, J. Hydrol. 328, 3–4, 647–658. https://doi.org/10.1016/j.jhydrol.2006.01.008. [Google Scholar]
 Bourdet D., Whittle T.M., Douglas A.A., Pirad Y.M. (1983) A new set of type curves simplifies well test analysis, World Oil 95–106. [Google Scholar]
 Boyadjiev L., Luchko Y. (2017) Mellin integral transform approach to analyze the multidimensional diffusionwave equations, Chaos Solit. Fract. 102, 127–134. https://doi.org/10.1016/j.chaos.2017.03.050. [Google Scholar]
 Cai M., Li C. (2019) On Riesz derivative, Fract. Calc. Appl. Anal. 22, 2, 287–301. https://doi.org/10.1515/fca20190019. [Google Scholar]
 Caputo M. (1967) Linear models of dissipation whose Q is almost Frequency IndependentII, Geophys. J. R. Astron. Soc. 13, 5, 529–539. [Google Scholar]
 Chang J., Yortsos Y.C. (1990) Pressuretransient analysis of fractal reservoirs, SPE Form. Eval. 5, 1, 31–39. [Google Scholar]
 Chang A., Sun H., Zhang Y., Zheng C., Min F. (2019) Spatial fractional Darcy’s law to quantify fluid flow in natural reservoirs, Phys. A Stat. Mech. Its Appl., 519, 119–126. [Google Scholar]
 Chen C., Raghavan R. (2015) Transient flow in a linear reservoir for spacetime fractional diffusion, J. Pet. Sci. Eng. 128, 194–202. https://doi.org/10.1016/j.petrol.2015.02.021. [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. [Google Scholar]
 Chu W., Pandya N., Flumerfelt R.W., Chen C. (2019a) 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. (2019b) A new technique for quantifying pressure interference in fractured horizontal shale wells, SPE Res. Eval. Eng. 23, 01, 143–157. https://doi.org/10.2118/191407PA. [Google Scholar]
 Cloot A., Botha J.F. (2006) A generalised groundwater flow equation using the concept of noninteger order derivatives, Water SA 32, 1, 1–7. www.wrc.org.za. [Google Scholar]
 Cooper H.H., Jacob C.E. (1946) A generalized graphical method for evaluating formation constants and summarizing wellfield history, Trans. AGU 27, 526–534. [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]
 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]
 Deng C.S., Schilling R.L. (2019) Exact asymptotic formulas for the heat kernels of space and timefractional equations, Fract. Calc. Appl. Anal. 22, 4, 968–989. https://doi.org/10.1515/fca20190052. [Google Scholar]
 Doe T., Shi C., Knitter C., Enachescu C. (2014) Discrete fracture network simulation of production data from unconventional wells, paper URTeC 1923802, in: Proceedings of The Unconventional Resources Technology Conference, Denver CO. [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]
 EstradaRodriguez G., Gimperlein H., Painter K.J., Stocek J. (2018) Spacetime fractional diffusion in cell movement models with delay, Math. Models Methods Appl. Sci. 29, 01, 1–2. https://doi.org/10.1142/S0218202519500039. [Google Scholar]
 FlamencoLopez F., CamachoVelazquez R. (2001) Fractal transient pressure behavior of naturally fractured reservoirs, in: Paper 71591 presented at the Annual Technical Conference and Exhibition, Soc. Pet. Eng., New Orleans, LA. https://doi.org/10.2118/71591MS. [Google Scholar]
 Fomin S., Chugunov V., Hashida T. (2011) Mathematical modeling of Anomalous Diffusion in Porous Media, Fractional Differential Calculus 1, 1–28. [Google Scholar]
 GarciaRivera J., Raghavan R. (1979) Analysis of short time pressure transient data dominated by wellbore storage and skin, J. Pet. Tech. 31, 5, 623–631. [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. [Google Scholar]
 Gehlhausen A.N. (2009) Evaluation of the fractional Theis solution, Master of Science Thesis, University of Nevada, Reno. http://hdl.handle.net/11714/4109. [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]
 Gorenflo R., Mainardi F. (1998) Random walk models for spacefractional diffusion processes, Fract. Calc. Appl. Anal. 1, 2, 167–191. http://www.math.bas.bg/fcaa/. [Google Scholar]
 Gorenflo R., Kilbas A.A., Mainardi F., Rogosin S.V. (2014) MittagLeffler Functions, Related Topics and Applications, Springer, Berlin Heidelberg, p. 443. [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. [Google Scholar]
 Henry B.I., Langlands T.A.M., Straka P. (2010) An introduction to fractional diffusion, in: Presented at the Conference: Complex Physical, Biophysical and Econophysical Systems – Proceedings of the 22nd Canberra International Physics Summer School, pp. 37–89. [Google Scholar]
 Huang F., Liu F. (2005) The fundamental solution of the spacetime fractional advectiondispersion equation, J. Appl. Math. Comput. 18, 1–2, 339–350. https://doi.org/10.1007/BF02936577. [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]
 Kwaśnicki M. (2017) Ten equivalent definitions of the fractional Laplace operator, Fract Calc Appl Anal. 20, 1, 7–51. [Google Scholar]
 Le Mehaute A., Crepy G. (1983) Introduction to transfer and motion in fractal media: The geometry of kinetics, Solid State Ionics 1, 9–10, 17–30. [Google Scholar]
 Lenormand R. (1992) On use of fractional derivatives for fluid flow in heterogeneous media, in: Proceedings 3rd European Conference on the Mathematics of Oil Recovery, Delft The Netherlands. [Google Scholar]
 Luchko Y. (2016) A new fractional calculus model for the twodimensional anomalous diffusion and its analysis, Math. Model. Nat. Phenom. 11, 3, 1–17. [Google Scholar]
 Luchko Y. (2017) On some new properties of the fundamental solution to the multidimensional space and timefractional diffusionwave equation, Mathematics 5, 4, 76. https://doi.org/10.3390/math5040076. [Google Scholar]
 Luchko Y. (2019) Subordination principles for the multidimensional spacetimefractional diffusionwave equation, Theory Probab. Math. Stat. 98, 121–141. [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, Microporous Mesoporous Mater. 178, 15, 39–43. https://doi.org/10.1016/j.micromeso.2013.02.054. [Google Scholar]
 Mainardi F., Luchko Y., Pagnini G. (2001) The fundamental solution of the spacetime fractional diffusion equation, Fract. Calc. Appl. Anal. 4, 2, 153–192. [Google Scholar]
 Mainardi F. (2010) Fractional calculus and waves in linear viscoelasticity, Imperial College Press, London, p. 344. [Google Scholar]
 Mattax C.C., Dalton R.L. (eds) (1990) Reservoir simulation, SPE Monograph Series 13, 187 pp. [Google Scholar]
 Meerschaert M.M., Zhang Y., Baeumer B. (2008) Tempered anomalous diffusion in heterogeneous systems, Geophys. Res. Lett. 35, L17403. https://doi.org/10.1029/2008GL034899. [Google Scholar]
 Metzler R., Glockle W.G., Nonnenmacher T.F. (1994) Fractional model equation for anomalous diffusion, Physica A 211, 1, 13–24. [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]
 Moodie T.B., Tait R. (1983) On thermal transients with finite wave speeds, J. Acta Mechanica 50, 1–2, 97–104. https://doi.org/10.1007/BF01170443. [Google Scholar]
 Nigmatullin R. (1984) To the theoretical explanation of the universal response, Phys. Status Solidi B Basic Res. 123, 2, 739–745. [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. [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. [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. [Google Scholar]
 Oldham K.B., Spanier J. (1974) The fractional calculus; theory and applications of differentiation and integration to arbitrary order, Academic Press, New York, p. 234. [Google Scholar]
 Park H.W., Choe J., Kang J.M. (2000) Pressure behavior of transport in fractal porous media using a fractional calculus approach, Energy Sources, Part A: Recovery, Utilization, and Environmental Effects 22, 10, 881–890. [Google Scholar]
 Podlubny I. (1998) Fractional differential equations. An introduction to fractional derivatives, fractional differential equations, some methods of their solution and some of their applications, Academic Press, New York, p. 340. [Google Scholar]
 Podlubny I. (2005) MittagLeffler function. www.mathworks.com/matlabcentral/fileexchange/8738mittaglefflerfunction. [Google Scholar]
 Raghavan R., Ozkan E. (1994) A method for computing unsteady flows in porous media, Pitman Research Notes in Mathematics Series (318), Longman Scientific & Technical, Harlow, UK, p. 188. [Google Scholar]
 Raghavan R., Chen C. (2015) Anomalous subdiffusion to a horizontal well by a subordinator, Transp. Porous Med. 107, 387–401. https://doi.org/10.1007/s112420140444y. [Google Scholar]
 Raghavan R., Chen C. (2018) 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) The Theis solution for subdiffusive flow in rocks, Oil Gas Sci. Technol.  Rev. IFP Energies nouvelles 74, 6. https://doi.org/10.2516/ogst/2018081. [Google Scholar]
 Raghavan R., Chen C. (2020) A study in fractional diffusion: Fractured rocks produced through horizontal wells with multiple, hydraulic fractures, Oil Gas Sci. Technol.  Rev. IFP Energies nouvelles 75, 68. https://doi.org/10.2516/ogst/2020062. [Google Scholar]
 Sapora A., Cornetti P., Chiaia B., Lenzi E.K. (2017) Nonlocal diffusion in porous media: a spatial fractional approach, J. Eng. Mech. 143, 5, D40160071–D40160077. [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, paper URTEC2154675, in: Proceedings of Unconventional Resources Technology Conference, San Antonio, Texas. https://doi.org/10.15530/URTEC20152154675. [Google Scholar]
 Sinkov K., Weng X., Kresse O. (2021) Modeling of proppant distribution during fracturing of multiple perforation clusters in horizontal wells, paper SPE204207MS, in: Presented at the SPE Hydraulic Fracturing Technology Conference and Exhibition. https://doi.org/10.2118/204207MS. [Google Scholar]
 Stanislavsky A., Weron K., Weron A. (2014) Anomalous diffusion with transient subordinators: A link to compound relaxation laws, J. Chem. Phys. 140, 5, 054113. https://doi.org/10.1063/1.4863995. [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. [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. [Google Scholar]
 Suzuki A.T., Hashida K.Li, 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. Natural Gas Geosci. 1, 6, 445–455. https://doi.org/10.1016/j.jnggs.2016.11.009. [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. Trans. AGU 16, 2, 519–524. [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 Reserv. Eval. Eng. 8, 3, 248–254. https://doi.org/10.2118/77452PA. [Google Scholar]
 Uchaikin V.V. (2012) Fractional derivatives for physicists and engineers; background and theory, Vol. 1, SpringerVerlag, New York, NY, p. 385. [Google Scholar]
 Wang Y. (2013) Anomalous transport in weakly heterogeneous geological porous media, Phys. Rev. E 87, 03214410321448. https://doi.org/10.1103/PhysRevE.87.032144. [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. [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. [Google Scholar]
 Zhokh A., Strizhak P. (2019) Investigation of the anomalous diffusion in the porous media: a spatiotemporal scaling. Heat Mass Transf. 55, 2693–2702. https://link.springer.com/article/10.1007/s00231019026024. [Google Scholar]
All Figures
Fig. 1 Derivative responses reflecting anomalous diffusion (super and subdiffusion) under the constraint 2α = β + 1, equation (22). Three values of r_{D} besides the wellbore shown as an unbroken line are considered to demonstrate that pressure responses are a function of . The dashed line reflects equation (29). For each value of r_{D}, in agreement with equation (29). 

In the text 
Fig. 2 Pressure and derivative distributions reflecting anomalous diffusion (super and subdiffusion) under the constraint 2α = β + 1, equation (22). Results shown apply for all values of r_{D}. Two sets of α and β values are considered; for one of the cases we specifically show that results may be correlated in terms of . For each value of r_{D}, both p_{D} (r_{D}, t_{D} → ∞) and . The topmost dashed line (α = 0.8, β = 0.6) reflects the trend noted by equation (31). 

In the text 
Fig. 3 Pressure and derivative responses at the wellbore, r_{D} = 1, reflecting the influence of pure anomalous, superdiffusion; that is, responses obtained under the presumption that α = 1. The downward trend in the derivative curves is to be expected as the pressure response at the well approaches a steady response as the predicted by equations (22) and (29). 

In the text 
Fig. 4 Exponents of the downward sections of the derivative curves for superdiffusion to determine β. 

In the text 
Fig. 5 Pressure and derivative responses at the wellbore, r_{D} = 1, reflecting influence of subdifussion near the classical limit of α of 1 under the presumption that both super and subdiffusion effects govern transients. Downward trend in the derivative curves for α close to the limit of classical diffusion suggests that the well pressure approaches a steady value. 

In the text 
Fig. 6 Pressure derivative distributions reflecting anomalous diffusion (super and subdiffusion) for the general case. The constraint, 2α = β + 1, does not apply. Results are presented in terms of principally for convenience. 

In the text 