© R. Raghavan & C.-C. Chen, published by IFP Energies nouvelles, 2020

Nomenclature

Aj: Constants; see equations (28)–(30)

a: Defined in equation (7)

B: Formation volume factor [L³/L³]

ct: Total compressibility [L T²]

d: Dimension

Eα,β (z): Mittag-Leffler function; see equation (9)

f_j (φ): See equations (18) and (19)

h: Thickness [L]

$H_{p,q}^{m,n}\left[x|\begin{array}{c}(a_1,A_1),...,(a_p,A_p)\\ (b_1,B_1),...,(b_q,B_q)\end{array}\right]$ : Fox function

Kν(z): Modified Bessel function; second kind of order ν

k: Permeability [L²]

$\tilde{k}$ : Subdiffusive permeability; see equation (11) [L²/T^1−α]

$\mathcal{l}$ : Reference length [L]

L: Distance between the interface and the well [L]

p: Pressure [M/L/T²]

p_i : Initial pressure [M/L/T²]

p′: Logarithmic derivative [M/L/T²]

q: Rate [L³/T]

$\mathbb{R}$ : Real part of a number

r: Distance [L]

s: Laplace transform variable [1/T]

t: Time [T]

T_D : See equation (46)

Tr: See equation (31)

v(x, t): Flux [L/T]

W(α, β; t): Wright function

W: Width [L]

x : Coordinates [L]

z: A complex function

α: Exponent; see equation (10)

γ_j (x; t): Source functions; equations (18) and (19)

Γ(·): Gamma function

γ: Euler’s constant (0.5772…)

η: Diffusivity; various

$\tilde{\eta }$ : “Diffusivity”; see equation (20) [L²/T^α]

η_D : See equation (48)

λ: Mobility [L T/M]

$\tilde{\lambda }$ : Equation (11)

µ: Viscosity [M/L/T]

ϕ: Porosity [L³/L³]

ψ(x): Digamma function

$\mathrm{{\rm Y}}$ : Contact resistance; see equation (16)

${\mathrm{{\rm Y}}}_{\mathrm{D}}$ : See equation (47)

${\mathrm{{\rm Y}}}_{\mathrm{r}}$ : See equation (32)

Subscripts

D: Dimensionless

F: Fault

i: Initial

j: Zone with distinct properties

w: Well, wellbore

Superscripts

−: Laplace transform

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 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

