Subdiffusive flow in a composite medium with a communicating (absorbing) interface

Two-dimensional 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.).


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 line-source 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 with r being the distance, t being time and where a 6 ¼ 1 is a constant. Equation (1) for a < 1, 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 fracture-matrix interfaces, are modeled in classical situations as skin regions Raghavan, 2013, 2014) (ignoring scaling laws), and the phenomenon is known as interporosity skin (Cinco-Ley 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 well-attested 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 steady-state 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, cðx; tÞ, 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): where the symbol S represents the source, the symbol S 0 represents the image sources that are distributed to reflect the interfaces, the symbol P 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 symbolr epresents convolution in both space and time.

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, X, in R 2 on either side of the interface is examined. A contact resistance at the interface governs fluid movement and the system is produced through a line-source well located at the point (0, y 0 ). 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 l. Initially, throughout the system the pressure, pðx; tÞ, is assumed to be uniform at the value of p i . The transient diffusion equation for this configuration is where the symbols vðx; tÞ and /ðxÞ 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, and also pðfx; yg ! 1; tÞ ¼ p i : Assuming flow to a line-source well, the condition at a wellbore for a specified withdrawal rate, q, is where B is the formation volume factor, and r is the distance between the point (x, y) and the well location (0, y 0 ). 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 v x 1 ðx; y; tÞ ¼ a½p 1 ðx; y; tÞ À p 2 ðx; y; tÞ; y ¼ 0; and v x 1 ðx; y; tÞ ¼ v x 2 ðx; y; tÞ; y ¼ 0; where the symbol a ¼ k F =ðlw F Þ 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, steady-state 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.

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 cross-sectional 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, Dx, over a fixed time, Dt.
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 power-law 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, g, follows the relation g $ r Àðdw À2Þ , 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 Continuous-Time 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, Dx, at intervals that are separated through a fixed time interval, Dt, the CTRW formulation uses an i.i.d. (independent, identically distributed) random variable for increments of the step length, Dx, and a waiting (pausing) time function, w(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 wðtÞ ¼ ðt is the Mittag-Leffler function discussed in many texts (Erdelyi et al., 1955;Mainardi, 2010;Povstenko, 2015;Uchaikin, 2013): where C represents the complex plane and where CðÁÞ is the Gamma function. The symbol x is a constant that depends on the properties of the fractal object; with 0 x 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, Dx, that is an even function and the other a power-law waiting time function, w(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, wðtÞ.
The CTRW formulation leads to a flux law of the following form: Here a is the subdiffusion coefficient with 0 < a < 1, l is the fluid viscosity,k is the permeability corresponding to a and o a f ðtÞ=ot a , is the fractional derivative defined in Caputo (1967) For a = 1, equation (10) reverts to the classical Darcy equation. The Laplace transformation of the Caputo fractional operator, c 0 D a t f ðtÞ, 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).
While on 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).

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, Gðr; tÞ, in R 2 . Accordingly, for Zone 1 and Zone 2 we write instantaneous solutions, c j (x, y); j = 1, 2, respectively, in terms of their Laplace transformation, c j ðx; yÞ with the source located at (0; y 0 ) to be and c 2 ðx; yÞ ¼ whereg j , the "diffusivity" of the medium, is given bỹ and 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 K 0 ðxÞ 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 (u) and f 2 (u), we obtain the following expressions: or and, where The functions A 1 , A 2 and A 3 are given, respectively, by and Here and T j ¼k j h 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 Dp (x, y, z; t) (Raghavan and Ozkan, 1994) is Here (x 0 , y 0 ) is the location of the source, (x, y) any location in the porous rock and the strength of the continuous linesource isqðx 0 ; y 0 ; tÞ=ð/cÞ with dS 0 being the element (a line or surface) through which fluid is extracted. For an instantaneous line-source with the source strength being constant, we have Thus on using equation (34) and the expressions for c j ðx; yÞ derived above, the expressions for Áp j ðx; y; tÞ in terms of their respective Laplace transformations, Áp j ðx; yÞ, are or in Zone 1 (y > 0), and in Zone 2 (y < 0).

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: and where ' is the reference length.
Using these definitions, we obtain the following expression for the pressure distributions, p D j ðx D Þ; j = 1 and 2, in terms of the Laplace transformation for Zones 1 and 2 to be, respectively, for y D > 0 and for y D < 0. 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 !,g r and T r where But the A j 's; that is, both q j 's, and the three variables, T r , g r , and Y r , must, of course, be rescaled. The scalings are: with and Scalings involving the ratio of a j 's enter into the picture because we consider subdiffusive flow; in accordance with equation (10); the dimension ofk is a function of a. These scalings also follow if the partial differential equations are considered. The groupings T D and g D may be replaced by T D and ð/c t hÞ D if desired, where ð/c t hÞ D ¼ ð/c t hÞ 1 ð/c t hÞ 2 : ð50Þ A number of approximate solutions may be derived if q 1 = q 2 , which impliesg 1 andg 2 and also a 1 = a 2 . In the following we use ' ¼ r w , where r w is the well radius. In this situation for Zone 1 and Zone 2, respectively, we have and The right-hand sides of equations (51) and (52) where the symbol H m;n p;q x ða 1 ;A 1 Þ;:::;ða p ;A p Þ ðb 1 ;B 1 Þ;:::;ðb q ;B q Þ " # is the Fox function or the H-function, Rð.Þ > 0, Rðz 2 Þ > 0, RðsÞ > 0 as noted in Saxena et al. (2006), and thus we have the following result for p D 1 ðx D ; t D Þ, namely and for p D 2 ðx D ; t D Þ, 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 long-time approximation for this specific situation by considering the following expression for K 0 ðxÞ as x ? 0 (implicitly we consider situations where s ? 0), namely using where c is Euler's constant. The long-time pressure response, p D (x D ; t D ?1), may be obtained by using the first two terms of the above expression, as well as the result of Oberhettinger and Badii (1973) s Àm ln s ¼ ½CðmÞ À1 t mÀ1 ½wðmÞ À ln t; where w() is the Digamma function, and we therefore have Also, the logarithmic derivative, p 0 D ðx D ; t D Þ, 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 a = 1 in an infinitely large system (T D = 1), respectively, correspond to the Theis (1935) solution and the semi-logarithmic approximation of the Theis (1935) solution, namely the Cooper and Jacob (1946) approximation. Along the same lines if both a 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 forg 2 different fromg 1 at this time. Although this is the case, we know that for classical diffusion (a j ¼ 1) the influence ofg r is small and equations (58) and (59) hold in most cases when the a j 's are 1; see Bixel et al. (1963). That being so, it is plausible that equations (58) and (59) may hold forg r 6 ¼ 1 provided that the a j 's were equal. This, of course, would have to be explored.

