A study in fractional diffusion: Fractured rocks produced through horizontal wells with multiple, hydraulic fractures

The spatiotemporal evolution of transients in fractured rocks often displays unusual characteristics and is traced to multifaceted origins such as micro-discontinuity in rock properties, rock fragmentation, long-range connectivity and complex flow paths. A physics-based 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 long-range 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.

Abstract. The spatiotemporal evolution of transients in fractured rocks often displays unusual characteristics and is traced to multifaceted origins such as micro-discontinuity in rock properties, rock fragmentation, long-range connectivity and complex flow paths. A physics-based 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 long-range 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.

Introduction
The vast array of multiple mathematical models proposed to evaluate fractured reservoirs attest to their complexity. It is generally agreed that fluid is transported along preferred pathways with the pathways supported by a medium of lower conductivity. Both the pathways and lower conductivity system are complex geologically and in terms of geological gradients. Some of the complexity is usually better understood subsequent to the drilling of new (additional) wells; see Raghavan (2004). This view is a result of both mapping and a well-defined, spatiotemporal evolution of responses at wells. Broadly speaking, the types of models used today may be classified in two ways: (i) those that are known as equivalent porous media models (Yeh et al., 2015) as discussed in Pruess and Narasimhan (1985), Haggerty and Gorelick (1995), and Noetinger et al. (2001) which are extensions of the work of Barenblatt et al. (1960) and Warren and Root (1963) and, (ii) those that are known as Distributed Fracture Network (DFN) schemes where fractures are introduced as discontinuities with specific properties (orientation, aperture, size, position, etc.) in a uniform continuum and explicitly model pathways representing natural fractures (Artus, 2020;Cacas et al., 1990a, b;Dershowitz et al., 2002;Karimi-Fard et al., 2004;Noetinger, 2015). The equivalence of the two approaches is compared in Dong et al. (2019) and the circumstances in which the DFN network may be treated as a multiple, interacting continuum is noted there. Irrespective of the scheme adopted, one may not be indifferent to the manifold aspects of the transport of fluids and treat the system as an undifferentiated continuum concerning properties. Raghavan et al. (2001) elucidate the complexities and intrinsic connections of addressing the origins of reservoir heterogeneity and demonstrate the prodigious amount of work involved in dataset curation. Here as noted in Cacas et al. (2001) we must integrate of matters relating to fracture permeability, fracture aperture and fracture connectivity to mention only a few aspects; see Belayneh et al. (2006), Kang et al. (2015) and Bisdom et al. (2016). One may attempt to address these manifold features through either of the models noted above. Here our focus is on the use of equivalent-porous-media concepts based on anomalous diffusion by fractional calculus (Albinali et al., 2016;Fomin et al., 2011;Holy and Ozkan, 2016;Raghavan, 2011). In this scheme, rather than following the classical time-distance relationship, r 2~t , with r being the distance and t being time; the transient follows the relation r 2~ta where a is a constant with a < 1 (subdiffusion) or with a > 1 (superdiffusion), or even with a = 1 (combined influences). Chu et al. (2019Chu et al. ( , 2020 discuss the advantages of addressing atypical behaviors employing a subdiffusive model to evaluate wells of the Permian basin produced through multiple hydraulic fractures (for both single and multiple well-tests). Quantitative estimates of the exponent a for the tests discussed in Chu et al. (2019) are given in Raghavan and Chen (2018b). Our experience, confirmed by Chu (2018), and illustrated with an example here indicates that one must also consider superdiffusive flows in the fracture system. Significant challenges exist in addressing superdiffusive problems analytically; in fact, this is a fraught exercise. For example, Defterli et al. (2015) observe that this is an open problem and Metzler et al. (2007) note that the method of images may not be used to address reflection conditions; that is, closed boundary effects may not be modelled by image wells. One way to circumvent this issue is by considering approximate or asymptotic cases as we demonstrate here. Similarly, obtaining solutions of interest for 2D problems is difficult. On the basis of O'Shaughnessy and Procaccia (1985a, b), many studies have proposed solutions for subdiffusive, transient flows based on rules for fractal objects (Chang and Yortsos, 1990). For vertical wells, such solutions lead to an expression where the pressure drop, Dp(r, t), may be expressed as a complementary incomplete gamma function: Cða; xÞ ¼ R 1 x f aÀ1 expðÀfÞdf with the principal conclusion suggesting that at long enough times transients display characteristic power-law features, namely, Dp~t # , where # is a constant, contrary to the widely known result Dp~lnt. The finding that Dp~C(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 power-law features when conduits are not characterized as fractals; a similar observation for fluvial deposits may be found in Thomas et al. (2005). Fractional laws proposed in this paper for both sub-and super-diffusion preclude conjuring counterfactual geological options.
Our objective is to address patterns of derivative responses shown in Figure 1, processed and presented by Doe et al. (2014). Early-time responses, as is usual in most wells producing unconventional reservoirs, reflect the lack of appropriate measurements particularly bottom-hole pressures; see Scott et al. (2015). The principal point to note here is that all responses at late times display a trend well above what would expect of models based on classical diffusion. None of the responses display slopes of 1/2, 1/4, or 1/8 and it is not appropriate to force a straight line fit with such slopes. Our point is that if we are trying to address heterogeneous systems, we are required to go beyond classical systems. We offer one choice to do so. Regardless of the approach taken, three elements, as incorporated here, are paramount: a matrix system, a fracture or fissure system and a hydraulic fracture.
As in the classical case, we develop the model by combining the conservation of mass principle with the equation of state for a slightly compressible fluid and a fractional flux law. It is needless to emphasize the importance of a flux law as it is the marrow of reservoir-engineering predictions. Even in the classical case, the deterministic situation, a model derived by the three principles just noted, the basic driving force is attributed to random motions namely the random walk (Doyle and Snell, 1984;Klafter and Sokolov, 2011) that consists of independent and identically distributed steps, X n ; the first moment (the conservation of mass principle) and second moment (the time-distance relationship) being finite. Inherent in the classical Darcy law is the concept that the porous rock is a uniform body as the overall size of the porous body is much, much larger than the sizes the individual pores that constitute the porous rock and that any potential surface is a notional surface in that a particular cross-sectional element is not unlike any other and the normal gradient reflects the overall flux. If, however, we postulate that this is no longer the case, the system is heterogeneous and distant points connected, then modification is required; one such option being to consider anomalous diffusion with random jumps having a power law probability tail in which case we use an equation of the form o t p ¼ o bþ1 x p for 0 < b < 1 1 with b reflecting the power law tail (PðjX n j > rÞ % r Àðbþ1Þ ); see Benson et al. (2000Benson et al. ( , 2004. Levy-stable-particles trace a path of dimension b + 1 namely the exponent of the x-derivative, o x , in o t p ¼ o bþ1 x p; in other words, the Hurst index, H, is 1/(b + 1) or hx 2 i~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 o a t p ¼ o 2 x p for 0 < a < 1 rather than o t p ¼ o bþ1 x p; in this case PðT n > tÞ % t Àa for 0 < a < 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 waiting-time solutions as documented in Henry et al. (2010). In fact, the expression o a t p ¼ o 2 x p may be arrived at in a number of ways; see Montroll and Weiss (1965), Kenkre et al. (1973), Hilfer and Anton (1995). The exponent a reflects the anomalous diffusion coefficient (random walk dimension), d w (Metzler et al., 1994) for diffusion on fractal objects. The earth-science literature abounds with examples of situations where subdiffusive models apply, namely in situations wherein cracks, fissures and cervices abound (Caine et al., 1996;Evans, 1988;Mitchell and Faulkner, 2009;Savage and Brodsky, 2011;Scholz et al., 1993). Mathematical models depicting flows in porous rocks may be found in Jourde et al. (2002), Cortis and Knudby (2006) and Suzuki et al. (2016). In combined situations, we consider o a t p ¼ o bþ1 x p for 0 < a, b < 1 while recognizing that for combinations of a and b 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 x p for 0 < b < 2. We simply note, however, that the range of b may have to be circumscribed to solve for certain auxiliary conditions. We shall return to this point later.
2 For a discussion of addressing superdiffusion through the exponent a see Metzler et al. (2014). We do not consider this perspective. It is an open problem for situations considered here. nonlocality by deriving a tractable, analytical solution. As noted in Ozkan and Raghavan (1991), the first step involves the development of point-source solutions that satisfy the specifications of the desired model. Should such a solution be available, as it is in many cases, prediction of the performance of a fractured well producing through a horizontal well through multiple fractures is tractable. Here, of course, we consider situations where fluids are transported in the fracture and matrix systems under anomalous diffusion, specifically, through fluid movement in the fracture system that permits for conditions instrumental in both super-and subdiffusive transport and rock properties in the matrix system conducive to subdiffusive flow.
Thus transient diffusion considered here is not construed to be diffusion in homogeneous systems in the manner envisaged by early models of naturally fractured reservoirs discussed in Barenblatt et al. (1960), Warren and Root (1963), Kazemi (1969) and de Swaan-O (1976) although we do consider responses at the identical scale. Probability density functions that result from solving fractional differential equations are different from those generated by the equation governing classical diffusion and are examined here.

The model
The mathematical model consists of three regions: the fracture system, X f , the matrix system, X m , and hydraulic fracture, X F . Figure 2 taken from Raghavan and Ohaeri (1981) is a portrayal of the model put forward in de Swaan-O (1976) which shows the colligation of the matrix and fracture elements where 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 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, X m , occurs in the usual way as it is in perfect contact with the fracture system X f (Cinco-Ley and Meng, 1988;Houze et al., 1988;Yeh et al., 1986). All fluid produced through the hydraulic fracture may enter the hydraulic fracture only through the fracture system; that is, there is no communication between the matrix system and the hydraulic fracture. Perfect contact is assumed at each interface where there is hydraulic communication. The fracture tip is assumed to be sealed. As noted earlier, our representation of the fractured rock incorporates both super-and subdiffusive effects to model flow in the fracture system while flow in the matrix is presumed to be subdiffusive. Darcy flow is assumed in the hydraulic fracture. To consider a system of multiple fractures intercepting a horizontal wellbore and producing against a common wellbore pressure, we follow the outlines of Chen and Raghavan (2013) through Duhamel's theorem.
The fractional flux law discussed in Fomin et al. (2011), Magin et al. (2013), and Chen and Raghavan (2015) uses the fractional derivative, o a /ot a f(t), in the form defined by Caputo (1967) o a ot a f ðtÞ ¼  Kazemi (1969) and de Swaan-O (1976); after Raghavan and Ohaeri (1981). and assumes to govern transients in the fractured rock by an expression of the form for the velocity, v(x, t): with the exponents a and b both reflecting the heterogeneities in the fractured rock. 3 Note that unlike most studies a one-sided form of the spatial fractional derivative is used. Denoting the exponents of the fracture system by the subscript f and that of the matrix system by the subscript m, a f and b f are presumed to be 1 as we consider both sub-and superdiffusive influences in the fracture system. For the matrix system, a m is presumed to be 1 while b m is taken to be 1 because only subdiffusive influences are considered.

The partial differential equations, associated conditions
The three partial differential equations to be solved to obtain the pressure distributions in the fracture system (Region X f ), matrix system (Region X m ) and the hydraulic fracture (Region X F ) for t > 0, are respectively: and Here p i the initial pressure is identical in all parts of the system at t = 0, and Dp 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 no-flow boundary; thus Further, for perfect contact we require as the requirement of continuity in pressure at the fracturematrix interface, X f À X 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 X f À X 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 and The conditions to be satisfied at the fracture-reservoir interface, X f -X F , are, again, the continuity of pressure and the continuity of flux; in consequence, respectively, where the symbolq represents the flux.
As is usual we follow the outline of the idealization given in Cinco- Ley and Meng (1988) to consider flow in the hydraulic fracture; that being so, the following assumptions apply: (i) the tips of the hydraulic fracture are sealed, (ii) the width of the fracture is negligibly small compared with its length; thus, one may use the pseudofunction approximation for flow the direction normal to its length, and (iii) the volume in the fracture is inconsequential compared with the reservoir volume and therefore we presume flow is steady in the fracture. The assumptions considered hitherto lead to the consideration that the fracture is a source term but of finite conductivity, 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.

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 3 The fractional flux law for space fractional diffusion may be stated in other forms. Often a pde in the form of o t p ¼ o b x p with 0 < b < 2 is considered to model superdiffusion. In such cases a fractional flux law of the form vðbÞ ¼ ðk 1;b =lÞo bÀ1 x p needs to be used; see Benson et al. (2004). For boundary conditions of the kind considered in this work the applicable range for b would be much less and happens to be identical to the one used in this work.
with ' being the reference length. The symbols k 0 and x 0 are the analogues of k and x, 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 k 0 and x 0 are defined, respectively, as and The solution of equation (14) and its associated boundary conditions may then be coupled with the solution of equation (6) given by where jx À x w j L F , andx refers to the x coordinate along the fracture plane. The coupling noted above occurs through equation (12) and equation (13) at the interface X f À X F if we note that Áp f ðx; y w Þ and " qðx 00 Þ 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, Dp 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, Dp 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).

The function f ðsÞ
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 andBoyrchuk (1971), de Swaan-O (1976), Chen (1982), Ozkan and Raghavan (1991), Cinco-Ley and Meng (1988), Noetinger et al. (2001Noetinger et al. ( , 2016 and Raghavan and Chen (2017 ; flow regime in the matrix system reflects transient conditions), and for s ? 0; late times specified as Flow Regime 3, (f ðsÞ ! 1 þ x 0 as tanh(x) % 1/x; flow regime in the matrix system reflects pseudosteady conditions). These three conditions in order of appearance will exhibit slopes of m (early times), m/2 (intermediate times), and m (late times) at the well provided that boundary effects are negligibly small. The term boundary, incidentally, in the context of the transient-flow-model, includes the dashed line shown in Figure 2. Should the boundaries of the fracture system dominate the well response and the flow regime in the matrix system reflect the transient state, then Flow Regime 4, the counterpart of Flow Regime 2 first identified in Raghavan and Ohaeri (1981), will be evident. Flow Regime 4 may be readily identified through its characteristic slope as we discuss below. Finally should pseudosteady conditions prevail everywhere, then the counterpart of Flow Regime 3, namely Flow Regime 5, will be reflected (Table 1).
In concluding this section, we should note that additional options exist to examine solutions; for example, Chen (1982) considers f s ð Þ % 1 þ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi 3x 0s gðg f Þ=k 0 q ; this aspect is not examined here.

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 b < 1; that is, if superdiffusive conditions are to be modeled. Here we use options that are available and examine an instantaneous plane-source-solution in an infinite system which incorporates both sub-and superdiffusive considerations. Let c(x) be the instantaneous source function, located at x 0 , that acts at t = 0, with c(x) = 0, for t < 0.

Plane source
The plane-source solution is particularly useful in obtaining early time responses for a fractured well. Thus, considering a linear system we may derive the plane-source solution along the lines outlined in Carslaw and Jaeger (1959). The required expression is given by In the above expression, E a,b (z) is the two-parameter Mittag-Leffler function (Erdelyi et al., 1955;Mainardi, 2010;Uchaikin, 2013;and Povstenko, 2015): where C represents the complex plane, and the symbol C(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 a,b (z). For b f = 1, the right-hand side of equation (23) reduces to an expression in terms of a stretched exponential given by wherek The right-hand side of equation (25) where R s ð Þ > 0, 1 > 0 as noted in Saxena et al. (2006), or through the Wright function, W(a, b; t), as indicated in Povstenko (2015) by the expression L À1 ½s Àb expðÀks a Þ ¼ t bÀ1 W ðÀa; b; Àkt Àa Þ; The Wright function is defined by

Probability density function
If f(s) = 1, the right-hand side of equation (23) is a probability density function. Denoting Gðjxj; tÞ to be the probability density function, we require that and  To get the integral of the right-hand side of equation (23) and thus on integration of equation (23) the result is: On using the asymptotic expansion in Haubold et al. , jzj ! 1; j arg zj n; 0 <ã < 2 and pã=2 < n < min ½p; pã; for z e C and n is a real number with b < 1, we have Therefore, equation (31) is satisfied. Consequently, the right-hand side of equation (23) is a pdf.
If both a f and b f equal to 1, we have the Gaussian as the inversion of the right-hand side of cðxÞ in equation (23) yields the expression for c(x, t) to be The moments of an even order which would also be of interest are (Uchaikin, 2013).

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, Dp(x, t), in terms of the Laplace transformation for a continuous, plane/ line source at location, M 0 , of strengthqðM 0 ; tÞ=/c is with the symbol S in the case under study representing the surface of the source. IfqðM 0 ; tÞ were constant, then equation (38) simplifies to For the system under consideration, equations (23) and (38) yield the following expression for pressure distribution: whereỹ ! 0 is the distance from the source. The flux distribution in the system corresponding to equation (40) is with The above expressions are also given in Chen and Raghavan (2015). For a bounded system, the pressure distribution, Áp f ðỹ; tÞ, in terms of the Laplace transformation may be obtained from where We may relate the fractional gradient to the pressure drop atỹ ¼ 0 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).

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 X N j¼1 q j ðtÞ ¼ qðtÞ: The assumption of an infinite-conductivity wellbore requires, of course, that the same pressure be maintained at every point inside the horizontal wellbore, namely, that for each fracture location, j = 1, 2, . . ., N we have The convolution equation to compute the pressure response at the wellbore for the multiple-fracture system is for k = 1, 2, . . ., N. Here, jk represents the effect of Fracture j on the pressure response at Fracture k and Áp 0 wjk ðtÞ ¼ dÁp wjk =dt.
In terms of the Laplace transformation, the above three expressions, respectively, are as follows: Áp w ¼ Áp w j¼1;n ; and Thus we are required to solve simultaneously for Dp 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 off-diagonal elements of the matrix would have reflected interactions among fractures. Each Áp wii element may be specified independent of other such elements; consequently, different fracture properties may be incorporated at each location.

Asymptotic solutions
Regardless of the method of evaluating measurements, asymptotic solutions are helpful, even essential, to meet or satisfy a number of objectives. For example, though approximate they provide the following: qualitative information as to whether initial expectations were appropriate, guidance on the appropriate physics that needs to be incorporated in evaluating the response, guidance on the type of the mathematical model to be used in quantitative evaluations, planning of future tests and the like. For example, if we consider a fractured well producing a fractured reservoir, we find that the characteristic responses expected of fractured wells, such as a slope of 1/4 on log-log coordinates during bilinear flow in the classical case, are not evident if transient effects in the matrix significantly influence the well response; see Yeh et al. (1986).
Expressions for responses at a horizontal well producing through multiple fractures that are yet to be presented in the literature are derived in the Appendix for the model considered here. Expressions are derived for all Flow Regimes considered. Table 2 is a summary of the exponents derived in the Appendix; we document results in terms of a f , b f and a m . Specifics of the exponents presented in Table 2 are discussed in the Results section. The exponents given do revert to the classical cases in the literature under identical conditions either for fissured or non-fissured idealizations. In fact, Table 2 reproduces all information available in the literature.

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 b f = 1, equation (53) 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 b 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.
In deriving an asymptotic solution for Flow Regimes 4 and 5, we show in the Appendix that two-term approximations for E a,b (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).

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, a f , b f and a 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 Mittag-Leffler functions are calculated in the manner indicated in Gorenflo et al. (2002).

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 a f , b f and a 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 a f and b f are met. A reduction in b 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 a f (greater restriction). Inspection of the role of the exponents for Flow Regime 2, however, suggests a more complicated picture: For fixed values of a f and b f , the slope increases as a m decreases. But changes in the slope are more complicated for fixed values of a f as b f changes, and depends on the magnitude of a m . If a m were <1, slopes changes depend on the values of a m and a f for fixed values of b f ; then slope changes depend on the values of a m , a f for fixed values of b f ; however, if a m = 1, then for fixed values of a f the slope decreases as b f decreases. The product a m b f in the numerator causes the differences just noted.
If b were < 1, then all other things being equal the exponent (slope of the line on log-log coordinates) will be even much smaller. Chu (2018), who framed his comments in terms of the Chow (1952) Pressure Group, Dp c , has noted that we need to consider situations when slopes are very, very small implying situations beyond subdiffusive effects. And super-connected networks resulting from b f < 1 may 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.
FR 4 1 À a m /2 FR 5, Constant rate 1 FR 5, Constant pressure 1 + a f ; a f < 1 result in slopes (exponents of t) of lines on log-log plots being much less than the slopes based on classical expectations given in Yeh et al. (1986); that is, being < 0.125. While on this subject, we should note that if the values of a f and b f were such that b f þ 1 À 2a f ¼ 0, then the exponents for the linear flow period, Flow Regimes 1 or 3, noted in Table 2 suggest that the well response is identical to the classical case; that is, the slope would be 1/2. Similarly, if b f þ 1 À 2a f À a m b f were equal to 0, then the models of anomalous diffusion for Flow Regime 2 in Table 2 suggest that the trend in behaviors for anomalous diffusion would be similar to those predicted by classical models. Essentially, many points on the (a f À b f ) or (a f À b f -a m ) phase planes yield slopes identical to those of classical diffusion because the points representing classical diffusion, namely, a f of 1, b f of 1, and a 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.

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 and 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 p 0 wD ðt D Þ, (t|dp 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 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 k 0 and x 0 (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 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, p 0 wD , 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).  Table 2. Finally, the computations of the single-fracture solution identified by the circles in Figure 6 again confirm that the single and multiple fracture solutions correlate in the manner indicated by equation (58).

Pressure tests: a and b
We noted in the Introduction that our intent is to consider patterns similar to those in Doe et al. (2014) shown in Figure 1. Figure 7 taken from Raghavan and Chen (2019) is one such case, and offers a plot of flowing and buildup responses in the standard way; that is, pressure and the corresponding derivative at an oil well with the filled-in circles being buildup responses (Dp vs. Dt) and the squares representing flowing responses IRNP(t e ) vs. t e , in terms of the nomenclature used in Rate-Transient-Analysis, RTA. The symbol IRNP(t e ) is the integral of rate-normalized pressure defined in Palacio and Blasingame (1993) as and the symbol t e represents material-balance time defined in Camacho-Velázquez and Raghavan (1989) as (In this example, no other information other than that shown here is available to us as we were asked to show a proof of the concept.) Here we consider issues beyond those identified in Raghavan and Chen (2019) who showed that the flowing and buildup responses correspond to Flow Regime 4 and Flow Regime 2, respectively. The alignment of the flowing and buildup responses is striking and establishes the need to buildup data which is not often measured although early-time flowing-well responses are often quite suspect. The basis for combining buildup and responses is discussed in Raghavan (1980). The slope at late times corresponding to Flow Regime 4 is much larger than what we would expect for classical diffusion and based on the exponent in Table 2 suggests that the value of a m is 0.6. By the same token the slope of the first straight line in Figure 7, namely, 0.1 reflective of Flow Regime 2 is much smaller than expectations based on classical diffusion, 1/4 or 1/8. The estimates provide a well-founded basis for a model of the kind proffered here. The value of a m along with the estimate of 0.1 for the early-time slope may be translated to a range of values for a f and b f . Scouring a map of possible ranges of a f and b f , we arrive at 0.8 a f 1 and 0.1 b f 0.4. This analysis suggests that long-range interactions are important to the performance of this well. Consequently, it behooves the operator to estimate values of a f and b 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.
with a 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. All curves shown here, and others, of course, including the Gaussian may be expressed through the similarity variable, 1,  by the expression where 1 is deduced by considering the arguments of the Mittag-Leffler, Fox and Wright functions on the right-hand sides of equations (24), (27), and (29), respectively, and the solutions for Flow Regime 1 in the Appendix; see equation (A.36). The ability to correlate solutions in terms of 1 is shown in Figure 9 where a typical plot of x D Gðx D ; t D Þ as function of 1 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 a and b.
For a 1 and b = 1, f(1) in equation (63) may be expressed analytically through the Wright function, W(a, b; t), as was noted in equation (28); thus noting that 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 a and b; for regardless of the approach used the following factors are common for predicting the performance of fractured rocks: High and low conductivity paths, discontinuous properties that encompass several orders of magnitude, convoluted flow-paths that serve as conduits and barriers and complex geology (see Graphical Abstract). None of these may be addressed in the classical context. One particular disadvantage of DFN models is that, to our knowledge, at this time no particular rules or guidelines of a general kind that are normally available in evaluating transient responses exist.
Our work provides an option to overcome this limitation. A recent contribution (Artus, 2020) that came to our attention as this work was being completed gives much credence to a model that includes sub-and superdiffusive options of the kind discussed here.

A.1 Asymptotic solutions, development
We first briefly document details pertaining to the finiteconductivity model for a single fracture in a rectangular drainage region and then sketch out the details of the multiple-fracture model in that drainage region. To obtain characteristic responses at early times, namely Flow Regimes 1, 2 and 3, we proceed under the assumption of 1D flow in an infinite reservoir. Modelling, for convenience, one fourth of the drainage region shown in Figure 4, namely x ! x w = 0 and assuming all distances are measured from the well center, equations (6), (12) and (13) subject to the assumptions stated below equation (13) lead to ðA:2Þ For a system produced at a constant rate, q, noting that the well tips are assumed to be sealed, equation (A.1) yields the pressure distribution in the fracture, Áp F ðy D ; sÞ, to be withq F given as 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 The above expression may be readily extended to the situation where a horizontal well is produced through multiple fractures because at early times each fracture acts individually; that is, each fracture acts as if the others do not exist, therefore it is possible to amalgamate responses. Converting equation (A.9) to the situation where the fracture produces at a constant pressure, noting that the total rate for constant-pressure-production may thus be obtained by summing each individual response as all fractures produce against a common well pressure and reconverting the resulting sum to the situation where the horizontal well produces at a constant rate, q, we have where the symbolq F' is now the right-hand side of equation ( ; ðA:14Þ withq F' given as On extracting the first term from the parentheses on the right-hand side of equation (A.15), expanding the resulting expression by the binomial theorem and dropping all but the leading terms we replaceq F' in equation (A.10) by the final result to obtain ðA:16Þ Now on using the expression for f ðsÞ 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, C FD' , to be The symbols p wD (t D ), t D , and C FD' , are defined, respectively, by the following expressions and 1Àa f a f f :

ðA:21Þ
If all three exponents (a f , b f and b 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 tanhq F %q F or tanh(x) % x À x 3 /3 for tanhq F (Camacho-Velázquez, 1984;Camacho-Velázquez et al., 1987). We follow the latter option, and using the right-hand side of equation ( withb f and a, respectively, given as ; ðA:26Þ and a ¼ 1 À a f ð1 Àb f Þ: ðA:27Þ With the result noted in equation (A.22) we may obtain well responses for linear flow, again, on substituting appropriate expressions forkb f and we do so on the assumption that X 2' ¼ 1. On using the right-hand side of equation ( 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 In this section we address our expectations of well responses for the linear-flow regime. As noted earlier, Flow Regime 3 (and also Flow Regime 5, see below) applies for time ranges when the argument of tanh(x) in equation (15) Consequently, the well response, p wD (t D ), is given by ðA:36Þ During Flow Regime 1 the well performs as if the matrix system were nonexistent; consequently, the response during that period may be obtained from equation (A.36) by setting x 0 = 0. In other words, much like the Barenblatt et al. (1960) formulation, responses during Flow Regime 3 parallel that of Flow Regime 1 with the displacement being proportional to ð1 þ x 0 Þb f .

A.1.5 Flow Regime 4 and Flow Regime 5
We consider responses to a single fracture for times that are long enough; also, if times are long enough the influence of the conductivity of the fracture may be ignored. The needed expressions are derived directly from the development given here rather than the indirect approach used in Raghavan and Chen (2018a). We then compare the derived expression with that obtained for subdiffusion. As subdiffusive solutions obtained are exact, again as in the text, parallels are drawn to the expression for the pressure distribution in 2D.
As the solutions involve the ratio of the Mittag-Leffler functions, for small values of s, the ratio of Mittag-Leffler Functions, E a 1 ;b 1 z ð Þ and E a 2 ;b 2 z ð Þ, may be written as