$$<r^2>\sim t^a,$$(1)with $r$ being the distance, $t$ being time and where $a$ $\ne 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 (Prats and 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, $\gamma (\mathit{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):

$$\gamma (x,y;t)=\mathrm{\Pi }⊛S+\mathrm{\Pi }⊛\mathrm{S\text{'}},$$(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 ${\mathcal{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′). 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, $p(x,t)$, is assumed to be uniform at the value of $p_i$. The transient diffusion equation for this configuration is

$$\nabla \cdot \left[v\right(\mathit{x},t\left)\right]=\varphi \left(\mathit{x}\right)c_t\frac{\mathrm{\partial }p(\mathit{x},t)}{\mathrm{\partial }t},$$(3)where the symbols $v(\mathit{x},t)$ and $\varphi \left(\mathit{x}\right)$ 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,

$$p(\mathit{x};0)=p_i, \hspace{1em}\mathit{x}\in \mathrm{\Omega };$$(4)and also

$$p\left(\right\{x,y\}\to \mathrm{\infty };t)=p_i.$$(5)

Assuming flow to a line-source well, the condition at a wellbore for a specified withdrawal rate, q, is

$$\underset{r\to 0}{\lim }\left[2\pi h\mathrm{rv}\left(r\right)\right]=\mathrm{qB},$$(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

$$v_{x_1}(x,y;t)=a\left[p_1\right(x,y;t)-p_2(x,y;t\left)\right],\hspace{1em} y=0,$$(7)and

$$v_{x_1}(x,y;t)=v_{x_2}(x,y;t), \hspace{1em}y=0,$$(8)where the symbol $a=k_{\mathrm{F}}/\left(\mu w_{\mathrm{F}}\right)$ 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.

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 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, Δ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 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, η, follows the relation $\eta \sim r^{-(d_{\mathrm{w}}-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, Δ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 $\psi \left(t\right)=(t^{\omega }/C)E_{\omega ,\omega }(-t^{\omega }/C)$, where $E_{\omega ,\omega }\left(x\right)$ is the Mittag-Leffler function discussed in many texts (Erdelyi et al., 1955; Mainardi, 2010; Povstenko, 2015; Uchaikin, 2013):

$$E_{\alpha ,\beta }\left(z\right)=\sum _{n = 0}^{\mathrm{\infty }}\mathrm{}\frac{z^n}{\mathrm{\Gamma }(\beta +\alpha n)},\hspace{1em} \alpha , \beta >0,\hspace{0.5em}z \in C,$$(9)where C represents the complex plane and where $\mathrm{\Gamma }(\cdot )$ 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 power-law 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, $\psi \left(t\right)$.

The CTRW formulation leads to a flux law of the following form:

$$v(\mathit{x},t)=-\tilde{\lambda }\frac{{\mathrm{\partial }}^{1-\alpha }}{\mathrm{\partial }{t}^{1-\alpha }}\left[\nabla p(\mathit{x},t)\right],$$(10)with

$$\tilde{\lambda }=\frac{\tilde{k}}{\mu }.$$(11)

Here α is the subdiffusion coefficient with 0 < α < 1, μ is the fluid viscosity, $\tilde{k}$ is the permeability corresponding to α and ${\mathrm{\partial }}^{\alpha }f\left(t\right)/\mathrm{\partial }{t}^{\alpha }$, is the fractional derivative defined in Caputo (1967) sense; thus

$$\frac{{\mathrm{\partial }}^{\alpha }}{\mathrm{\partial }{t}^{\alpha }}f\left(t\right)=\frac{1}{\mathrm{\Gamma }\left(1-\alpha \right)}{\int }_0^t\mathrm{}\mathrm{d}t(t-\mathrm{t\text{'}}{)}^{-\alpha }\frac{\mathrm{\partial }}{\mathrm{\partial }\mathrm{t\text{'}}}f(\mathrm{t\text{'}}).$$(12)

For α = 1, equation (10) reverts to the classical Darcy equation. The Laplace transformation of the Caputo fractional operator, ${{}_0{}^{c}D}_t^{\alpha }f\left(t\right)$, is given by

$${\int }_0^{\mathrm{\infty }}\mathrm{}{\mathrm{e}}^{-\mathrm{st}} {{}_0{}^{c}D}_t^{\alpha }f\left(t\right)\mathrm{d}t=s^{\alpha }f\left(s\right)-s^{\alpha -1}f\left(0\right),\hspace{1em} 0\le \alpha <1.$$(13)

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 $\Delta p_j(\mathit{x},t)=p_i-p_j(\mathit{x},t)$, 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

$$\nabla \cdot [{\tilde{\lambda }}_j\nabla \Delta p_j(\mathit{x},t\left)\right]={\varphi }_j{c}_{\mathrm{tj}}\frac{{\mathrm{\partial }}^{{\alpha }_j}}{\mathrm{\partial }{t}^{{\alpha }_j}}\left[\Delta p_j(\mathit{x},t)\right].$$(14)

Similarly, equations (6)–(8) lead, respectively, to

$${\left[r\frac{{\mathrm{\partial }}^{1-{\alpha }_1}}{\mathrm{\partial }{t}^{1-{\alpha }_1}}\left(\frac{\mathrm{\partial }p}{\mathrm{\partial }r}\right)\right]}_{r\to 0}=\frac{\mathrm{qB\mu }}{2\pi {\tilde{k}}_{j}h},\hspace{1em} t>0,$$(15)

$$-\mathrm{{\rm Y}}\frac{{\mathrm{\partial }}^{1-{\alpha }_1}}{\mathrm{\partial }{t}^{1-{\alpha }_1}}\left[\frac{\mathrm{\partial }\mathrm{\Delta }{p}_1}{\mathrm{\partial }x}(x,y;t)\right]=\Delta p_2(x,y;t)-\Delta p_1(x,y;t),\hspace{1em} y=0,$$(16) $\frac{{\tilde{k}}_1{w}_{\mathrm{F}}}{k_F\mathcal{l}}$where $\mathrm{{\rm Y}}={\tilde{k}}_1{w}_F/k_F$, and

$${\tilde{\lambda }}_1\frac{{\mathrm{\partial }}^{1-{\alpha }_1}}{\mathrm{\partial }{t}^{1-{\alpha }_1}}\left(\frac{\mathrm{\partial }\mathrm{\Delta }{p}_1}{\mathrm{\partial }x}\right)={\tilde{\lambda }}_2\frac{{\mathrm{\partial }}^{1-{\alpha }_2}}{\mathrm{\partial }{t}^{1-{\alpha }_2}}\left(\frac{\mathrm{\partial }\mathrm{\Delta }{p}_2}{\mathrm{\partial }x}\right), \hspace{1em}y=0.$$(17)

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₀ (x), that reflects the free subdiffusive propagator, $\mathcal{G}(\mathit{r};t)$, in ${\mathcal{R}}^2$. Accordingly, for Zone 1 and Zone 2 we write instantaneous solutions, γ_j (x, y); j = 1, 2, respectively, in terms of their Laplace transformation, ${\bar{\gamma }}_j(x,y)$ with the source located at ($0,y\mathrm{\prime }$) to be

$$\begin{array}{ll}& {\bar{\gamma }}_1(x,y)=\frac{1}{2\pi s^{1-{\alpha }_1}{\tilde{\eta }}_1}{K}_0\left(\sqrt{s^{{\alpha }_1}/{\tilde{\eta }}_1} r_1\right)\\ & +{\int }_0^{\mathrm{\infty }}\mathrm{}{f}_1\left(\phi \right){\mathrm{e}}^{-y{q}_1}\cos \left(\phi x\right)\mathrm{d}\phi ,\end{array}$$(18)and

$${\bar{\gamma }}_2(x,y)={\int }_0^{\mathrm{\infty }}\mathrm{}{f}_2\left(\phi \right){\mathrm{e}}^{y{q}_2}\cos \left(\phi x\right)\mathrm{d}\phi ,$$(19)where ${\tilde{\eta }}_j$, the “diffusivity” of the medium, is given by

$${\tilde{\eta }}_j=\frac{{\tilde{k}}_j}{{\varphi }_j\mathrm{c\mu }},$$(20)and

$$r_1^2=x^2+(y-\mathrm{y\text{'}}{)}^2$$(21)with the symbol q_j given by

$$q_j=\sqrt{{\phi }^2+\frac{s^{{\alpha }_j}}{{\tilde{\eta }}_j}}.$$(22)

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\left(x\right)$ is

$$K_0\left[a(x^2+y^2{)}^{\frac{1}{2}}\right]={\int }_0^{\mathrm{\infty }}\mathrm{}\exp \left[-y(u^2+a^2{)}^{\frac{1}{2}}\right]\frac{\cos \left(\mathrm{xu}\right)}{(u^2+a^2{)}^{\frac{1}{2}}}\mathrm{d}u.$$(23)

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₁ (φ) and f₂(φ), we obtain the following expressions:

$$\begin{array}{ll}& {\bar{\gamma }}_1\left(\mathit{x}\right)=\frac{1}{2\pi s^{1-{\alpha }_1}{\tilde{\eta }}_1}[K_0\left(\sqrt{s^{{\alpha }_1}/{\tilde{\eta }}_1} r_1\right)\\ & +{\int }_0^{\mathrm{\infty }}\mathrm{}{A}_1{\mathrm{e}}^{-q_1(y+y_1^{\mathrm{\prime }})}\cos \left(\phi x\right)\mathrm{d}\phi ]\end{array}$$(24)or

$$\begin{array}{ll}& {\bar{\gamma }}_1\left(\mathit{x}\right)=\frac{1}{2\pi s^{1-{\alpha }_1}{\tilde{\eta }}_1}[K_0\left(\sqrt{s^{{\alpha }_1}/{\tilde{\eta }}_1} r_1\right)-K_0\left(\sqrt{s^{{\alpha }_1}/{\tilde{\eta }}_1} r_2\right)\\ & +2{\int }_0^{\mathrm{\infty }}\mathrm{}{A}_2 {\mathrm{e}}^{-q_1(y+y_1^{\mathrm{\prime }})}\cos \left(\phi x\right)\mathrm{d}\phi ],\end{array}$$(25)and,

$${\bar{\gamma }}_2\left(\mathit{x}\right)=\frac{1}{\pi s^{1-{\alpha }_1}{\tilde{\eta }}_1}{\int }_0^{\mathrm{\infty }}\mathrm{}{A}_3 {\mathrm{e}}^{(q_{2}y-q_1\mathrm{y\text{'}})}\cos \left(\phi x\right)\mathrm{d}\phi ,$$(26)where

$$r_2^2=x^2+(y+\mathrm{y\text{'}}{)}^2.$$(27)

The functions $A_1$, $A_2$ and $A_3$ are given, respectively, by

$$A_1=\frac{q_1-(1-{\mathrm{{\rm Y}}}_r{q}_1)T_r{q}_2}{q_1[q_1+(1+{\mathrm{{\rm Y}}}_r{q}_1\left)T_r{q}_2\right]},$$(28)

$$A_2=\frac{1+{\mathrm{{\rm Y}}}_r{T}_r{q}_2}{q_1+(1+{\mathrm{{\rm Y}}}_r{q}_1)T_r{q}_2},$$(29)and

$$A_3=\frac{1}{q_1+(1+{\mathrm{{\rm Y}}}_r{q}_1)T_r{q}_2}.$$(30)

Here

$$T_r=\frac{T_2}{T_1}{s}^{{\alpha }_1-{\alpha }_2},$$(31)

$${\mathrm{{\rm Y}}}_r=\frac{{\tilde{k}}_1{w}_{\mathrm{F}}}{k_F}{s}^{1-{\alpha }_1}$$(32)and $T_j={\tilde{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 Δp (x, y, z; t) (Raghavan and Ozkan, 1994) is

$$\overline{\Delta p}\left(\mathit{x}\right)=\frac{1}{\varphi c}{\int }_{\tilde{S}}\mathrm{}\overline{\stackrel{̃}{q}}\left(x^{\text{'}},y^{\text{'}}\right)\bar{\gamma }\left(x-y^{\text{'}},y-y^{\text{'}}\right)\mathrm{d}\mathrm{S\text{'}}.$$(33)

Here (x′, y′) is the location of the source, (x, y) any location in the porous rock and the strength of the continuous line-source is $\tilde{q}(x\mathrm{\prime }, y\mathrm{\prime };t)/\left(\varphi c\right)$ with dS′ 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

$$\overline{\Delta p}\left(\mathit{x}\right)=\frac{1}{s}\frac{\tilde{q}}{\varphi c}{\int }_{\tilde{S}}\mathrm{}\bar{\gamma }\left(x-\mathrm{x\text{'}},y-\mathrm{y\text{'}},z-\mathrm{z\text{'}}\right)\mathrm{d}\mathrm{S\text{'}}.$$(34)

Thus on using equation (34) and the expressions for ${\bar{\gamma }}_j(x,y)$ derived above, the expressions for $\Delta p_j(x,y;t)$ in terms of their respective Laplace transformations, ${\overline{\Delta p}}_j(x,y)$, are

$$\begin{array}{ll}& {\overline{\Delta p}}_1(x,y)=\frac{1}{s^{2-{\alpha }_1}}\frac{q_{\mathrm{w}}\mathrm{B\mu }}{2\pi T_1}\left[K_0\left(\sqrt{s^{{\alpha }_1}/{\tilde{\eta }}_1}\hspace{0.5em}{r}_1\right)+{\int }_0^{\mathrm{\infty }}\mathrm{}{A}_1 {\mathrm{e}}^{-q_1\left(y+{y\mathrm{\prime }}_1\right)}\cos \left(\phi x\right)\mathrm{d}\phi \right]\\ & \end{array},$$(35)or

$$\begin{array}{ll}& {\overline{\Delta p}}_1(x,y)=\frac{1}{s^{2-{\alpha }_1}}\frac{\mathrm{qB\mu }}{2\pi T_1}\left[K_0\left(\sqrt{s^{{\alpha }_1}/{\tilde{\eta }}_1} r_1\right)- K_0\left(\sqrt{s^{{\alpha }_1}/{\tilde{\eta }}_1} r_2\right)+2{\int }_0^{\mathrm{\infty }}\mathrm{}{A}_2 {\mathrm{e}}^{-q_1(y+{y\mathrm{\prime }}_1)}\cos \left(\phi x\right)\mathrm{d}\phi \right]\\ & \end{array},$$(36)in Zone 1 (y > 0), and

$${\overline{\Delta p}}_2(x,y)=\frac{1}{s^{2-{\alpha }_1}}\frac{\mathrm{qB\mu }}{\pi T_1}{\int }_0^{\mathrm{\infty }}\mathrm{}{A}_3 {\mathrm{e}}^{(q_{2}y-q_1{y\mathrm{\prime }}_1)}\cos \left(\phi x\right)\mathrm{d}\phi ,$$(37)in Zone 2 (y < 0).

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:

$$p_{\mathrm{D}}({\mathit{x}}_{\mathbf{D}};t_{\mathrm{D}})=\frac{2\pi {\tilde{k}}_{1}h}{\mathrm{qB\mu }}{\left(\frac{{\tilde{\eta }}_1}{{\mathcal{l}}^2}\right)}^{\frac{1-{\alpha }_1}{{\alpha }_1}}\Delta p(\mathit{x};t),$$(38)

$$t_{\mathrm{D}}=\frac{{\tilde{\eta }}_1}{{\mathcal{l}}^2}{t}^{{\alpha }_1},$$(39)and

$${\mathit{x}}_{\mathbf{D}}=\frac{\mathit{x}}{\mathcal{l}},$$(40)where $\mathcal{l}$ is the reference length.

Using these definitions, we obtain the following expression for the pressure distributions, ${\bar{p}}_{{\mathrm{D}}_j}\left({\mathit{x}}_{\mathbf{D}}\right)$; j = 1 and 2, in terms of the Laplace transformation for Zones 1 and 2 to be, respectively,

$$\begin{array}{ll}& {\bar{p}}_{{\mathrm{D}}_1}\left({\mathit{x}}_{\mathbf{D}}\right)={\left(\frac{{\tilde{\eta }}_1}{{\mathcal{l}}^2}\right)}^{\frac{1-{\alpha }_1}{{\alpha }_1}}\frac{1}{s^{2-{\alpha }_1}}\left[K_0\left(\sqrt{s^{{\alpha }_1}\mathcal{l}/{\tilde{\eta }}_1} r_{\mathrm{D}1}\right)-K_0\left(\sqrt{s^{{\alpha }_1}{\mathcal{l}}^2/{\tilde{\eta }}_1} r_{{\mathrm{D}}_2}\right)+2{\int }_0^{\mathrm{\infty }}\mathrm{}{A}_2 {\mathrm{e}}^{-q_1\left(y_{\mathrm{D}}+y_{\mathrm{D}}^{\mathrm{\prime }}\right)}\cos \left(\phi x_{\mathrm{D}}\right)\mathrm{d}\phi \right]\\ & \end{array},$$(41)for y_D > 0 and

$${\bar{p}}_{{\mathrm{D}}_2}\left({\mathit{x}}_{\mathbf{D}}\right)={\left(\frac{{\tilde{\eta }}_1}{{\mathcal{l}}^2}\right)}^{\frac{1-{\alpha }_1}{{\alpha }_1}}\frac{2}{s^{2-{\alpha }_1}}{\int }_0^{\mathrm{\infty }}\mathrm{}{A}_3 {\mathrm{e}}^{(q_2{y}_{\mathrm{D}}-q_1{\mathrm{y\text{'}}}_{\mathrm{D}})}\cos \left(\phi x_{\mathrm{D}}\right)\mathrm{d}\phi ,$$(42)for $y_{\mathrm{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 $\mathrm{{\rm Y}}$, ${\tilde{\eta }}_r$ and T_r where

$${\eta }_r=\frac{{\tilde{\eta }}_1}{{\tilde{\eta }}_2}.$$(43)

But the A_j’s; that is, both q_j’s, and the three variables, T_r, ${\tilde{\eta }}_r$, and ${\mathrm{{\rm Y}}}_r$, must, of course, be rescaled. The scalings are:

$$\begin{array}{ll}& T_r=\frac{T_2}{T_1}{s}^{{\alpha }_1-{\alpha }_2}\\ & =T_{\mathrm{D}}{\left(\frac{{\mathcal{l}}^2{s}^{{\alpha }_1}}{{\tilde{\eta }}_1}\right)}^{1-{\alpha }_{\mathrm{D}}},\end{array}$$(44)

$${\mathrm{{\rm Y}}}_r={\mathrm{{\rm Y}}}_{\mathrm{D}}{\left(\frac{{\mathcal{l}}^2{s}_1^{\alpha }}{{\tilde{\eta }}_1}\right)}^{\frac{1-{\alpha }_1}{{\alpha }_1}},$$(45)with

$$T_{\mathrm{D}}=\frac{T_2}{T_1}{\left(\frac{{\tilde{\eta }}_1}{{\mathcal{l}}^2}\right)}^{\left(1-{\alpha }_{\mathrm{D}}\right)},$$(46)

$${\mathrm{{\rm Y}}}_{\mathrm{D}}=\frac{{\tilde{k}}_1{w}_{\mathrm{F}}}{k_F\mathcal{l}}{\left(\frac{{\tilde{\eta }}_1}{{\mathcal{l}}^2}\right)}^{\frac{1-{\alpha }_1}{{\alpha }_1}},$$(47)

$${\eta }_{\mathrm{D}}=\frac{{\tilde{\eta }}_1^{{\alpha }_{\mathrm{D}}}}{{\tilde{\eta }}_2}{\mathcal{l}}^{2\left(1-{\alpha }_{\mathrm{D}}\right)},$$(48)and

$${\alpha }_{\mathrm{D}}=\frac{{\alpha }_2}{{\alpha }_1}.$$(49)

Scalings involving the ratio of α_j’s enter into the picture because we consider subdiffusive flow; in accordance with equation (10); the dimension of $\tilde{k}$ 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 $(\varphi c_{t}h{)}_{\mathrm{D}}$ if desired, where

$$(\varphi c_{t}h{)}_{\mathrm{D}}=\frac{(\varphi c_{t}h{)}_1}{(\varphi c_{t}h{)}_2}.$$(50)

A number of approximate solutions may be derived if q₁ = q₂, which implies ${\tilde{\eta }}_1$ and ${\tilde{\eta }}_2$ and also α₁ = α₂. In the following we use $\mathcal{l}=r_{\mathrm{w}}$, where r_w is the well radius. In this situation for Zone 1 and Zone 2, respectively, we have

$$\begin{array}{ll}{\bar{p}}_{{\mathrm{D}}_1}\left({\mathit{x}}_{\mathbf{D}}\right)={\left(\frac{{\tilde{\eta }}_1}{r_{\mathrm{w}}^2}\right)}^{\frac{1-{\alpha }_1}{{\alpha }_1}}\frac{1}{s^{2-{\alpha }_1}}\left[K_0\left(\sqrt{s^{{\alpha }_1}{r}_{\mathrm{w}}^2/{\tilde{\eta }}_1} r_{{\mathrm{D}}_1}\right)-\frac{1-T_{\mathrm{D}}}{1+T_{\mathrm{D}}}{K}_0\left(\sqrt{s^{{\alpha }_1}{r}_{\mathrm{w}}^2/{\tilde{\eta }}_1} r_{{\mathrm{D}}_2}\right)\right],& \end{array}$$(51)and

$${\bar{p}}_{{\mathrm{D}}_2}\left({\mathit{x}}_{\mathbf{D}}\right)={\left(\frac{{\tilde{\eta }}_1}{r_{\mathrm{w}}^2}\right)}^{\frac{1-{\alpha }_1}{{\alpha }_1}}\frac{2}{(1+T_r)}\frac{1}{s^{2-{\alpha }_1}}{K}_0\left(\sqrt{s^{{\alpha }_1}{r}_{\mathrm{w}}^2/{\tilde{\eta }}_1} r_{{\mathrm{D}}_1}\right).$$(52)

The right-hand sides of equations (51) and (52) may be inverted with the expression

$$2L^{-1}\left[s^{-\mathrm{\varrho }}{K}_{\nu }\right(z{s}^{\varsigma }\left)\right]=t^{\varrho -1}{H}_{\mathrm{1,2}}^{\mathrm{2,0}}\left[{\frac{z^2{t}^{-2\varsigma }}{4}|}_{\left(\frac{\nu }{2},1\right), \left(-\frac{\nu }{2},1\right)}^{\varrho ,2\varsigma }\right],$$(53)where the symbol $H_{p,q}^{m,n}\left[{x|}_{(b_1,B_1),...,(b_q,B_q)}^{(a_1,A_1),...,(a_p,A_p)}\right]$ is the Fox function or the H-function, $\mathbb{R}\left(\varrho \right)>0$, $\mathbb{R}\left(z^2\right)>0$, $\mathbb{R}\left(s\right)>0$ as noted in Saxena et al. (2006), and thus we have the following result for $p_{{\mathrm{D}}_1}({\mathit{x}}_{\mathbf{D}};t_{\mathrm{D}})$, namely

$$p_{{\mathrm{D}}_1}({\mathit{x}}_{\mathbf{D}};t_{\mathrm{D}})=\frac{1}{2}{t}_{\mathrm{D}}^{\frac{1-{\alpha }_1}{{\alpha }_1}}\sum _{j=0}^{j=1}\mathrm{}\left\{{\left(\frac{1-T_{\mathrm{D}}}{1+T_{\mathrm{D}}}\right)}^j{H}_{\mathrm{1,2}}^{\mathrm{2,0}}\left[{\frac{r_{{\mathrm{D}}_{j+1}}^2}{4t_{\mathrm{D}}}|}_{\left(\mathrm{0,1}\right),\left(\mathrm{0,1}\right)}^{2-{\alpha }_1,{\alpha }_1}\right]\right\},$$(54)and for $p_{{\mathrm{D}}_2}({\mathit{x}}_{\mathbf{D}};t_{\mathrm{D}})$, namely

$$p_{{\mathrm{D}}_2}({\mathit{x}}_{\mathbf{D}};t_{\mathrm{D}})=\frac{1}{(1+T_r)}{t}_{\mathrm{D}}^{\frac{1-{\alpha }_1}{{\alpha }_1}}\left\{H_{\mathrm{1,2}}^{\mathrm{2,0}}\left[{\frac{r_{{\mathrm{D}}_1}^2}{4t_{\mathrm{D}}}|}_{\left(\mathrm{0,1}\right),\left(\mathrm{0,1}\right)}^{2-{\alpha }_1,{\alpha }_1}\right]\right\}.$$(55)

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\left(x\right)$ as x → 0 (implicitly we consider situations where s → 0), namely using

$$K_0\left(z\right)=-\left(\ln \frac{z}{2}+\gamma \right)+\frac{z^2}{4}\left[1-\left(\gamma +\ln \frac{z}{2}\right)\right]+O(z^4\ln z),$$(56)where γ is Euler’s constant. The long-time 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)

$$s^{-\nu }\ln s=\left[\mathrm{\Gamma }\right(\nu \left){]}^{-1}{t}^{\nu -1}\right[\psi \left(\nu \right)-\ln t],$$(57)where ψ(⋅) is the Digamma function, and we therefore have

$$\begin{array}{ll}& p_{{\mathrm{D}}_1}({\mathit{x}}_{\mathbf{D}},t_{\mathrm{D}})=\frac{1}{2}\frac{1}{\mathrm{\Gamma }(2-{\alpha }_1)}{t}_{\mathrm{D}}^{\frac{1-{\alpha }_1}{{\alpha }_1}}\{\left[\ln \frac{4t_{\mathrm{D}}}{{\mathrm{e}}^{2\gamma }{r}_{\mathrm{D}1}^2}-{\alpha }_1\psi (2-{\alpha }_1)\right]\\ & +\left(\frac{1-T_{\mathrm{D}}}{1+T_{\mathrm{D}}}\right)\left[\ln \frac{4t_{\mathrm{D}}}{{\mathrm{e}}^{2\gamma }{r}_{\mathrm{D}2}^2}-{\alpha }_1\psi (2-{\alpha }_1)\right]\}.\end{array}$$(58)

Also, the logarithmic derivative, ${p\mathrm{\prime }}_{\mathrm{D}}({\mathit{x}}_{\mathbf{D}},t_{\mathrm{D}})$, is given by the expression

$$\begin{array}{ll}& {\mathrm{p\text{'}}}_{{\mathrm{D}}_1}({\mathit{x}}_{\mathbf{D}},t_{\mathrm{D}})=\frac{1}{2}\frac{1}{\mathrm{\Gamma }(1-{\alpha }_1)}{t}_{\mathrm{D}}^{\frac{1-{\alpha }_1}{{\alpha }_1}}\{\left[\ln \frac{4t_{\mathrm{D}}}{{\mathrm{e}}^{2\gamma }{r}_{\mathrm{D}1}^2}-{\alpha }_1\psi (1-{\alpha }_1)\right]\\ & +\left(\frac{1-T_{\mathrm{D}}}{1+T_{\mathrm{D}}}\right)\left[\ln \frac{4t_{\mathrm{D}}}{{\mathrm{e}}^{2\gamma }{r}_{\mathrm{D}2}^2}-{\alpha }_1\psi (1-{\alpha }_1)\right]\}.\end{array}$$(59)

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 semi-logarithmic 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 ${\tilde{\eta }}_2$ different from ${\tilde{\eta }}_1$ at this time. Although this is the case, we know that for classical diffusion (${\alpha }_j=1$) the influence of ${\tilde{\eta }}_r$ 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 ${\tilde{\eta }}_r\ne 1$ provided that the ${\alpha }_j$’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, ${p\mathrm{\prime }}_{\mathrm{D}}({\mathit{x}}_{\mathbf{D}};t_{\mathrm{D}})$, as these curves best demonstrate the characteristics of the rock under consideration. For each of the solutions presented, the pressure response, $p_{\mathrm{D}}({\mathit{x}}_{\mathbf{D}};t_{\mathrm{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.

2.1 Perfect contact

Figure 1 presents derivative, $p_{\mathrm{wD}}^{\mathrm{\prime }}\left(t_{\mathrm{D}}\right)$, responses for a number of situations; the intent is to demonstrate the structure of the solutions and the influence of properties, specifically the exponent ${\alpha }_j$ and the ratio $T_{\mathrm{D}}$. As one of our interests here is to demonstrate the accuracy of the solutions we assume the diffusivity ratio, ${\eta }_{\mathrm{D}}$, to be 1. All properties considered are noted in Figure 1. The symbol $L_{\mathrm{D}}$ that is the dimensionless distance to the fault is defined by

$$L_{\mathrm{D}}=\frac{L}{\mathcal{l}}.$$(60)

Simulations obtained through equation (36) for ${\mathrm{{\rm Y}}}_{\mathrm{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 α_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 α₂ 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 α₂ 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 ${\eta }_{\mathrm{D}}$ and ${\alpha }_j$. The solutions are for $T_{\mathrm{D}}$ of $100$ and ${\alpha }_1$ of $0.9$. The top set corresponds to dimensionless pressure, $p_{\mathrm{wD}}$, response, and the bottom set represents the corresponding derivative, ${p\mathrm{\prime }}_{\mathrm{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 ${\eta }_{\mathrm{D}}$ and/or ${\alpha }_2$ as soon as the region beyond the interface influences the well response. As in Figure 1, we assume $L_{\mathrm{D}}$ 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 α₂ of 0.9 and 0.5 and three values of η_D: 5 (squares), 1 (unbroken lines) and 0.2 (squares). Considering the $p_{\mathrm{wD}}\left(t_{\mathrm{D}}\right)$ solutions for α₂ of 0.5, namely Set A, we see that the influence of ${\eta }_{\mathrm{D}}$ 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 $p_{\mathrm{wD}}^{\mathrm{\prime }}\left(t_{\mathrm{D}}\right)$ 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.

2.2 Influence of a contact resistance

Figure 3 presents the influence of a fault, ${\mathrm{{\rm Y}}}_{\mathrm{D}}=10^3$, located at a distance of $L_{\mathrm{D}}$ of $100$ at an observation point, $x_{\mathrm{D}}=60,y_{\mathrm{D}}=40$ with both the well and the observation location on the same side of the interface. The influences of T_D, η_D and α₂ for fixed values of Zone 1 properties with α₁ 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 ${\eta }_{\mathrm{D}}$ of $8$; the one identified as A is for an ${\eta }_{\mathrm{D}}=1$. Dimensionless times in Figure 3 are normalized with respect $L_{\mathrm{D}}=100$; that is, all solutions shown here are expressed as $t_{\mathrm{D}}/L_{\mathrm{D}}^2$. 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_{\mathrm{D}}$; we consider three values of $T_{\mathrm{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 $(\varphi c_t{)}_1/(\varphi c_t{)}_2$. Ultimately all lines become flat and reflect $2/(1+T_{\mathrm{D}})$ as both $\alpha$’s are 1 (see Eq. (59)).

The situation is different if ${\alpha }_2<1$; the dashed lines in Figure 3 reflect solutions for ${\alpha }_2$ of $0.9$ for identical values of $T_{\mathrm{D}}$. These solutions deviate from the unbroken lines with the upward trend representing the fact of an additional restriction, namely that ${\alpha }_2$ is now $0.9$. A further reduction in ${\alpha }_2$ will result in an earlier deviation and the magnitude of the upward trend would be larger; see line with small dashes for $T_{\mathrm{D}}$ of $10$ which represents ${\alpha }_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 ${\eta }_{\mathrm{D}}$. The basis of this observation may be understood by comparing the curve labelled A which corresponds to ${\eta }_{\mathrm{D}}=1$ and $T_{\mathrm{D}}=1$ with its counterpart the solution for ${\eta }_{\mathrm{D}}=8$ and $T_{\mathrm{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 ${\alpha }_1=0.9$ and ${\alpha }_2=1.0$; that is, subdiffusive flow in Zone $1$ and classical flow in Zone 2 – the opposite of situations considered in Figure 3 although a slightly different value of ${\eta }_{\mathrm{D}}$ is used and the range of $T_{\mathrm{D}}$ 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 ${\alpha }_1$ from $1.0$ is $0.9$ results in a significant change in the magnitude of $p_{\mathrm{wD}}\left(t_{\mathrm{D}}\right)$. All solutions considered here except for the bottommost curve assume subdiffusive flow in Zone $1$, and, again all solutions for ${\alpha }_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 ${\alpha }_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 ${\alpha }_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 ${\alpha }_2$ decreases. In practice such steeper slopes would not reflect additional boundary effects.

3 Discussion

Before closing with a few general observations, for practicality, we note that the exponent $\alpha$ 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:

$$\mathcal{G}(\mathit{r};t)=\frac{1}{4\pi \tilde{\eta }{t}^{\alpha }}{H}_{\mathrm{1,2}}^{\mathrm{2,0}}\left[\frac{r^2}{4\tilde{\alpha }{t}^{\nu }}{|}_{\left(\mathrm{0,1}\right),\left(\mathrm{0,1}\right)}^{1-\alpha ,\alpha }\right],$$(61)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 $\mathcal{G}(\mathit{r};t)$ in the following form,

$$\mathcal{G}\left(\mathit{r};t\right)\sim \frac{{\alpha }^{(\mathrm{d\alpha }-1)/2}}{\sqrt{2-\alpha }}{\left(4\pi \tilde{\eta }{t}^{\alpha }\right)}^{-d/2}{\left(\frac{r^2{\alpha }^{\alpha }}{4\tilde{\eta }{t}^{\alpha }}\right)}^{-\frac{d(1-\alpha )}{2(2-\alpha )}}\exp \left[-(2-\alpha ){\left(\frac{r^2{\alpha }^{\alpha }}{4\tilde{\eta }{t}^{\alpha }}\right)}^{\frac{1}{(2-\alpha )}}\right],$$(62)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. $\mathcal{G}(\mathit{r};t)$ exhibits power-law behavior at long-enough times for $\alpha <1$. This expression agrees with those in Carslaw and Jaeger (1959) for $d\le 3$ if $\alpha =1$. Grebenkov (2010) arrives at the same expression as that given in equation (62) but for $\mathcal{G}(\mathit{r};t\to 0)$ by using the result

$$L^{-1}\left[s^{-\kappa }\exp \right(-r\sqrt{s^{\alpha }/\tilde{\eta }}\left)\right]\simeq \frac{{\alpha }^{\kappa +1/2}}{\sqrt{\pi (2-\alpha )}}{\left(\frac{r^2{\alpha }^{\alpha }}{4\tilde{\eta }{t}^{\alpha }}\right)}^{\frac{(\kappa +1/2)}{(2-\alpha )}}{t}^{-(\kappa +1)}\exp \left[-(2-\alpha ){\left(\frac{r^2{\alpha }^{\alpha }}{4\tilde{\eta }{t}^{\alpha }}\right)}^{\frac{1}{(2-\alpha )}}\right],$$(63)with $\kappa =(d+1)/4-1$. Incidentally, the left-hand side of equation (63) may also be expressed in terms of the Wright function, $W(\alpha ,\beta ;t)$, as indicated in Povstenko (2015) by

$$L^{-1}\left[s^{-\beta }\exp \right(-\lambda s^{\alpha }\left)\right]=t^{\beta -1}W(-\alpha ,\beta ;-\lambda t^{-\alpha }), \hspace{1em}0<\alpha <1, \lambda >0.$$(64)

The moments of an even order which would also be of interest are (Uchaikin, 2013); see equation (1):

$${\int }_{{\mathcal{R}}^d}\mathrm{}\left|\mathit{x}{|}^{2n}\mathcal{G}\right(\left|\mathit{x}\right|,t)\mathrm{d}\mathit{x}=\frac{\mathrm{\Gamma }(n+d/2)\mathrm{\Gamma }(n+1)}{\mathrm{\Gamma }(d/2)\mathrm{\Gamma }(\alpha n+1)}\eta t^{\mathrm{n\alpha }}, \hspace{1em}n=\mathrm{0,1},2,....$$(65)

Variations of $\mathcal{G}(\mathit{r};t)$ as a function of r are shown in Figure 5 for three values of $\alpha$: $1$, $0.9$, and 0.8. These values were computed independent of the expressions (62) and (63). The characteristic bell shape is evident only for $\alpha =1$. Note that $\mathcal{G}(\mathit{r};t)$ may be estimated by equation (62) for $r=0$ for $\alpha$ of $1$.

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 ${\alpha }_j$ (or at least ${\alpha }_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, $\Delta p\mathrm{\prime }$, as this curve when combined with the pressure response, $\Delta 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.

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 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.

References

All Figures