Results
Our main concern in the following is the understanding of the characteristics of the derivative curve, p 0 D ðx D ; t D Þ, as these curves best demonstrate the characteristics of the rock under consideration. For each of the solutions presented, the pressure response, p D ðx D ; t D Þ, 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 well-test problems. Figure 1 presents derivative, p 0 wD ðt D Þ, responses for a number of situations; the intent is to demonstrate the structure of the solutions and the influence of properties, specifically the exponent a j and the ratio T D . As one of our interests here is to demonstrate the accuracy of the solutions we assume the diffusivity ratio, g D , to be 1. All properties considered are noted in Figure 1. The symbol L D that is the dimensionless distance to the fault is defined by

Perfect contact
Simulations obtained through equation (36) for Ç D 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 a j = 0.9; j = 1, 2 with L D of 100; the long-time 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 a 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, long-time responses reflect a system of lower effective permeability as a 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 pressure-time 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 g D and a j . The solutions are for T D of 100 and a 1 of 0:9. The top set corresponds to dimensionless pressure, p wD , response, and the bottom set represents the corresponding derivative, p 0 wD . Again the properties of Zone 1 are held constant; thus all solutions at later times emanate from a single curve, the branching depends on g D and/or a 2 as soon as the region beyond the interface influences the well response. As in Figure 1, we assume L D to be 100. The dashed line (AN) reflects derivative responses given in equation (59) when both the diffusivities, g's and a's are equal. Two sets of solutions, Set A and Set B, are considered for a 2 of 0.9 and 0.5 and three values of g D : 5 (squares), 1 (unbroken lines) and 0.2 (squares). Considering the p wD ðt D Þ solutions for a 2 of 0.5, namely Set A, we see that the influence of g D is smallsimulations for g D = 0.2, the circles, and that for g D = 5, the squares, straddle the g D = 1 solution. If the a j 's are equal (0.9 case), however, the influence of g D is insignificant, Set B. These observations carry over to p 0 wD ðt D Þ curves. In essence, these results, and others not discussed here, show that the Bixel et al. (1963) observation that the influence of g D is insignificant carries over to the subdiffusive solutions but the a j 's must be equal. Figure 3 presents the influence of a fault, Ç D ¼ 10 3 , located at a distance of L D of 100 at an observation point, x D ¼ 60; y D ¼ 40 with both the well and the observation location on the same side of the interface. The influences of T D , g D and a 2 for fixed values of Zone 1 properties with a 1 of 1, which result in all curves starting from the same point, are presented. All solutions except the line marked as Label A are for g D of 8; the one identified as A is for an g D ¼ 1. Dimensionless times in Figure 3 are normalized with respect L D ¼ 100; that is, all solutions shown here are expressed as t D =L 2 D . The unbroken lines identical to the ones given in Raghavan and Ozkan (2011) correspond to the classical case with the variable of interest being T D ; we consider three values of T D : 10, 1 and 0:5. 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 ð/c t Þ 1 =ð/c t Þ 2 . Ultimately all lines become flat and reflect 2=ð1 þ T D Þ as both a's are 1 (see Eq. (59)).

