Regular Article
Subdiffusive flow in a composite medium with a communicating (absorbing) interface
^{1}
The Raghavan Group, PO Box 52756, Tulsa, OK 74152, USA
^{2}
Kappa Engineering, 11767 Katy Freeway, #500, Houston, TX 77079, USA
^{*} Corresponding author: raghavan.raj@gmail.com
Received:
23
October
2019
Accepted:
2
March
2020
Twodimensional subdiffusion in media separated by a partially communicating interface is considered. Starting with the appropriate Green’s functions, solutions are developed in terms of the Laplace transformation reflecting two circumstances at the interface: situations where there is perfect contact and situations where the interface offers a resistance. Asymptotic solutions are derived; limiting forms of the expressions reduce to known solutions for both classical diffusion and subdiffusion. Specifics are analyzed in depth with reference to flow in porous media with potential applications to the evaluation of the role of subsurface faults and flow in fractured rocks. Characteristics of the derivative responses are documented extensively as they are the linchpin for evaluation of pressure tests. Results given here may be used for evaluation at the Theis (1935; Eos Trans. AGU 2, 519–524) scale along with geological and geophysical properties, and production statistics. Yet a subdiffusive model does not imply a single value for properties. The method presented here may be extended to multiple contiguous media and to subdiffusive transport in many contexts (complex wellbores such as inclined, fractured and horizontal wells, situations such as sequestration, production of geothermal systems, etc.).
© R. Raghavan & C.C. Chen, published by IFP Energies nouvelles, 2020
This is an Open Access article distributed under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
Nomenclature
Aj: Constants; see equations (28)–(30)
a: Defined in equation (7)
B: Formation volume factor [L^{3}/L^{3}]
ct: Total compressibility [L T^{2}]
Eα,β (z): MittagLeffler function; see equation (9)
f_{j} (φ): See equations (18) and (19)
Kν(z): Modified Bessel function; second kind of order ν
: Subdiffusive permeability; see equation (11) [L^{2}/T^{1−α}]
L: Distance between the interface and the well [L]
p_{i} : Initial pressure [M/L/T^{2}]
p′: Logarithmic derivative [M/L/T^{2}]
s: Laplace transform variable [1/T]
T_{D} : See equation (46)
Tr: See equation (31)
α: Exponent; see equation (10)
γ_{j} (x; t): Source functions; equations (18) and (19)
: “Diffusivity”; see equation (20) [L^{2}/T^{α}]
η_{D} : See equation (48)
: Equation (11)
: Contact resistance; see equation (16)
: See equation (47)
: See equation (32)
Subscripts
j: Zone with distinct properties
Superscripts
1 Introduction
Interfaces are ubiquitous in many fields of endeavor; extensive literature exists in a range of technologies (biophysics, geophysics, energy recovery, tomography, acoustics and electromagnetism to mention a few). Developments in heat transmission made through Green’s function techniques are often the starting point in many studies; see Carslaw and Jaeger (1959), in both one and two dimensions for single and multilayered media. Often point or linesource solutions are employed as in Mandelis et al. (2001) and Shendeleva (2004). More recently, subdiffusive flow in linear, composite media have been explored (Kosztoowicz, 2017; Povstenko, 2013; Povstenko and Kyrylych, 2019; Raghavan and Chen, 2018). The motivation for such studies, as in all other cases, is to evaluate the consequences of finite speeds of propagation of the transient in linear, composite systems of the form
(1)with being the distance, being time and where is a constant. Equation (1) for , namely subdiffusion, has drawn much interest in recent times because of diffusion in fractal structures (Gefen et al., 1983; Metzler et al., 1994; Nigmatullin, 1984a, b), although interest in finite speeds of propagation has been discussed for a long time; see Gurtin and Pipkin (1968). More recent examples of subdiffusion from an application vantage point particularly from the perspective of flow in porous media, the focus of this work, may be found in Caputo (1999), Jourde et al. (2002), Cortis and Knudby (2006) and Raghavan (2011). In this work we discuss influence of an interface under subdiffusion in 2D as no such studies exist. We examine the problem in the context of application to earth sciences such as the flow of groundwater and extraction of hydrocarbons by considering the interface in the form of a partial hydrologic barrier. Partial hydrologic barriers are of particular interest as they could serve as conduits at times or seals in other times (Bense and Person, 2006; Miller, 2005; Raghavan and Ozkan, 2011). There is evidence in the earth science literature that the permeability changes significantly around faults and that damaged regions with cracks, fissures and cervices abound around most faults. As a result, scaling laws that systematically reflect the damage around faults have been developed (Caine et al., 1996; Evans, 1988; Mitchell and Faulkner, 2009; Savage and Brodsky, 2011; Scholz et al., 1993). Liesegang bands, regions of reduced permeability around fracturematrix interfaces, are modeled in classical situations as skin regions (Prats and Raghavan, 2013, 2014) (ignoring scaling laws), and the phenomenon is known as interporosity skin (CincoLey et al., 1985); photographs of Liesegang bands may be found in Fu et al. (1994), Sharp et al. (1996) and Neville et al. (1998). As we employ the method of sources and sinks, our results may be directly transposed to apply to other disciplines by a change in nomenclature. The source functions to determine the needed solutions are constructed in the manner of Raghavan and Ozkan (1994) in terms of the Laplace transformation with the method proposed in Sommerfeld (1909). The effectiveness of his method is wellattested by the number of researchers who have employed Sommerfeld’s technique in a wide range of disciplines such as electromagnetic radiation (Lindell and Alanen, 1984; Lindell et al., 1986), geophysics (Hänninen et al., 2002) and flow in porous media (Raghavan, 2010) in addition to heat conduction (Shendeleva, 2004).
We presume two options to represent conditions at the interface: perfect contact or a steadystate resistance as modeled in Carslaw and Jaeger (1959). The latter observation implies that the resistance either causes or produces no transients. We may readily expand the solutions derived for other well configurations if desired; see Raghavan (2012) or address anisotropy as we arrive at expressions for the instantaneous source function, , very similar in structure to that used to represent various forms of linear barriers with image wells of the form described in Tualle et al. (2000):
(2)where the symbol S represents the source, the symbol S′ represents the image sources that are distributed to reflect the interfaces, the symbol Π is the solution for a source placed in an infinite medium with the properties of the region where the source S is located and the symbol ⊛ represents convolution in both space and time.
1.1 The Mathematical model
Transient diffusion in a porous rock of thickness, h, structured in the form of a composite system with an interface located along the plane y = 0 and extending to infinity in both the x and y directions is considered. The influence of the distinct properties of the rock in a connected domain, Ω, in on either side of the interface is examined. A contact resistance at the interface governs fluid movement and the system is produced through a linesource well located at the point (0, y′). In and of itself the interface has no storage capacity; consequently no transient effects are produced as the result of its existence. As is usual we consider a slightly compressible liquid of compressibility, c, and viscosity μ. Initially, throughout the system the pressure, , is assumed to be uniform at the value of . The transient diffusion equation for this configuration is
(3)where the symbols and represent the flux and porosity, respectively. This expression is solved for the following constraints: Because the system is quiescent initially and also infinite in areal extent we require,
Assuming flow to a linesource well, the condition at a wellbore for a specified withdrawal rate, q, is
(6)where B is the formation volume factor, and r is the distance between the point (x, y) and the well location (0, y′).
To satisfy the conditions at the interface, y = 0, we require that the pressure drop across the interface be proportional to the normal flux and also the normal fluxes be continuous both for all time, t (Yaxley, 1987); thus, we require that
(8)where the symbol is a constant and the subscript “F” signifies the fault. If “a” were 0, then the condition specified in equation (7) reflects perfect contact and thus the continuity in pressure for all time, t. In essence, steadystate Darcy flow is assumed across the fault and the fault causes no transients in of itself as it contains no pore volume. The subscripts 1 and 2 in the above expressions represent the regions on either side of the interface; henceforth, we presume Zone 1 represents the region where the well is located. We have treated the existence of the discontinuity at y = 0 in a manner very similar to that of a resistance discussed in Carslaw and Jaeger (1959).
In the following we work in terms of the Laplace transformation and use the method of sources and sinks; see Raghavan and Ozkan (1994). As already mentioned, the influence of the interface is incorporated by image sources located symmetrically about the interface with the resulting pressure distribution determined through convolution.
1.2 The Darcy law, subdiffusive flow
The experiments that Darcy (1856) conducted to arrive at the flux law that the superficial velocity is proportional to the potential difference by considering the flow of water in porous columns may be derived from general energy considerations as shown in Miller (1954). For his experiments, a number of considerations are implicit as set out in Childs (1969). Specifically, as noted in Raghavan (1993) five factors are inherent: (i) the porous medium is an uniform body as it is appreciably larger than the pores that constitute the porous rock, (ii) the structure of the porous rock and the potentials considered are notional surfaces; in essence, the different crosssectional elements that permit fluid flow are so alike that the fluxes are identical and the potential surfaces reflect the overall flux, (iii) the flux is inversely proportional to the length of the porous column, (iv) motion is independent of time (Philip, 1957), and (v) inertia terms are nonzero but are of little consequence (Lindquist, 1933). In the context of random walks, the flux law such as the Darcy law is equivalent to a particle travelling a fixed distance, Δx, over a fixed time, Δt.
Should the structure of the porous medium be considered to be a fractal porous object (Gefen et al., 1983; O’Shaughnessy and Procaccia, 1985a, b), then the nature of the flux law changes; both permeability and porosity are powerlaw functions of the Euclidean dimension, d, the fractal or Hausdorff dimension d_{f}, and the random walk dimension or anomalous diffusion coefficient d_{w} (Chang and Yortsos, 1990) because the diffusivity of the porous rock, η, follows the relation , where r is distance; see Gefen et al. (1983). Other fractal models such as those proposed in Metzler et al. (1994), though never specifically noted, require a different flux law. The easiest way to model their work is to consider a time dependent diffusivity which leads to the fractional differential equation proposed by them.
Anomalous diffusion, the focus of this work, is best described in terms of a general random walk, namely a ContinuousTime Random Walk (CTRW) proposed by Montroll and Weiss (1965); also see Noetinger and Estebenet (2000) and Henry et al. (2010). Unlike standard random walks that use a fixed step length, Δx, at intervals that are separated through a fixed time interval, Δt, the CTRW formulation uses an i.i.d. (independent, identically distributed) random variable for increments of the step length, Δx, and a waiting (pausing) time function, ψ(t), that determines the time between jumps. Kenkre et al. (1973) show that a CTRW formulation leads to an equation (part of the family of master equations) that is equivalent to the Montroll–Weiss formulation. Hilfer and Anton (1995) connect the CTRW formulation to fractional differential equations through the waiting time function , where is the MittagLeffler function discussed in many texts (Erdelyi et al., 1955; Mainardi, 2010; Povstenko, 2015; Uchaikin, 2013):
(9)where C represents the complex plane and where is the Gamma function. The symbol ω is a constant that depends on the properties of the fractal object; with 0 ≤ ω ≤ 1. In a similar manner, Henry et al. (2010) demonstrate that if we were to use two probability density (pdfs) functions, one a step length density function, Δx, that is an even function and the other a powerlaw waiting time function, ψ(t), that is scale invariant with temporal memory, then such an option also leads to a fractional flux law. For instance, the Pareto waiting time density (see Feller, 1971) with a mean waiting time that is infinite could be used for, .
The CTRW formulation leads to a flux law of the following form:
Here α is the subdiffusion coefficient with 0 < α < 1, μ is the fluid viscosity, is the permeability corresponding to α and , is the fractional derivative defined in Caputo (1967) sense; thus
For α = 1, equation (10) reverts to the classical Darcy equation. The Laplace transformation of the Caputo fractional operator, , is given by
Implicit in this definition is the fact that the flux at any instant, t, may not be obtained directly from the instantaneous pressure distribution, p(x, t).
For subdiffusive flows, fractional flux laws of the form proposed here have been used for many years in a variety of contexts concerning diffusion: in biophysics (Magin et al., 2013), fractal structures (Dassas and Duby, 1995; Metzler et al., 1994), flow in porous media (Albinali et al., 2016; Lẽ Mehautẽ and Crepy, 1983; Nigmatullin, 1984a, b; Płociniczak, 2015; Raghavan, 2011; Raghavan and Chen, 2017; Su, 2014; Su et al., 2015), contaminant transport (Fomin et al., 2011; Kim et al., 2015; Molz et al., 2002; Suzuki et al., 2016), heat transmission (Angulo et al., 2000; Gurtin and Pipkin, 1968; Moodie and Tait, 1983; Norwood, 1972). Experimental aspects are discussed in Tao et al. (2016), Zhokh and Strizhak (2018, 2019) and Yanga et al. (2019). Additional options are available to modify equation (10) should one address phenomena such as truncation (Odling et al., 1999) to model structures that are discontinuous resulting in the variation of permeability around faults. Also scalings that result in power law variations yield cutoffs which may result in changes in flow regimes (classical to subdiffusion). In this work we, however, limit our discussion to equation (10).
While the subject of pdfs, we simply note that a standard random walk may be modeled by considering the Markovian exponential waiting time density (memorylessness; see Feller, 1971).
1.3 Formulation
Using , where the subscript j = 1, 2 reflects Zone 1 and Zone 2, respectively, we may combine equations (3) and (10) to obtain the transient diffusion partial differential equations for subdiffusive flow in Zones 1 and 2, respectively, to be
Similarly, equations (6)–(8) lead, respectively, to
1.4 General solution
Instantaneous line sources of the kind needed in this study may be constructed by choosing expressions to satisfy the basic partial differential equation and the auxiliary conditions; for line sources considered in this work one such expression is the modified Bessel function of the second kind of order 0, K_{0} (x), that reflects the free subdiffusive propagator, , in . Accordingly, for Zone 1 and Zone 2 we write instantaneous solutions, γ_{j} (x, y); j = 1, 2, respectively, in terms of their Laplace transformation, with the source located at () to be
(19)where , the “diffusivity” of the medium, is given by
(21)with the symbol q_{j} given by
The motivation for the above expressions is the observation given in Erdélyi et al. (1954, see p. 17); Gradshteyn and Ryzhik (1994, p. 532; 3.961–3.9612) that the integral representation of is
In addition, the first term is the Green’s function for a medium that is infinite in lateral extent. Our procedure, first proposed by Sommerfeld (1909) has been used by many; see, for example, Mandelis et al. (2001), Shendeleva (2004) and Raghavan (2010).
On using equations (16) and (17) to solve for f_{1} (φ) and f_{2}(φ), we obtain the following expressions:
The functions , and are given, respectively, by
Here
(32)and with the source functions on hand we may write the pressure distributions through the procedure followed in Raghavan and Ozkan (1994). The pressure in terms of the Laplace transform of Δp (x, y, z; t) (Raghavan and Ozkan, 1994) is
Here (x′, y′) is the location of the source, (x, y) any location in the porous rock and the strength of the continuous linesource is with dS′ being the element (a line or surface) through which fluid is extracted. For an instantaneous linesource with the source strength being constant, we have
Thus on using equation (34) and the expressions for derived above, the expressions for in terms of their respective Laplace transformations, , are
1.5 Solutions for perfect contact
To proceed further, we define dimensionless variables for pressure, p_{D} (x _{ D }; t_{D}), time, t_{D}, and distance, x _{ D }, respectively, as follows:
(40)where is the reference length.
Using these definitions, we obtain the following expression for the pressure distributions, ; j = 1 and 2, in terms of the Laplace transformation for Zones 1 and 2 to be, respectively,
(42)for . If we consider the integrals in equations (36) and (37), we see that the response at any point M for the interface located at a distance, L from the well is a function of , and T_{r} where
But the A_{j}’s; that is, both q_{j}’s, and the three variables, T_{r}, , and , must, of course, be rescaled. The scalings are:
Scalings involving the ratio of α_{j}’s enter into the picture because we consider subdiffusive flow; in accordance with equation (10); the dimension of is a function of α. These scalings also follow if the partial differential equations are considered. The groupings T_{D} and η_{D} may be replaced by T_{D} and if desired, where
A number of approximate solutions may be derived if q_{1} = q_{2}, which implies and and also α_{1} = α_{2}. In the following we use , where r_{w} is the well radius. In this situation for Zone 1 and Zone 2, respectively, we have
The righthand sides of equations (51) and (52) may be inverted with the expression
(53)where the symbol is the Fox function or the Hfunction, , , as noted in Saxena et al. (2006), and thus we have the following result for , namely
Although many discussions of Fox functions are available; see Mathai and Saxena (1978) and Mainardi et al. (2005), such functions are best evaluated through packages such as Mathematica. Also, if j = 0; that is, one well in a system that is infinite in its extent, equation (55) reduces to the expression, as expected, in Raghavan and Chen (2019).
Along the same lines, we may obtain a longtime approximation for this specific situation by considering the following expression for as x → 0 (implicitly we consider situations where s → 0), namely using
(56)where γ is Euler’s constant. The longtime pressure response, p_{D} (x _{ D }; t_{D}→∞), may be obtained by using the first two terms of the above expression, as well as the result of Oberhettinger and Badii (1973)
(57)where ψ(⋅) is the Digamma function, and we therefore have
Also, the logarithmic derivative, , is given by the expression
Focussing on the first term of equation (58) for now and assuming T_{D} = 1, we note that for evaluating fractured reservoirs, interestingly, Bernard et al. (2006) propose that responses be evaluated as a function of t ln t; this work perhaps gives a foundation for their empirical recommendation. Both expressions, equations (55) and (58) for α = 1 in an infinitely large system (T_{D} = 1), respectively, correspond to the Theis (1935) solution and the semilogarithmic approximation of the Theis (1935) solution, namely the Cooper and Jacob (1946) approximation. Along the same lines if both α_{j}’s are unity, then both these expressions correspond to the ones given in Raghavan (2010). Most importantly, it is to be noted that unlike classical diffusion, for subdiffusive flows the influence of Zone 1 never stabilizes; that is, one must consider transient effects caused by Zone 1 for all times. This difference is exceedingly important and is a direct result of the observation made earlier that the flux, v(x, t), is not be obtained from the instantaneous pressure distribution, p(x, t).
At long enough times the slope of the second line would reflect 2/(1 + T_{D}), if both lines corresponding to Zones 1 and 2 were evident, which would require values of L to be neither too small nor too large, where L is the distance between the interface and well. Then (1 − T_{D})/(1 + T_{D}) may be estimated by extending the line corresponding to Zone 1 and determining the separation between the two lines.
Unfortunately, it does not appear to be possible to extend the above results for different from at this time. Although this is the case, we know that for classical diffusion () the influence of is small and equations (58) and (59) hold in most cases when the α_{j}’s are 1; see Bixel et al. (1963). That being so, it is plausible that equations (58) and (59) may hold for provided that the ’s were equal. This, of course, would have to be explored.
2 Results
Our main concern in the following is the understanding of the characteristics of the derivative curve, , as these curves best demonstrate the characteristics of the rock under consideration. For each of the solutions presented, the pressure response, , is also computed and examined. Our focus is responses at wells on the same side of the interface. Simulations shown were obtained through equation (36) through the aid of the Stehfest algorithm (1970a, b) as is normal for most welltest problems.
2.1 Perfect contact
Figure 1 presents derivative, , responses for a number of situations; the intent is to demonstrate the structure of the solutions and the influence of properties, specifically the exponent and the ratio . As one of our interests here is to demonstrate the accuracy of the solutions we assume the diffusivity ratio, , to be 1. All properties considered are noted in Figure 1. The symbol that is the dimensionless distance to the fault is defined by
Fig. 1 Responses at a well producing a composite system under subdiffusive flow. All solutions generated through equation (36) for of illustrate derivative responses, . The influences of , and the exponent, , are considered; the properties of Zone are held constant. The two dashed lines marked AN reflect calculations obtained through equation (59) and are in excellent agreement with the rigorous solutions which meet the two stipulations (equal ’s and ’s) stated. The upward rise in the responses as the magnitude of the exponent decreases reflects the lower effective permeability of Zone 2. 
Simulations obtained through equation (36) for of 0 are shown as an unbroken lines. Because the properties of Zone 1 are held constant, all solutions originate from a single curve and branch out at later times as the result of the influence of Zone 2. Considering the specifics of the results shown here the bottommost curve labelled A and topmost curve labelled D correspond to values of T_{D} of 100 and 0.01, respectively, for α_{j} = 0.9; j = 1, 2 with L_{D} of 100; the longtime approximations shown as dashed lines (marked AN) for identical conditions are in excellent agreement as they meet the stipulations in equation (59). Curves A and D may also be considered proxies for constant pressure and closed boundaries, respectively. Curves B and C correspond to α_{2} of 0.7 and 0.5, respectively, each for T_{D} of 100; all three curves, A, B and C deviate from each other as the influence of Zone 2 becomes dominant. Not unexpectedly, longtime responses reflect a system of lower effective permeability as α_{2} is smaller. It is interesting that Curves C and D essentially merge at long enough times. This result, as usual, suggests the difficulties one encounters in addressing the uniqueness of results obtained by evaluating the pressuretime curve; not unexpectedly supplementary information helpful (see Sect. 4). In Figure 2 we continue the exploration of the influence of properties, specifically the influences of and . The solutions are for of and of . The top set corresponds to dimensionless pressure, , response, and the bottom set represents the corresponding derivative, . Again the properties of Zone are held constant; thus all solutions at later times emanate from a single curve, the branching depends on and/or as soon as the region beyond the interface influences the well response. As in Figure 1, we assume to be 100. The dashed line (AN) reflects derivative responses given in equation (59) when both the diffusivities, η’s and α’s are equal. Two sets of solutions, Set A and Set B, are considered for α_{2} of 0.9 and 0.5 and three values of η_{D}: 5 (squares), 1 (unbroken lines) and 0.2 (squares). Considering the solutions for α_{2} of 0.5, namely Set A, we see that the influence of is small – simulations for η_{D} = 0.2, the circles, and that for η_{D} = 5, the squares, straddle the η_{D} = 1 solution. If the α_{j}’s are equal (0.9 case), however, the influence of η_{D} is insignificant, Set B. These observations carry over to curves. In essence, these results, and others not discussed here, show that the Bixel et al. (1963) observation that the influence of η_{D} is insignificant carries over to the subdiffusive solutions but the α_{j}’s must be equal.
Fig. 2 The well response for production at a constant rate for a composite system; pressure, , and the corresponding derivative, , responses, the top and bottom sets of curves, respectively, are shown. The intent is to demonstrate the influence of diffusivity ratio, , and the exponent, . The unbroken lines presume that the ’s are equal; the symbols denoted as circles and squares presume that . If the ’s are equal, then the influence of η_{D} is negligibly small; see Curves labelled B; the corresponding derivative curve shows that the analytical solution shown as a dashed line (marked AN) predicts values in excellent agreement. The circles and squares illustrate the influence of for . In general, solutions depend on both and . 
2.2 Influence of a contact resistance
Figure 3 presents the influence of a fault, , located at a distance of of at an observation point, with both the well and the observation location on the same side of the interface. The influences of T_{D}, η_{D} and α_{2} for fixed values of Zone 1 properties with α_{1} of , which result in all curves starting from the same point, are presented. All solutions except the line marked as Label A are for of ; the one identified as A is for an . Dimensionless times in Figure 3 are normalized with respect ; that is, all solutions shown here are expressed as . The unbroken lines identical to the ones given in Raghavan and Ozkan (2011) correspond to the classical case with the variable of interest being ; we consider three values of : , and . The peak in the derivative curve, not unexpectedly, reflects the restriction to flow as a result of the fault, with the slight concave upward behavior being the influence of the ratio . Ultimately all lines become flat and reflect as both ’s are 1 (see Eq. (59)).
Fig. 3 Influence of contact resistance on derivative responses at an observation well, . Both wells are on the same side of the fault, hence the sharp rise in the derivative curve. The principal goal is to demonstrate the influence of , and . Properties of Zone are held constant and classical diffusion prevails (); unbroken lines are for α_{2} = 1; dashed lines for and small dashes for . 
The situation is different if ; the dashed lines in Figure 3 reflect solutions for of for identical values of . These solutions deviate from the unbroken lines with the upward trend representing the fact of an additional restriction, namely that is now . A further reduction in will result in an earlier deviation and the magnitude of the upward trend would be larger; see line with small dashes for of which represents of . It is clear that subdiffusive characteristics of Zone 2 are reflected across the fault.
Earlier the influence of the porositycompressibility ratio was addressed and it was noted that the concave upward portion was reflective of the value of . The basis of this observation may be understood by comparing the curve labelled A which corresponds to and with its counterpart the solution for and .
Results that are the complement of the ones in Figure 3 are discussed in Figure 4. Incidently, here, we also draw attention to particular features of subdiffusion. The unbroken lines are responses for and ; that is, subdiffusive flow in Zone and classical flow in Zone 2 – the opposite of situations considered in Figure 3 although a slightly different value of is used and the range of values is larger. While comparing the responses in Figures 3 and 4 the significant change in the values of the ordinate needs to be borne in mind because a change in the magnitude of from is results in a significant change in the magnitude of . All solutions considered here except for the bottommost curve assume subdiffusive flow in Zone , and, again all solutions for of start out together with the steep rise reflecting the conductivity of the fault; it is clear that subdiffusive effects of Zone are evident even when the properties of Zone 2 are dominant unlike what we would expect if were 1.0 as shown by the bottommost curve in Figure 3. The two dashed lines presume subdiffusive flow in Zone ; initially responses fall below that for of reflecting the speed of propagation as discussed in Section 1. With time, however, the slopes of the curves become steeper reflecting a lower effective permeability as decreases. In practice such steeper slopes would not reflect additional boundary effects.
Fig. 4 Influence of contact resistance on derivative responses at an observation well, . Again, both wells are on the same side of the fault, the change in from 1.0 to 0.9 produces a significant change in the magnitude of the response. 
3 Discussion
Before closing with a few general observations, for practicality, we note that the exponent may be rendered as a Probability Density Function (PDF) through the instantaneous Green’s function, G(r, t), through the expression in equation (18). The first term yields:
(61)where is the distance from the source. By considering a series of transformations of the Fox function Schneider and Wyss (1989) arrive at an expression for in the following form,
(62)where the symbol is the number of dimensions; that is, equation (62) addresses a far more general problem than the one we consider for in our case. exhibits powerlaw behavior at longenough times for . This expression agrees with those in Carslaw and Jaeger (1959) for if . Grebenkov (2010) arrives at the same expression as that given in equation (62) but for by using the result
(63)with . Incidentally, the lefthand side of equation (63) may also be expressed in terms of the Wright function, , as indicated in Povstenko (2015) by
The moments of an even order which would also be of interest are (Uchaikin, 2013); see equation (1):
Variations of as a function of r are shown in Figure 5 for three values of : , , and 0.8. These values were computed independent of the expressions (62) and (63). The characteristic bell shape is evident only for . Note that may be estimated by equation (62) for for of .
4 Commentary
Although our goal is to place a focus on the mathematical aspects and characteristics of the responses under subdiffusive flow, we briefly mention a few practical aspects we follow. Evaluation of well responses is to be done in the same manner as that for classical diffusion. The additional complication is the determination of the exponents (or at least ) which would be the first step in identifying subdiffusive flow; the determination of the exponent would depend, among other things, on the duration of the test; for in the case of other factors such as the properties of the rock or well locations we may not have much control. The exponent is best determined through the derivative curve, , as this curve when combined with the pressure response, , better constrains the “uncertainty in the solution space”; see Raghavan (2004). It is for this reason we have, as mentioned, discussed the derivative extensively in the text. Given that the process is to be data driven we should garner as much information as possible regarding reservoir architecture (maps, cores, logs, etc.); see Raghavan et al. (2001) for a comprehensive review of combining, geology, geophysics and production statistics. In addition, it is essential to know completion information; see Raghavan (1993). We should note that although the analysis is conducted on the Theis (1935) scale, the use of a subdiffusive model does not imply a single value for properties. The end product is to arrive at a reservoir description that reproduces future performance. For all of this the derivative curve should be the starting point.
5 Concluding remarks
A 2D source function solution under subdiffusion in two contiguous media separated by a partially absorbing or communicating interface is documented. In all other aspects the media are infinite in extent. Two dimensional subdiffusion problems, to our knowledge, are yet to be considered. The role of the interface exerting a resistance is detailed. For the situation described exact and asymptotic solutions are derived. Computational results are discussed in the context of applications to the earth sciences; specifically the solutions may be applied to hydrocarbon or geothermalfluid extraction, flow of groundwater, roles of faults, and sequestration. Extensions to other situations related to flow in porous media are readily possible: variablerate problems through Duhamel’s theorem including situations involving wellbore storage effects, multiple, contiguous media, complex wellbores such as inclined, fractured and horizontal wells and equivalent porous media models such as those described in Barenblatt et al. (1960) and de SwaanO (1976). Quantitative descriptions of subdiffusive flow and their effects on estimating properties of the rock and permeable interfaces are discussed. Should subdiffusive effects exist it is not possible to ignore the existence of their timedependent influences at long enough times as in the case of classical diffusion.
In summary, the Green’s function formulation to address the role of a partially communicating interface detailed here is a general one; we may tap into the potential of the solutions discussed here in a number of contexts such as biophysics, fractal media, and the like, as Laplace transformations and the method of sources and sinks are used. As mentioned in the Introduction, the expressions derived here are directly transferable to other fields; only a change in Nomenclature is needed.
References
 Albinali A., Holy R., Sarak H., Ozkan E. (2016) Modeling of 1D anomalous diffusion in fractured nanoporous media, Oil Gas Sci. Technol. Rev. IFP Energies nouvelles 71, 56. [CrossRef] [Google Scholar]
 Angulo J.M., RuizMedina M.D., Anh V.V., Grecksch W. (2000) Fractional diffusion and fractional heat equation, Adv. Appl. Probab. 32, 4, 1077–1099. [Google Scholar]
 Barenblatt G.I., Zheltov YuP, Kochina I.N. (1960) Basic concepts in the theory of seepage of homogeneous liquids in fissured rocks, J. Appl. Math. Mech. 24, 5, 1286–1303. [CrossRef] [Google Scholar]
 Bense V., Person M. (2006) Faults as conduitbarrier system to fluid flow in silicilastic sedimentary aquifers, Water Resour. Res. 42, W05421. doi: 10.1029/2005WR004480. [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, 647–658. doi: 10.1016/j.jhydrol.2006.01.008. [CrossRef] [Google Scholar]
 Bixel H.C., Larkin B.K., Van Poollen H.K. (1963) Effect of linear discontinuities on pressure buildup and drawdown behavior, J. Pet. Tech. 15, 8, 885–895. doi: 10.2118/611PA. [CrossRef] [Google Scholar]
 Caine J.S., Evans J.P., Forster C.B. (1996) Fault zone architecture and permeability structure, Geology 24, 11, 1025–1028. [Google Scholar]
 Caputo M. (1967) Linear models of dissipation whose Q is almost Frequency IndependentII, Geophys. J. R. Astron. Soc. 13, 5, 529–539. [NASA ADS] [CrossRef] [Google Scholar]
 Caputo M. (1999) Diffusion of fluids in porous media with memory, Geothermics 28, 1, 113–130. [Google Scholar]
 Carslaw H.S., Jaeger J.C. (1959) Conduction of heat in solids, 2nd edn., Vol. 22, Clarendon Press, Oxford, pp. 319–326, 327–352. [Google Scholar]
 Chang J., Yortsos Y.C. (1990) Pressuretransient analysis of fractal reservoirs, SPE Form. Eval. 5, 1, 31–39. [CrossRef] [Google Scholar]
 Childs E.C. (1969) An Introduction to the Physical Basis of Soil Water Phenomena, John Wiley and Sons Ltd, London, pp. 153–178. [Google Scholar]
 CincoLey H., SamaniegoV. F., Kuchuk F. (1985) The pressure transient behavior for naturally fractured reservoirs with multiple block size. In: Paper SPE 14168 Presented at SPE Annual Technical Conference and Exhibition, 22–26 September, Las Vegas, Nevada. [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. [CrossRef] [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]
 Darcy H. (1856) Les Fontaines publiques de la ville de Dijon, Exposition et application des principes à suivre et des formules à employer dans les questions de distribution d’eau ; ouvrage terminé par un appendice relatif aux fournitures d’eau de plusieurs villes au filtrage des eaux et à la fabrication des tuyaux de fonte, de plomb, Victor Dalmont, Paris, pp. 590–594. [Google Scholar]
 Dassas Y., Duby Y. (1995) Diffusion toward fractal interfaces, potentiostatic, galvanostatic, and linear sweep voltammetric techniques, J. Electrochem. Soc. 142, 12, 4175–4180. [Google Scholar]
 de SwaanO A. (1976) Analytical solutions for determining naturally fractured reservoir properties by Well testing, Soc. Pet. Eng. J. 16, 3, 117–122. doi: 10.2118/5346PA. [CrossRef] [Google Scholar]
 Erdélyi A., Magnus W., Oberhettinger F., Tricomi F.G. (1954) Tables of integral transforms, based, in part, on notes left by Harry Bateman and Compiled by the Staff of the Bateman Manuscript Project, vol. 1, McGrawHill, New York, 17. [Google Scholar]
 Erdelyi A., Magnus W.F., Oberhettinger F., Tricomi F.G. (1955) Higher transcendental functions, vol. 3, McGrawHill, New York, pp. 206–227. [Google Scholar]
 Evans J.P. (1988) Deformation mechanisms in granitic rocks at shallow crustal levels, J. Struct. Geol. 10, 5, 437–443. [Google Scholar]
 Feller W. (1971) An introduction to probability theory and its applications. II, 2nd edn., Wiley, New York, pp. 8–10, 50. [Google Scholar]
 Fomin S., Chugunov V., Hashida T. (2011) Mathematical modeling of anomalous diffusion in porous media, Fractional Differ. Calc. 1, 1–28. [Google Scholar]
 Fu L., Milliken K.L., Sharp J.M. Jr. (1994) Porosity and permeability variations in fractured and liesegangbanded Breathitt sandstones (Middle Pennsylvanian), eastern Kentucky: Diagenetic controls and implications for modeling dualporosity systems, J. Hydrol. 154, 1–4, 351–381. [CrossRef] [Google Scholar]
 Gefen Y., Aharony A., Alexander S. (1983) Anomalous diffusion on percolating clusters, Phys. Rev. Lett. 50, 1, 77–80. doi: 10.1103/PhysRevLett.50.77. [Google Scholar]
 Gradshteyn I.S., Ryzhik I.M. (1994) Table of integrals, series and products, 5th edn., Vol. 532, Academic Press Inc., Orlando, pp. 3.961–3.962. [Google Scholar]
 Grebenkov D.B. (2010) Subdiffusion in a bounded domain with a partially absorbingreflecting boundary, Phys. Rev. E: Stat. Phys. Plasmas Fluids 81, 1, 021128. [CrossRef] [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. doi: 10.1007/BF00281373. [CrossRef] [MathSciNet] [Google Scholar]
 Hänninen J.J., Pirjola R.J., Lindell I.V. (2002) Application of the exact image theory to studies of ground effects of space weather, Geophys. J. Int. 151, 2, 534–542. [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]
 Hilfer R., Anton L. (1995) Fractional master equations and fractal time random walks, Phys. Rev. E 51, 2, R848–R851. [Google Scholar]
 Jourde H., Pistrea S., Perrochet P., Droguea C. (2002) Origin of fractional flow dimension to a partially penetrating well in stratified fractured reservoirs, new results based on the study of synthetic fracture networks, Adv. Water Res. 25, 4, 371–387. [CrossRef] [Google Scholar]
 Kenkre V.M., Montroll E.W., Shlesinger M.F. (1973) Generalized master equations for continuoustime random walks, J Stat Phys 9, 1, 45–50. [Google Scholar]
 Kim S., Kavvas M.L., Ercan A. (2015) Fractional ensemble average governing equations of transport by timespace nonstationary stochastic fractional advective velocity and fractional dispersion, II: Numerical investigation, J. Hydrol. Eng. 20, 2, 04014040. [Google Scholar]
 Kosztoowicz K. (2017) Subdiffusion in a system consisting of two different media separated by a thin membrane, Int. J. Heat Mass Transfer 111, 1322–1333. [CrossRef] [Google Scholar]
 Lẽ Mehautẽ A., Crepy G. (1983) Introduction to transfer and motion in fractal media: The geometry of kinetics, Solid State Ion. 1, 9–10, 17–30. [Google Scholar]
 Lindell I., Alanen E. (1984) Exact image theory for the Sommerfeld halfspace problem, part III: General formulation, IEEE Trans. Antennas Propag. 32, 10, 1027–1032. doi: 10.1109/TAP.1984.1143204. [Google Scholar]
 Lindell I., Alanen E., von Bagh H. (1986) Exact image theory for the calculation of fields transmitted through a planar interface of two media, IEEE Trans. Antennas Propag. 34, 2, 29–137. doi: 10.1109/TAP.1986.1143788. [Google Scholar]
 Lindquist E. (1933) On the flow of water through porous soil, 1st Congres des Grands Barrerges, Stockholm, 5, 91–101. [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. doi: 10.1016/j.micromeso.2013.02.054. [CrossRef] [PubMed] [Google Scholar]
 Mainardi F. (2010) Fractional calculus and waves in linear viscoelasticity, Imperial College Press, London, pp. 211–236. [Google Scholar]
 Mainardi M., Pagnini G., Saxena R.K. (2005) Fox H functions in fractional diffusion, J. Comput. Appl. Math. 178, 1–2, 321–331. [Google Scholar]
 Mandelis A., Nicolaides L., Chen Y. (2001) Structure and the reflectionless/refractionless nature of parabolic diffusionwave fields, Phys. Rev. Lett. 87, 2, 020801. doi: 10.1103/PhysRevLett.87.020801. [Google Scholar]
 Mathai A.M., Saxena R.K. (1978) The Hfunction with applications in statistics and other disciplines, Wiley, New Delhi, India, pp. 1–19. [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]
 Miller B. (2005) The Baton Rouge Fault: Conduit or impediment to groundwater flow? in: Paper Presented at 54th Annual Meeting Southeast, Sect. Geol. Soc. Am., Biloxi Miss. [Google Scholar]
 Miller F.G. (1954) Multiphase flow theory and the problem of spacing oil wells, United States Bureau of Mines 529, 8–10. [Google Scholar]
 Mitchell T.M., Faulkner D.R. (2009) The nature and origin of offfault damage surrounding strikeslip fault zones with a wide range of displacements: A field study from the Atacama fault system, northern Chile, J. Struct. Geol. 31, 8, 802–816. doi: 10.1016/j.jsg.2009.05.002. [Google Scholar]
 Molz F.J. III, Fix G.J. III, Lu S.S. (2002) A physical interpretation for the fractional derivative in Levy diffusion, Appl. Math. Lett. 15, 7, 907–911. [Google Scholar]
 Montroll E.W., Weiss G.H. (1965) Random walks on lattices II, J. Math. Phys. 6, 167–181. [Google Scholar]
 Moodie T.B., Tait R. (1983) On thermal transients with finite wave speeds, J. Acta Mech. 50, 1–2, 97–104. doi: 10.1007/BF01170443. [CrossRef] [Google Scholar]
 Neville I.R., Sharp J.M. Jr., Kreisel I. (1998) Contaminant transport in sets of parallel finite fractures with fracture skins, J. Contam. Hydrol. 31, 1–2, 83–109. [Google Scholar]
 Nigmatullin R. (1984a) On the theory of relaxation with “remnant” memory, Phys. Stat. Sol. B 124, 1, 389–393. doi: 10.1002/pssb.2221240142. [CrossRef] [Google Scholar]
 Nigmatullin R.R. (1984b) To the theoretical explanation of the universal response, Phys. Status Solidi B Basic Res. 123, 2, 739–745. [CrossRef] [Google Scholar]
 Noetinger B., Estebenet T. (2000) Upscaling of double porosity fractured media using continuoustime random walks methods, Transp. Porous Med. 39, 3, 315–337. [CrossRef] [Google Scholar]
 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. doi: 10.1115/1.3422771. [CrossRef] [Google Scholar]
 Oberhettinger F., Badii L. (1973) Tables of Laplace transforms, Springer Verlag, Berlin, p. 268. [Google Scholar]
 Odling N.E., Gillespie P., Bourgine B., Castaing C., Chiles J.P., Christiansen N.P., Fillion E., Albert G., Olsen C., Lena T., Trice R., Aarseth E.S., Walsh J.J., Watterson J. (1999) Variations in fracture system geometry and their implications for fluid flow in fractured hydrocarbon reservoirs, Petrol. Geosci. 5, 373–384. [CrossRef] [Google Scholar]
 O’Shaughnessy B., Procaccia I. (1985a) Analytical solutions for diffusion on fractal objects, Phys. Rev. Lett. 54, 5, 455–458. [Google Scholar]
 O’Shaughnessy B., Procaccia I. (1985b) Diffusion on fractals, Phys. Rev. A: At. Mol. Opt. Phys 32, 5, 3073–3083. [CrossRef] [MathSciNet] [PubMed] [Google Scholar]
 Philip J.R. (1957) Transient fluid motions in saturated porous media, Aust. J. Phys. 10, 43–53. [CrossRef] [Google Scholar]
 Płociniczak Ł. (2015) Analytical studies of a timefractional porous medium equation. Derivation, approximation and applications, Commun. Nonlinear Sci. Numer. Simul. 24, 1–3, 169–183. [Google Scholar]
 Povstenko Y. (2015) Linear fractional diffusionwave equation for scientists and engineers, Birkhäuser 24–30, 32, 34. [Google Scholar]
 Povstenko Y.Z. (2013) Fractional heat conduction in infinite onedimensional composite medium, J. Therm. Stresses 36, 4, 351–363. doi: 10.1080/01495739.2013.770693. [CrossRef] [Google Scholar]
 Povstenko Y., Kyrylych T. (2019) Fractional heat conduction in solids connected by thin intermediate layer: Nonperfect thermal contact, Continuum Mech. Thermodyn. 31, 1719–1731, doi: 10.1007/s0016101900750w. [CrossRef] [Google Scholar]
 Prats M., Raghavan R. (2013) Finite horizontal well in a uniformthickness reservoir crossing a natural fracture, SPE J. 18, 5, 982–992. doi: 10.2118/163098PA. [CrossRef] [Google Scholar]
 Prats M., Raghavan R. (2014) Properties of a natural fracture and its skins from reservoir well tests, SPE J. 19, 3, 390–397. doi: 10.2118/167262PA. [CrossRef] [Google Scholar]
 Raghavan R. (1993) Well test analysis, Prentice Hall, Englewoods Cliffs, NJ, pp. 6–8, 13. [Google Scholar]
 Raghavan R. (2004) A review of applications to constrain pumptest responses to improve on geological description and uncertainty, Rev. Geophys. 42, 1–29. [Google Scholar]
 Raghavan R. (2010) A composite system with a planar interface, J. Pet. Sci. Eng. 70, 3–4, 229–234. doi: 10.1016/j.petrol.2009.11.015. [Google Scholar]
 Raghavan R. (2011) Fractional derivatives: Application to transient flow, J. Pet. Sci. Eng. 80, 1, 7–13. doi: 10.1016/j.petrol.2011.10.003. [Google Scholar]
 Raghavan R. (2012) A horizontal well in a composite system with planar interfaces, Adv. Water Res. 38, 38–47. doi: 10.1016/j.advwatres.2011.12.009. [CrossRef] [Google Scholar]
 Raghavan R., Chen C. (2017) Addressing the influence of a heterogeneous matrix on well performance in fractured rocks, Transp. Porous Med. 117, 1, 69–102. doi: 10.1007/s1124201708205. [CrossRef] [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., Ozkan E. (2011) Flow in composite slabs. SPE Journal 16, 2, 374–387. doi: 10.2118/140748PA. [CrossRef] [Google Scholar]
 Raghavan R., Chen C. (2018) Time and space fractional diffusion in finite systems, Transp. Porous Med. 123, 173–193. [CrossRef] [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. doi: 10.2516/ogst/2018081. [CrossRef] [Google Scholar]
 Raghavan R., Dixon T.N., Phan V.Q., Robinson S.W. (2001) Integration of geology, geophysics, and numerical simulation in the interpretation of a well test in a fluvial reservoir, SPE Reserv. Evalu. Eng. 4, 3, 201–208. doi: 10.2118/72097PA. [CrossRef] [Google Scholar]
 Savage H.M., Brodsky E.E. (2011) Collateral damage: Evolution with displacement of fracture distribution and secondary fault strands in fault damage zones, J. Geophys. Res. 116, B03405. doi: 10.1029/2010JB007665. [Google Scholar]
 Saxena R.K., Mathai A.M., Haubold H.J. (2006) Fractional reactiondiffusion equations, Astrophys. Space Sci. 305, 3, 289–296. [Google Scholar]
 Schneider W.R., Wyss W. (1989) Fractional diffusion and wave equations, J. Math. Phys. 30, 134. doi: 10.1063/1.528578. [Google Scholar]
 Scholz C.H., Dawers N.H., Yu J.Z., Anders M.H., Cowie P.A. (1993) Fault growth and fault scaling laws: Preliminary results, J. Geophys. Res. 98, B12, 21951–21961. [Google Scholar]
 Sharp J.M. Jr., Kreisel I., Milliken K.L., Mace R.E., Robinson N.I. (1996) Fracture skin properties and effects on solute transport: Geotechnical and environmental implications, in: M. Aubertin, F. Hassam, H. Mitri (eds), Rock Mechanics, Tools and Techniques, Balkema, Rotterdam, pp. 1329–1335. [Google Scholar]
 Shendeleva M.L. (2004) Instantaneous line heat source near a plane interface, J. Appl. Phys. 95, 5, 2839–2845. [Google Scholar]
 Sommerfeld A. (1909) Uber die Ausbreitung der Wellen in der drahtlosen Telegraphie, Ann. Phys. 28, 665–736. [Google Scholar]
 Stehfest H. (1970a) Algorithm 368: Numerical inversion of Laplace transforms [D5], Commun. ACM 13, 1, 47–49. [Google Scholar]
 Stehfest H. (1970b) Remark on algorithm 368: Numerical inversion of Laplace transforms, Commun. ACM 13, 10, 624. [Google Scholar]
 Su N., 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. doi: 10.1016/j.jhydrol.2015.09.033. [CrossRef] [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. doi: 10.1016/j.jhydrol.2014.09.021 [CrossRef] [Google Scholar]
 Suzuki A.T., Hashida K.L., Horne R.N. (2016) Experimental tests of truncated diffusion in fault damage zones, Water Resour. Res. 52, 8578–8589. doi: 10.1002/2016WR019017. [Google Scholar]
 Tao S., Gao X., Li C., Zeng J., Zhang X., Yang C., Zhang J., Gong Y. (2016) The experimental modeling of gas percolation mechanisms in a coalmeasure tight sandstone reservoir: A case study on the coalmeasure tight sandstone gas in the Upper Triassic XuJiahe Formation, Sichuan Basin, China, J. Nat. Gas Geosci. 1, 6, 445–455. doi: 10.1016/j.jnggs.2016.11.009. [CrossRef] [Google Scholar]
 Theis C.V. (1935) The relationship between the lowering of the piezometric surface and the rate and duration of discharge of a well using groundwater storage, Eos Trans. AGU 2, 519–524. [CrossRef] [Google Scholar]
 Tualle J.M., Tinet E., Prat J., Avrillier S. (2000) Light propagation near turbidturbid planar interfaces, Opt. Commun. 183, 5–6, 337–346. [Google Scholar]
 Uchaikin V.V. (2013) Fractional derivatives for physicists and engineers. Volume I: Background and Theory, Springer, New York, 151, 296. [Google Scholar]
 Yanga S., Zhoua H.W., Zhang S.Q., Ren W.G. (2019) A fractional derivative perspective on transient pulse test for determining the permeability of rocks, Int. J. Rock Mech. Min. Sci. 113, 92–98. [CrossRef] [Google Scholar]
 Yaxley L.M. (1987) Effect of a partially communicating fault on transient pressure behavior, SPE Form. Eval. 2, 4, 590–598. doi: 10.2118/14311PA. [CrossRef] [Google Scholar]
 Zhokh A., Strizhak P. (2018) NonFickian transport in porous media: Always temporally anomalous? Transp. Porous Med. 124, 2, 309–323. doi: 10.1007/s1124201810666. [CrossRef] [Google Scholar]
 Zhokh A., Strizhak P. (2019) Investigation of the anomalous diffusion in the porous media: A spatiotemporal scaling, Heat Mass Transfer. 55, 1–10. doi: 10.1007/s00231019026024. [CrossRef] [Google Scholar]
All Figures
Fig. 1 Responses at a well producing a composite system under subdiffusive flow. All solutions generated through equation (36) for of illustrate derivative responses, . The influences of , and the exponent, , are considered; the properties of Zone are held constant. The two dashed lines marked AN reflect calculations obtained through equation (59) and are in excellent agreement with the rigorous solutions which meet the two stipulations (equal ’s and ’s) stated. The upward rise in the responses as the magnitude of the exponent decreases reflects the lower effective permeability of Zone 2. 

In the text 
Fig. 2 The well response for production at a constant rate for a composite system; pressure, , and the corresponding derivative, , responses, the top and bottom sets of curves, respectively, are shown. The intent is to demonstrate the influence of diffusivity ratio, , and the exponent, . The unbroken lines presume that the ’s are equal; the symbols denoted as circles and squares presume that . If the ’s are equal, then the influence of η_{D} is negligibly small; see Curves labelled B; the corresponding derivative curve shows that the analytical solution shown as a dashed line (marked AN) predicts values in excellent agreement. The circles and squares illustrate the influence of for . In general, solutions depend on both and . 

In the text 
Fig. 3 Influence of contact resistance on derivative responses at an observation well, . Both wells are on the same side of the fault, hence the sharp rise in the derivative curve. The principal goal is to demonstrate the influence of , and . Properties of Zone are held constant and classical diffusion prevails (); unbroken lines are for α_{2} = 1; dashed lines for and small dashes for . 

In the text 
Fig. 4 Influence of contact resistance on derivative responses at an observation well, . Again, both wells are on the same side of the fault, the change in from 1.0 to 0.9 produces a significant change in the magnitude of the response. 

In the text 
Fig. 5 The probability density function as a function of α in Figure 5; η = 1 and t = 1. 

In the text 