An application of a multiindex, time-fractional differential equation to evaluate heterogeneous, fractured rocks

. A multiindex, distributed fractional differential equation is derived and solved in terms of the Laplace transformation. Potential applications of the proposed model include the study of ﬂ uid ﬂ ow in heterogeneous rocks, the examination of bimodal ﬂ uid exchange between mobile-immobile regions in groundwater systems, the incorporation of the existence of liesegang bands in fractured rocks, and addressing the in ﬂ uences of faulted and other skin regions at interfaces, among others. Asymptotic solutions that reveal the structure of the resulting solutions are presented; in addition, they provide for ensuring the accuracy of the numerical computations. Fractional ﬂ ux laws based on Continuous Time Random Walks (CTRW) serve as a linchpin to account for complex geological considerations that arise in the ﬂ ow of ﬂ uids in heterogeneous rocks. Results are intended to be applied at the Theis scale when combined with geological/geophysical models and production statistics to all aspects of subsurface ﬂ ow: production of geothermal and hydrocarbon ﬂ uids, injection of ﬂ uids into aquifers, geologic sequestration and hazardous waste disposal. Results may be extended to study the role of complex wellbores such as horizontal and fractured wells and more complex geological considerations such as faulted systems.

Ratio of fracture to total volumes; see equation (12) Subscripts D Dimensionless f Fracture system 1 Introduction Time-Fractional Flow Equations (t-FFEs) are particularly useful in the evaluation of systems of complex geology to incorporate the internal architecture and internal imprint involving barriers, intercalations and the like Chu et al. (2017Chu et al. ( , 2018, Ingle et al. (2020), Raghavan and Chen (2018). Such blockages result in interconnected highconductivity pathways interspersed with low permeability blockages. t-FFEs represent models incorporating sub-diffusive propagation of the pressure fronts at the Theis scale and result in temporally non-Darcy, non-stationary flow to account for the complex circuitry. As these models are based on motions that depend on a Continuous Time Random Walks (CTRW) (or a waiting-time), such models yield power-law declines in flux (Raghavan and Chen, 2017a;Xia et al., 2021). Examples of subdiffusive behaviors are reported in a number of sources in the earth sciences (Bisdom et al., 2016;Caine et al., 1996;Chu et al., 2017Chu et al., , 2018Cortis and Knudby, 2006;Doe, 1991;Evans, 1988;Jourde et al., 2002;Knudby et al., 2002;Mitchell and Faulkner, 2009;Noetinger and Estebenet, 2000;Noetinger et al., 2001Noetinger et al., , 2016Raghavan and Chen, 2018;Savage and Brodsky, 2011;Scholz et al., 1993;Suzuki et al., 2016;Xia et al., 2021). This paper examines applications of distributive fractional differential equations in a manner similar to those discussed in Luchko (2011) in two-porosity rocks with high and low permeability along the lines in Barenblatt et al. (1960). Our choice is dictated because of the model's ubiquitousness and its ability to apply to situations beyond fractured rocks such as the Mobile-IMmobile-zones (MIM) model of Silva et al. (2009) which is similar to the Multi Rate Mass Transfer (MRMT) model discussed in Haggerty et al. (2000); see Xia et al. (2021). In MIM models the volume of the region with higher connectivity is much larger than the one with lower connectivity. Solutions applicable to the evaluation of single and multiple (interference, tomographic) well tests are presented. The motivation for this study is to evaluate the consequences of finite speeds of propagation of the transient in complex systems of the form where conventional representations of heterogeneity through log-normal representations of permeability are inadequate (Xia et al., 2021). Two exponents that reflect the heterogeneity in permeability and denoted by the symbol a j are used to study the rate of propagation of pressure fronts. These exponents reflect the random walk dimension of fractal objects (Metzler et al., 1994).
2 Problem statement, flux law, and differential equations We consider the classic problem of the flow of a slightly compressible, constant viscosity fluid to a cylindrical wellbore subject to the idealization in Barenblatt et al. (1960) and Warren and Root (1963). Initially the pressure, p i , is assumed to be constant everywhere. Fluid compressibility and viscosity are denoted by the symbols c and l, respectively. As noted in the Introduction, flow in the porous rock of constant porosity, /, and thickness, h, is governed by the CTRW process through two exponents a f and a m where the subscripts f and m stand for the fissure and matrix systems, respectively, in which case the volume of the region corresponding to m is assumed to be much larger than that of f. The opposite would be true for the MIM model. This nomenclature is used merely for convenience and because of familiarity; for example, in the mobile-immobile context the subscripts f and m represent the mobile and immobile regions, respectively. The center of the well is assumed to be the origin with the system assumed to be infinite in its areal extent and as is usual we work in the cylindrical coordinate system through the Laplace transformation. Subdiffusive flow is represented through a fractional derivative defined by Caputo (1967) that may be deduced through the fractional engine discussed in Metzler and Klafter (2000). Detailed and elaborate citations that include theoretical considerations and experimental observations of the fractional flux law may be found in Raghavan and Chen (2020). This flux law is also particularly suitable for considering problems in 2D such as those considered in Beier (1994). The well is assumed to penetrate the rock completely and we consider constant-rate-production at a rate q; wellbore storage effects are not considered.

The Flux law
The flux law under a CTRW derived in Montroll and Weiss (1965), Noetinger and Estebenet (2000) and Henry et al. (2010) used in this study is given by where the symbol p(x, t) represents fluid pressure, and where o a f(t)/ot a is the temporal fractional derivative defined in the Caputo (1967) with CðÁÞ being the Gamma function. Equation (1) implies that under CTRW or subdiffusive flow, the flux at time, t, is not directly obtainable from the instantaneous pressure distribution, p(x, t). The flux law presumes an independent, identically distributed random variable for the step length, Dx, and a waiting time function, w(t), for the time between jumps. Many such time functions have been discussed; see, for example, Feller (1971), Kenkre et al.(1973) and Hilfer and Anton (1995). Expressions similar to equation (1) have been used to describe diffusion in fractal objects, contaminant transport, heat transmission, flow in porous media and for biophysical applications. The exponent a in the expression of the fractional derivative follows a probability density function, Gðr; tÞ, that exhibits power-law behavior at long enough times for a < 1 given by Schneider and Wyss (1989): Here d is the number of dimensions and agrees with the expressions in Carslaw and Jaeger (1959) for d 3 when a = 1.

The differential equations
On combining the equation for the conservation of mass with the flux law and the equation of state for a slightly compressible liquid, the Barenblatt et al. (1960) expressions that govern flow in a fractured rock under CTRW are: for the fissure system and for the matrix system by or, in other words, Here Dp is the pressure drop, p i À pðr; tÞ, r is the shape factor, ' is the reference length and the subscript D reflects nondimensionalization of the space variable with respect to '. On noting that the Laplace transformation of the Caputo fractional operator, C 0 D a t f ðtÞ, in terms of the Laplace variable, s, is given by Equations (5) and (7), respectively, in terms of the Laplace transformation, become whereg f , the diffusivity of the fracture system is based on fracture pore volume; that is, The bar symbol denotes the Laplace transformation. On now combining equations (9) and (10), we have and Here and The variables x and k, the storativity ratio and dimensionless transfer coefficient, respectively, are the wellknown standard variables used to describe naturally fissured rocks; in this case k is defined for subdiffusive flow.

Solution
As is well known, the principal advantage of expressing the differential equation in the form noted is that if f(s) were set to unity, then the expression is identical to that for flow in a porous rock for subdiffusive systems and, further, if a f were equal to unity then the differential equation reduces to that of a homogeneous rock considered in the literature in great detail. Thus it is immediately obvious that the solution for flow to a well in the form of a cylindrical system is given as (Theis, 1935) where, Dp f (r, t) = [p i À p f (r, t)] is the magnitude of the drop in pressure below the initial pressure, p i . The symbols I 0 (x) and K 0 (x) represent the modified Bessel functions of the first and second kind of order 0, respectively. The kernel, u, is given by Equation (19), of course, presumes that the initial condition is The simple consequence of using s a f f ðsÞ in the kernels is that the solutions in Raghavan and Chen (2020) that were derived for the single-porosity case may be used directly to write expressions for the pressure distributions when pressures are governed primarily by the fissures (at early times) and as a result of both the fissures and matrix (total system) at later times. For example, we note that for early times (s ? +1), the faster leg is dominant reflecting the fracture pore volume and we have f ðs ! þ1Þ ¼ 1; at intermediate times kx=½ð1 À xÞs am gðg f Þ < 1 < k=½ð1 À xÞs am gðg f Þ; ð23Þ the lower-transmissivity, slower-region plays a significant role, and under such circumstances we have and at long enough times (s ! 0þ), the faster region again dominates but reflects the total pore volume and The approximation given in equation (24) results in pressure distributions that are yet to be reported in the literature. As we see both exponents, a f and a m play a role in the transitional period from one behaving like a system governed only by the fast leg, equation (22), to one governed by both regions, equation (25). As we will see the responses are distinctly different from those of the classical case. This type of coupling is not evident in Xia et al. (2021); the two legs do not operate in parallel. We refer to the three approximations given in equations (22), (24) and (25) as Flow Regime 1, Flow Regime 2 and Flow Regime 3, respectively. We now seek a solution to equation (19) in an infinitely large system subject to the following boundary conditions: Equation (27) specifies the condition for a well of finite radius, r w , which we take to be the reference length; thus the dimensionless distance r D is given as The pressure distribution, Áp f r D ; s ð Þ, subject to the above conditions may be obtained in a straight forward manner to be If we are to understand the structure of the solutions, however, we need to approximate equation (29) by considering a line-source well; that is, assume x K 1 (x) ? 1 as x ? 0 and in this case we may derive solutions for the three f(s) approximations considered.

Solutions for the situation xK 1 (x) as x ? 0
Under this assumption, we may write equation (29) as To generalize and explore approximate solutions subject to the line-source approximation, we define dimensionless pressure, p Df (r D , t Df ), and dimensionless time, t Df , respectively, by and t Df ¼g With the asymptotic representation of K 0 (z) as z ? + 0 given by Carslaw and Jaeger (1959) where c is Euler's constant, and the relationship in Oberhettinger and Badii (1973) given by L À1 s Àm ln s ð Þ¼½CðmÞ À1 t mÀ1 ½wðmÞ À ln t; where w(Á) is the Digamma function and L is the Laplace transformation, the dimensionless pressures for the three approximate expressions of f(s), equations (22), (24) and (25) are, respectively, given by for early times, for intermediate times, and for long times.
As noted earlier the structure of the solution during Flow Regime 2 is significantly different from that reported in the literature; Streltsova-Adams (1978) now both a f and a m influence pressure. In addition, we see that during Flow Regime 3 two time scales, namely, those involving the fracture volume and the total volume affect the response with the fracture volume being, perhaps, more dominant. As expected, all three expressions given above reduce to the classical case for a f = 1 and a m = 1. Also note that as x ? 1, the role of Flow Regime 2 diminishes. The derivative or logarithmic derivative, p 0 D ðr D ; t Df Þ, corresponding to equations (35), (36) and (37) are, respectively, and p 0 D ðr D ; t Df Þ ¼ ð40Þ Again for a f and a m = 1, these results correspond to the semi-logarithmic approximation of the Theis (1935) solution, namely the Cooper and Jacob (1946) approximation. In concluding this discussion it must be clear that two requirements are to be satisfied. The value of the Laplace variable, s, needs to such that those required of the approximations of f(s) hold and those required for s to be small enough for the logarithmic approximation to apply also hold.

Computational results
The principal goal is to illustrate the nature of responses for subdiffusive flows in distributed time-fractional flow systems and illustrate the consequences of subdiffusive flow by elucidating the roles of the exponents, a f and a m . The interaction between the exponents, a f and a m , produces a wide range of responses and we document the range of expectations. Both pressure and derivative responses of the rigorous solution are compared with the asymptotic expressions derived here in order to discern the structure of the solutions. Our focus is on Flow Regime 2. The roles played by the constants k and x are similar to those for the classical solutions and are not addressed in any detail as they are well documented in Warren and Root (1963) and Kazemi (1969). All solutions of equation (19) and associated boundary conditions, equations (26) and (27), are computed in the usual way by the Stehfest algorithm (1970aStehfest algorithm ( , 1970b. The examples presented provide evidence for both the viability of the model and the accuracy of the computations. Figure 1 presents pressure responses at the wellbore, p wD ðt Df Þ, for the simplest system one may envisage, a f = 1 with a m being the parameter of interest; three values of a m are considered. The complement of these results, namely the corresponding derivative curves, p 0 wD ðt Df Þ, are shown in Figure 2. In both figures responses for the classical case are included to serve as a backdrop. For a f = 1 as noted in equations (35) and (37) straight lines with slopes of 1.151 on semilog coordinates will be evident for Flow Regimes 1 and 3, respectively, and this point is confirmed in Figure 1. During intermediate times, Flow Regime 2, the solution for a m = 1 is flat as suggested in equation (36) whereas the solutions for a m < 1 are distinctly different displaying broad sweeping patterns. As shown in Figure 2, the magnitude of the responses becomes muted as a m decreases; further, the duration for which Flow Regime 2 exists increases. Both these observations essentially reflect the reduction in effective transmissivity for lower values of a m .
The results in Figure 3 consider the situation when a f < 1. Results for both dimensionless pressure, p wD ðt Df Þ, and its derivative, p 0 wD ðt Df Þ, are shown on log-log coordinates to account for the existence of the t f ða f Þ Df term in equations (35) and (37). The circles and the squares display the duration for which Flow Regimes 1 and 3, respectively, exist. It is striking that the transitional period like that shown in Figure 2 during which Flow Regime 2 would be evident is entirely nonexistent. Conditions for which Flow Regime 2 would exist are not readily apparent for t-FFEs describing subdiffusive flow in two-porosity rocks when a f is much different from a m . Most important to note here is that the time of transition from Flow Regime 1 to Flow Regime 3 is short and power-law responses with slopes reflecting the exponent of t Df , namely, (1 À a f )/a f , are not evident. This result is markedly different from what we see when we consider Cartesian coordinates; see Raghavan (2012) who showed that at early times the well response to a fractured well in Cartesian coordinates is given by Chu et al. (2017Chu et al. ( , 2018 have used equation (41) to identify subdiffusive systems in convincing, compelling ways. Figure 4 explores the transitional period further and two situations where a f ? a m are considered. Results that are complements of the ones discussed in Figures 1 and 3 are considered here. The top curve represents an a f of 0.75 and a m of 0.6 and this response may be compared to the derivative curve in Figure 3. Unlike that situation, here Flow Regime 2 dominates the responses shown for several log cycles; the other two Flow Regimes are not evident. The responses for the bottom curve may be compared to the situation where a f and a m are equal in Figure 2. Although muted, here again, we note the existence of a ushaped derivative curve with a rather flat bottom; the circles corresponding to Flow Regime 2 do yield ð1 À a f Þ=a f because the logarithmic term in equation (39)

Discussion and concluding remarks
A multiindex, distributed fractional differential equation is developed to examine flow in rocks that may be characterized by two-conductivity systems, one a heterogeneous system consisting of dead ends, intercalations and the like wherein fluid movement occurs rapidly and the second, a slower forcing system which is also presumed to be heterogeneous. The interfaces between these regions are not necessarily distinct and may contain entities such as mineral deposits, skin regions, etc. The net effect is that pressure propagation is delayed and other complex phenomena in the form of flow reversals etc. need to be addressed; Raghavan (2004). A two-index time-fractional, subdiffusive system offers a potential path to address such situations. The time-fractional equation follows the outlines of the Barenblatt et al. (1960) model as it is a convenient vehicle for presenting the essential elements we desire; furthermore, it is both ubiquitous and well understood. This model opens up a wealth of options for examining nonlocal phenomena where diffusion is central; it is intended for situations where the geology is complex in fields such as the evaluation of geothermal systems, planning of contaminant transport, ascertaining waste disposal options and the like. Unlike conventional formulations, Gaussian statistics is not obeyed. Situations where this is the case are noted in the Introduction. Unlike DFN formulations, wherein the rules of Gaussian statistics are observed and details of smaller subunits are needed to predict performance; fractional diffusion assumes that such details may be dispensed with. It has occurred to us that incorporation of fractional diffusion within the subunits would be a useful exploration.
The principal contribution of this work is that subdiffusive features are introduced in a rigorous way for the first time. Three specific flow periods may be identified during transient flow and the intermediate flow period may be absent when a f is much different from a m . When evident, the characteristics of this period may be used to identify subdiffusive flow as it is much different from the classical case. The duration of this period increases as the magnitude of a m decreases all other things being constant and as x ? 1.
Potential applications of t-FFE models include fractured rocks where heterogeneity such as variations in permeability, variations and/or the width of fissure apertures are dominant, existence of liesegang bands (mineral deposits abutting fissures) and complex interconnectedness (Belayneh et al., 2006;Cacas et al., 2001;Kang et al., 2014;Reiss, 1980;Sharp et al., 1996), lower-conductivity skin regions in sandstones (Fu et al., 1994, solute retention, Zhang et al., 2020 and "water pile-up" in bimodal transport that results in blockages at low permeability regions and consequent early arrivals in regions of higher permeability (Xia et al., 2021). Nonlocal pressure fronts of the kind discussed here may have to be incorporated into field studies even if other information such as solute transport suggests otherwise (Zhang et al., 2020). As the model is developed in terms of sources and sinks and the Laplace transformation, numerical implementation is straight forward and extensions to other well configurations and boundary conditions is also straight forward (Raghavan and Chen, 2017b).
We do emphasize that tests must be supplemented with information on reservoir architecture (geological model, maps, logs, cores, etc.) for a complete analysis.