Influence of a contact resistance
The situation is different if a 2 < 1; the dashed lines in Figure 3 reflect solutions for a 2 of 0:9 for identical values of T D . These solutions deviate from the unbroken lines with the upward trend representing the fact of an additional restriction, namely that a 2 is now 0:9. A further reduction in a 2 will result in an earlier deviation and the magnitude of the upward trend would be larger; see line with small dashes for T D of 10 which represents a 2 of 0:8. It is clear that subdiffusive characteristics of Zone 2 are reflected across the fault.
Earlier the influence of the porosity-compressibility ratio was addressed and it was noted that the concave upward portion was reflective of the value of g D . The basis of this observation may be understood by comparing the curve labelled A which corresponds to g D ¼ 1 and T D ¼ 1 with its counterpart the solution for g D ¼ 8 and T D ¼ 1.
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 a 1 ¼ 0:9 and a 2 ¼ 1:0; that is, subdiffusive flow in Zone 1 and classical flow in Zone 2the opposite of situations considered in Figure 3 although a slightly different value of g D is used and the range of T D values is larger. While comparing the responses  (59) and are in excellent agreement with the rigorous solutions which meet the two stipulations (equalg j 's and a j 's) stated. The upward rise in the responses as the magnitude of the exponent a 2 decreases reflects the lower effective permeability of Zone 2. Fig. 2. The well response for production at a constant rate for a composite system; pressure, p wD , and the corresponding derivative, p 0 wD , responses, the top and bottom sets of curves, respectively, are shown. The intent is to demonstrate the influence of diffusivity ratio, g D , and the exponent, a 2 . The unbroken lines presume that the g D 's are equal; the symbols denoted as circles and squares presume that g D 6 ¼ 1. If the a j 's are equal, then the influence of g 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 g D for a 2 ¼ 0:5. In general, solutions depend on both g D and a D .
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 a 1 from 1:0 is 0:9 results in a significant change in the magnitude of p wD ðt D Þ. All solutions considered here except for the bottommost curve assume subdiffusive flow in Zone 1, and, again all solutions for a 2 of 0:9 start out together with the steep rise reflecting the conductivity of the fault; it is clear that subdiffusive effects of Zone 1 are evident even when the properties of Zone 2 are dominant unlike what we would expect if a 1 were 1.0 as shown by the bottommost curve in Figure 3. The two dashed lines presume subdiffusive flow in Zone 2; initially responses fall below that for a 2 of 1:0 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 a 2 decreases. In practice such steeper slopes would not reflect additional boundary effects.

Discussion
Before closing with a few general observations, for practicality, we note that the exponent a may be rendered as a Probability Density Function (PDF) through the instantaneous Green's function, G(r, t), through the expression in equation (18) where r 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 Gðr; tÞ in the following form, where the symbol d is the number of dimensions; that is, equation (62) addresses a far more general problem than the one we consider for d ¼ 2 in our case. Gðr; tÞ exhibits power-law behavior at long-enough times for a < 1. This expression agrees with those in Carslaw and Jaeger (1959) for d 3 if a ¼ 1. Grebenkov (2010) arrives at the same expression as that given in equation (62) with j ¼ ðd þ 1Þ=4 À 1. Incidentally, the left-hand side of equation (63) may also be expressed in terms of the Wright function, W ða; b; tÞ, as indicated in Povstenko (2015) by L À1 ½s Àb expðÀks a Þ ¼ t bÀ1 W ðÀa; b; Àkt Àa Þ; The moments of an even order which would also be of interest are (Uchaikin, 2013); see equation (1): n ¼ 0; 1; 2; :::: Variations of Gðr; tÞ as a function of r are shown in Figure 5 for three values of a: 1, 0:9, and 0:8. These values were Properties of Zone 1 are held constant and classical diffusion prevails (a 1 ¼ 1:0); unbroken lines are for a 2 = 1; dashed lines for a 2 ¼ 0:9 and small dashes for a 2 ¼ 0:8. computed independent of the expressions (62) and (63). The characteristic bell shape is evident only for a ¼ 1. Note that Gðr; tÞ may be estimated by equation (62) for r ¼ 0 for a of 1.

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 a j (or at least a 1 ) 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, Áp 0 , as this curve when combined with the pressure response, Áp, 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.

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 geothermal-fluid extraction, flow of groundwater, roles of faults, and sequestration. Extensions to other situations related to flow in porous media are readily possible: variable-rate 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 Swaan-O (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 time-dependent 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.