ogst
Home All issues Volume 76 (2021) Oil Gas Sci. Technol. – Rev. IFP Energies nouvelles, 76 (2021) 78 Full HTML
Open Access
Issue
Oil Gas Sci. Technol. – Rev. IFP Energies nouvelles
Volume 76, 2021
Article Number 78
Number of page(s) 10
DOI https://doi.org/10.2516/ogst/2021058
Published online 20 December 2021

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

Licence Creative CommonsThis 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

c : Compressibility [L T2/M]

Eα,β (z): Mittag–Leffler function; see (11)

h : Thickness [L]

k : Permeability [L2]

kα,β : Fractional permeability; see (2) [Lβ+1T1−α]

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

L $ \mathcal{L}$ : Laplace transform

n : Number of dimensions

p : Pressure [M/L/T2]

p′: Logarithmic derivative [M/L/T2]

q : Rate [L3/T]

t : Time [T]

ℝ: Real part of a number

rw : Radius of well bore [L]

z : A complex function

α : Exponent

β : Exponent

Γ(x): Gamma function

φα,β (x,t): Green’s function

Ψ(s): α-stable subordinator

η : Diffusivity; various

η ̃ $ \stackrel{\tilde }{\eta }$ : Fractional diffusivity; see (6) [Lβ+1T1−α]

λ : Mobility [L T/M]

μ : Viscosity [M/L/T]

ϕ : Porosity [L3/L3]

Subscripts

D : Dimensionless

i : Coordinate, initial

w: Well bore

Superscript

–: Laplace transform

1 Introduction

The ubiquity of heterogeneity in geologic media hardly needs an introduction – every occasional gardener knows firsthand. Technical ideas of the concept of heterogeneity preceded the introduction of Darcy’s flux-law with the introduction of the definition of facies, we are told, in 1838. Most working models of well performance, in practice, often degenerate to diffusion in homogeneous systems where fluid movement is presumed to follow a velocity-jump random walk following the working hypothesis that the mean square displacement or mean square deviation, 〈(Δx)2〉, is proportional to time, t; that is, 〈(Δx)2〉 ∝ t. Over time this idea has been expanded by resorting either to numerical means as summarized in Mattax and Dalton (1990) or analytically by considering Continuous Time Random Walks, CTRW, (Cortis and Knudby, 2006; Noetinger and Estebenet, 2000; Noetinger et al., 2016), fractal objects (Ali et al., 2014; Chang and Yortsos, 1990; Flamenco-Lopez and Camacho-Velazquez, 2001; Metzler et al., 1994; Park et al., 2000; Wang, 2013), waiting-time solutions (Henry et al., 2010) and Levy-stable random walks (Benson et al., 2000, 2004; Cloot and Botha, 2006; Gehlhausen, 2009; Lenormand, 1992). Many field evaluations of well responses in recent times (Bernard et al., 2006; Chu et al., 2019a, b; Cloot and Botha, 2006; Doe et al., 2014; Meerschaert et al., 2008; Raghavan and Chen, 2019; Scott et al., 2015, Sinkov et al., 2021; Thomas et al., 2005) note that the proper assessment of responses in systems producing with complex geology requires going beyond the classical transient diffusion; that is, beyond the premises of the Theis (1935), solution. They have proposed incorporation of the internal topology of the complex porous network along with actual paths where distant locations may be connected rather well whereas nearby locations are isolated because of localized barriers and dead ends which lead to mean square displacements of the form 〈(Δx)2〉 ∝ ta with a>=<1$ a\begin{array}{c}>\\ =\\ < \end{array}1$. Such an approach also requires reconsideration of the flux law – a new flux law in terms of fractional derivatives; see, for example, Fomin et al. (2011) and Chen and Raghavan (2015).

In this study we explore for the first time, to our knowledge, the combined effects of super- and subdiffusion by considering pressure distributions in 2D around a line-source well by examining solutions of the fractional differential equation, ∂tαp=1/r[∂r(r∂rβp)]$ {\mathrm{\partial }}_t^{\alpha }p=1/r[{\partial }_r(r{\mathrm{\partial }}_r^{\beta }p)]$ for 0 < α, β < 1. The properties of statistical distributions that govern the influence of obstructions, restrictions, impediments, intercalations, fault-damage zones, etc. are reflected through the exponent α, and the connectivity between distant points is reflected through the exponent β. Specifically, α reflects the waiting time in a CTRW (Henry et al., 2010; Noetinger and Estebenet, 2000), P(Tn>t)≈t-α$ \mathbb{P}({T}_n>t)\approx {t}^{-\alpha }$ for 0 < α < 1, where Tn is the random waiting time at which the jump occurs, and β reflects the path taken by Levy-stable-particles through the power law tail (P(|Xn|>r)≈r-(β+1)$ \mathbb{P}(|{X}_n|>r)\approx {r}^{-(\beta +1)}$) (Benson et al., 2004). Considered individually, the exponent, α, represents the reciprocal of the anomalous diffusion coefficient (random walk dimension), dw, (Metzler et al., 1994) for diffusion on fractal objects, <r2>∼t2/dw$ < {r}^2>\sim {t}^{2/{d}_{\mathrm{w}}}$, and the exponent, β, represents the Hurst index, H, for Levy-stable diffusion; H is 1/(β + 1) or <r2> ~ tH. Elaborations of these ideas with respect to transient, anomalous diffusion in porous solids are given in Raghavan and Chen (2020). For α = β = 1, the fractional equation reverts to the classical form with alternating sequences of fixed-velocity runs and reorientations. The analytical expressions of the Greens functions in Luchko (2016), Bazhlekova (2019) and Bazhlekova and Bazhlekov (2019) form the cornerstone of this work. The analysis of numerical experiments of expressions for the pressure distribution to a line-source well sheds light on behaviors yet to be reported; unusual features of pure superdiffusion in 2D which one may not expect are reported and differences with respect to 1D systems are documented. The results of the numerical experiments also establish the viability of the solutions for practical use and form the centrepiece of more complex well completions to be reported elsewhere. Recall that we consider both super- and subdiffusion behaviors.

2 The differential equation, auxillary conditions

This work, as is usual for single-phase flow, concentrates on the flow of a slightly compressible liquid to a cylindrical well which is treated as a line source; the well is considered to be located in a system that is infinite in extent. The structure considered is described by a set of scaling exponents, α and β, and we work in terms of pressure, p(r, t), permeability, kα,β, porosity, ϕ, thickness, h, fluid compressibility, c, and viscosity, μ in the cylindrical coordinate system assuming the axis of the well passes through the origin. Initially, at time t = 0, the pressure distribution is distributed uniformly at the value, pi; that is, equilibrium prevails at t = 0. Wellbore storage and skin effects are presumed to be zero.

The conservation equation of mass principle to a control volume in cylindrical coordinates for a slightly compressible, constant viscosity fluid is

1r∂∂rv(r,t)=ϕc∂∂tp(r,t),$$ \frac{1}{r}\frac{\mathrm{\partial }}{\mathrm{\partial }r}v(r,t)={\phi c}\frac{\mathrm{\partial }}{\mathrm{\partial }t}p(r,t), $$(1)

where v(r, t) is the velocity. The fractional flux law we use in this work is of the form

v(r,t)=-λα,β∂1-α∂t1-α[∂β∂rβp (r,t)],$$ v(r,t)=-{\lambda }_{\alpha,\beta }\frac{{\mathrm{\partial }}^{1-\alpha }}{\mathrm{\partial }{t}^{1-\alpha }}\left[\frac{{\mathrm{\partial }}^{\beta }}{\mathrm{\partial }{r}^{\beta }}p\enspace (r,t)\right], $$(2)

where λα,β = kα,β/μ, with kα,β being the effective “fractional permeability” with dimensions of [Lβ+1T1-α]$ [{L}^{\beta +1}{T}^{1-\alpha }]$. Both the exponents α and β are <1 with the temporal fractional derivative, ∂αf(t)/∂tα, defined in the Caputo (1967) sense, namely

∂α∂tαf(t)=1Γ(1-α)∫0tdt′(t-t′)-α∂∂t′f(t′),$$ \frac{{\mathrm{\partial }}^{\alpha }}{\mathrm{\partial }{t}^{\alpha }}f(t)=\frac{1}{\mathrm{\Gamma }\left(1-\alpha \right)}{\int }_0^t \mathrm{d}t\mathrm{\prime}(t-t\mathrm{\prime}{)}^{-\alpha }\frac{\mathrm{\partial }}{\mathrm{\partial }t\mathrm{\prime}f(t\mathrm{\prime}), $$(3)

and the spatial fractional derivative, ∂α f(r)/∂rα, in the Riemann–Liouville form, namely

∂α∂rαf(r)=1Γ(n-α)∂n∂rn[∫0rdξ(r-ξ)1-α+nf(ξ)],$$ \frac{{\mathrm{\partial }}^{\alpha }}{\mathrm{\partial }{r}^{\alpha }}f(r)=\frac{1}{\mathrm{\Gamma }\left(n-\alpha \right)}\frac{{\mathrm{\partial }}^n}{\mathrm{\partial }{r}^n}\left[{\int }_0^r \mathrm{d}\xi (r-\xi {)}^{1-\alpha +n}f(\xi )\right], $$(4)

where 0<n–1<R(α)<n$ 0<n\endash 1<\mathbb{R}(\alpha ) < n$.

As is well known fractional derivatives are convolutions; that is, nonlocal operators; for example, in [Oldham and Spanier (1974), Podlubny (1998)] equation (4) the value of the fractional derivative at the point r is a function of all points to the left of r. Expressions of the form given in equations (3) and (4) have been used previously (Fomin et al., 2011; Chen and Raghavan, 2015) but the manner in which fractional derivatives are used here is slightly different. For superdiffusive flows; that is, α = 1, Benson et al. (2004) use an expression similar to that given here. Many discussions of fractional flux laws are available; see Noetinger and Estebenet (2000), Suzuki et al. (2016), Henry et al. (2010), Su (2014), and Su et al. (2015) for subdiffusive flow and Molz et al. (2002), Garra and Salusti (2013), Kim et al. (2015) and Sapora et al. (2017) for superdiffusive flow. A number of experimental works also exist: Tao et al. (2016), Yanga et al. (2019), Zhokh and Strizhak (2018, 2019) and Chang et al. (2019). Fractional flux laws are also needed to address flow in fractal systems; see Le Mehaute and Crepy (1983), Nigmatullin (1984) and Beier (1990). Many discourses on fractional Fourier laws are also found in the literature (Angulo et al., 2000; Dassas and Duby, 1995; Gurtin and Pipkin, 1968; Moodie and Tait, 1983; Norwood, 1972). Raghavan and Chen (2019) summarize estimates of α based on a number of field studies for rocks of different types; they report the estimates in the range: 0.56 ≤ α ≤ 0.91 for β = 1. Expressions for the flux distributions are not always addressed as the preference is to use pseudo-differential operators (Balakrishnan, 1960).

On combining the conservation of mass principle and the flux law after noting that we consider flow to a well located at the origin we, as usual, arrive at the differential equation given by:

1r∂∂r(r ∂βp∂rβ)=1η̃∂αp∂tα,$$ \frac{1}{r}\frac{\mathrm{\partial }}{\mathrm{\partial }r}\left(r\enspace \frac{{\mathrm{\partial }}^{\beta }p}{\mathrm{\partial }{r}^{\beta }}\right)=\frac{1}{\stackrel{\tilde }{\eta }}\frac{{\mathrm{\partial }}^{\alpha }p}{\mathrm{\partial }{t}^{\alpha }}, $$(5)

where η̃$ \stackrel{\tilde }{\eta }$, the “fractional diffusivity” of the medium, is given as

η̃=k̃ϕctμ.$$ \stackrel{\tilde }{\eta }=\frac{\mathop{k}\limits^\tilde}{\phi {c}_t\mu }. $$(6)

We seek a solution to equation (5) subject to the following initial and boundary conditions:

p(r,t)=pi,  t→0,$$ p\left(r,t\right)={p}_i,\hspace{1em}\enspace t\to 0, $$(7)

limr→∞p(r,t)=pi,  t>0,$$ \underset{r\to \mathrm{\infty }}{\mathrm{lim}}p(r,t)={p}_i,\enspace \hspace{1em}t>0, $$(8)

and

{r∂1-α∂t1-α[∂β∂rβp(r,t)]}rw→0=qμ2πk̃α,βh,  t>0,$$ {\left\{r\frac{{\mathrm{\partial }}^{1-\alpha }}{\mathrm{\partial }{t}^{1-\alpha }}\left[\frac{{\mathrm{\partial }}^{\beta }}{\mathrm{\partial }{r}^{\beta }}p(r,t)\right]\right\}}_{{r}_w\to 0}=\frac{{q\mu }}{2\pi {\mathop{k}\limits^\tilde}_{\alpha,\beta }h},\enspace \hspace{1em}t>0, $$(9)

where q is the rate. Equation (9) specifies the condition for a line-source well; that is, lim rw → 0.

The analog of an expression corresponding to the left-hand side of equation (5) in 1D and 2D Cartesian coordinates given by the Riesz fractional derivative, RZDxαf(x)$ {{}_{{RZ}}D}_x^{\alpha }f(x)$, has been studied extensively as in Mainardi et al. (2001) and Huang and Liu (2005) for 1D, and as a fractional Laplacian, −(−Δ)αp(x, t); 0 < α < 1, as in Luchko (2019), Bazhlekova (2019) and Bazhlekova and Bazhlekov (2019) for 2D. The correspondence with classical diffusion is usually established through Fourier transforms (Gorenflo and Mainardi, 1998); it is given in a point-wise sense in Cai and Li (2019) and the fractional flux law becomes Darcy’s law in the limit. Kwaśnicki (2017) notes that 10 equivalent definitions of the fractional Laplacian have been used.

Incidentally,

∇⋅∇α-1p(x,t)=-2πcos(πα2)Γ(α+12)Γ(α+22)(-Δ)α/2p(x,t); 1<α<2.$$ \nabla \cdot {\nabla }^{\alpha -1}p({x},t)=-2\sqrt{\pi }\mathrm{cos}\left(\frac{{\pi \alpha }}{2}\right)\frac{\mathrm{\Gamma }(\frac{\alpha +1}{2})}{\mathrm{\Gamma }(\frac{\alpha +2}{2})}(-\Delta {)}^{\alpha /2}p({x},t);\hspace{1em}1<\alpha < 2. $$(10)

see Estrada-Rodriguez et al. (2018).

3 The Mittag–Leffler functions, Eα (z) and Eα,β(z)

As is inevitable in most analytical studies of anomalous diffusion, this work requires the use of the Mittag–Leffler functions, Eα (z) and Eα,β (z), which are briefly reviewed. Extensive discussions may be found in Erdelyi et al. (1955), Mainardi (2010), Uchaikin (2012) and Gorenflo et al. (2014) to name a few sources.

The function Eα,β(z)$ {E}_{\alpha,\beta }(z)$ is defined by the representation

Eα,β(z)=∑n=0∞znΓ(αn+β),  α>0, z ∈C,$$ {E}_{\alpha,\beta }(z)=\sum_{n=0}^{\mathrm{\infty }} \frac{{z}^n}{\mathrm{\Gamma }(\alpha n+\beta )},\enspace \hspace{1em}\alpha >0,\hspace{0.5em}z\enspace \in \mathbb{C}, $$(11)

where C$ \mathbb{C}$ denotes the complex plane. Implicit in all of the following is that Eα,1(z)=Eα(z)$ {E}_{\alpha,1}(z)={E}_{\alpha }(z)$.

Of particular interest to us is the Laplace transformation of Eα,β(z)$ {E}_{\alpha,\beta }(z)$. From the references cited, we have

L[tβ-1Eα,β(±atα)]=sα-βsα∓a; Re s>0, |as-α|<1,$$ \mathcal{L}\left[{t}^{\beta -1}{E}_{\alpha,\beta }\left(\pm a{t}^{\alpha }\right)\right]=\frac{{s}^{\alpha -\beta }}{{s}^{\alpha }\mp a};\hspace{1em}\mathrm{Re}\enspace s>0,\hspace{0.5em}|a{s}^{-\alpha }| < 1, $$(12)

where L$ \mathcal{L}$ denotes the Laplace transformation, s is the Laplace transform variable, and both α and β are real. Thus, for β = 1 we have

L[Eα(±axα)]=sα-1sα∓a; Re s>0, |as-α|<1.$$ \mathcal{L}[{E}_{\alpha }\left(\pm a{x}^{\alpha }\right)]=\frac{{s}^{\alpha -1}}{{s}^{\alpha }\mp a};\hspace{1em}\mathrm{Re}\enspace s>0,\hspace{0.5em}|a{s}^{-\alpha }| < 1. $$(13)

The following asymptotic expressions for Eα,β(z)$ {E}_{\alpha,\beta }(z)$ and Eα,β(-z)$ {E}_{\alpha,\beta }\left(-z\right)$ are often needed: For 0 < α < 2 and πα/2<ξ<min [π,πα]$ \pi \alpha /2<\xi < \mathrm{min}\enspace [\pi,\pi \alpha ]$ with z ϵ C$ z\enspace \epsilon \enspace \mathbb{C}$ we have

Eα,β(z)=1αz1-βαexp(z1α)-∑k=1∞z-kΓ(β-kα),  |z|→∞, |arg z|<ξ,$$ {E}_{\alpha,\beta }(z)=\frac{1}{\alpha }{z}^{\frac{1-\beta }{\alpha }}\mathrm{exp}({z}^{\frac{1}{\alpha }})-\sum_{k=1}^{\mathrm{\infty }} \frac{{z}^{-k}}{\mathrm{\Gamma }(\beta -k\alpha )},\hspace{1em}\enspace |z|\to \mathrm{\infty },\hspace{0.5em}|\mathrm{arg}\enspace z| < \xi, $$(14)

and

Eα,β(-z)=z-2Γ(β-2α); β-α=0,-1,-2… .$$ {E}_{\alpha,\beta }\left(-z\right)=\frac{{z}^{-2}}{\mathrm{\Gamma }(\beta -2\alpha )};\hspace{1em}\beta -\alpha =0,-1,-2\dots \enspace. $$(15)

We compute the Mittag–Leffler functions by the algorithm in Gorenflo et al. (2002) and Podlubny (2005).

4 Subordination, solution

The principle of subordination is a standard procedure for solving problems of the type considered in this work (Bazhlekova and Bazhlekov, 2019; Benson et al., 2000; Luchko, 2019; Raghavan and Chen, 2015; Stanislavsky et al., 2014). The procedure recasts the classical Green’s (or source) function in the classical diffusion equation, often called the parent process as in Benson et al. (2000), by expressing it in terms of an operational time, τ; that is, the classical solution is now expressed as ζ [x,Sα(t)]$ \zeta \enspace [{x},{S}_{\alpha }(t)]$ (Sα(t)≡τ)$ {S}_{\alpha }(t)\equiv \tau )$ rather than ζ [x,t]$ \zeta \enspace [{x},t]$. Operational time, on the other hand, is defined through a subordinator, Tα(τ)$ {T}_{\alpha }(\tau )$, a nondecreasing and right continuous function which ensures that its inverse subordinator, Sα(t)$ {S}_{\alpha }(t)$, exists. In essence, as outlined in Raghavan and Chen (2015) if the pressure distribution reflecting the process under consideration were Δp (x,t)$ \Delta p\enspace ({x},t)$, then Δp (x,t)$ \Delta p\enspace ({x},t)$ and Δpc (x,τ)$ \Delta {p}_c\enspace ({x},\tau )$ are related by the integral transform

Δp (x,t)=∫0∞dτ Δpc (x,τ)ψ(τ,t),$$ \Delta p\enspace ({x},t)={\int }_0^{\mathrm{\infty }} \mathrm{d}\tau \enspace \Delta {p}_c\enspace ({x},\tau )\psi (\tau,t), $$(16)

where ψ(τ,t)$ \psi (\tau,t)$ is the probability density of the inverse subordinator, Sα(t)$ {S}_{\alpha }(t)$; it contains all the requirements for the process under consideration (subdiffusion, superdiffusion or a combination). For example, Raghavan and Chen (2015) show that for subdiffusive processes Ψ(s)=sα$ \mathrm{\Psi }(s)={s}^{\alpha }$ where ψ̅(τ,s)$ \overline{\psi }(\tau,s)$ is

ψ̅(τ,s)=Ψ(s)se-τΨ(s).$$ \overline{\psi }(\tau,s)=\frac{\mathrm{\Psi }(s)}{s}{\mathrm{e}}^{-\tau \mathrm{\Psi }(s)}. $$(17)

Fortunately, for the combined case, the needed expressions, ψ̅(τ,s)$ \overline{\psi }\left(\tau,s\right)$ and Δpc (x,τ)$ \Delta {p}_c\enspace \left({x},\tau \right)$ are available (Bazhlekova, 2019; Bazhlekova and Bazhlekov, 2019; Luchko, 2019). Analogous to equation (16) for the combined process, we may show that the Green’s function, φ(β + 1)/2,α$ {\phi }_{({\beta \enspace }+\enspace 1)/2,\alpha }$ (x, t) is

φ(β+1)/2,α(x,t)=∫0∞dτ φ2,1 (x,τ)ψ(τ,t),$$ {\phi }_{(\beta +1)/2,\alpha }({x},t)={\int }_0^{\mathrm{\infty }} \mathrm{d}\tau \enspace {\phi }_{\mathrm{2,1}}\enspace ({x},\tau )\psi (\tau,t), $$(18)

where (Bazhlekova, 2019; Bazhlekova and Bazhlekov, 2019)

ψ̅(τ,s)=sα-1τ(β-1)/2Eα,β(sατ(β+1)/2),$$ \overline{\psi }(\tau,s)={s}^{\alpha -1}{\tau }^{(\beta -1)/2}{E}_{\alpha,\beta }\left({s}^{\alpha }{\tau }^{(\beta +1)/2}\right), $$(19)

with φ2,1 (x,t)$ {\phi }_{\mathrm{2,1}}\enspace ({x},t)$ being the classical expression given by

φ2,1(x,t)=14πηtexp(-x24ηt).$$ {\phi }_{\mathrm{2,1}}({x},t)=\frac{1}{4\pi \eta t}\mathrm{exp}\left(-\frac{{{x}}^2}{4\eta t}\right). $$(20)

As noted in Raghavan and Ozkan (1994) and Raghavan and Chen (2018), Green’s (or source) functions in transient diffusion in porous media reflect the instantaneous change in pressure, d(Δp)/dt, and on using the results of Bazhlekova and Bazhlekov (2019) given by

φ(β+1)/2,α(|x|,t)=12π∫0∞σJ0(|x|σ)Eα(σβ+1tα) dσ,$$ {\phi }_{(\beta +1)/2,\alpha }\left(|{x}|,t\right)=\frac{1}{2\pi }{\int }_0^{\mathrm{\infty }} \sigma {J}_0(|{x}|\sigma ){E}_{\alpha }\left({\sigma }^{\beta +1}{t}^{\alpha }\right)\enspace \mathrm{d}\sigma, $$(21)

we arrive at the expression for the logarithmic-pressure-derivative or pressure derivative, p′(xD,tD)≡td/dt[pD(xD,tD)]$ ({{x}}_{{D}},{t}_D)\equiv t\mathrm{d}/\mathrm{d}t\left[{p}_D\left({{x}}_{{D}},{t}_D\right)\right]$,

tddt[pD(xD,tD)]=tD1α∫0∞σJ0(|xD|σ)Eα(σβ+1tD) dσ.$$ t\frac{\mathrm{d}}{\mathrm{d}t}\left[{p}_D\left({{x}}_{{D}},{t}_D\right)\right]={t}_D^{\frac{1}{\alpha }}{\int }_0^{\mathrm{\infty }} \sigma {J}_0(|{{x}}_{{D}}|\sigma ){E}_{\alpha }\left({\sigma }^{\beta +1}{t}_D\right)\enspace \mathrm{d}\sigma. $$(22)

The symbols pD(xD,tD)$ {p}_D({{x}}_{{D}},{t}_D)$, tD, and xD, reflect the dimensionless pressure, dimensionless time and dimensionless distance, respectively, defined by the following expressions

pD(xD,tD)=2πkα,βhqμη̃1-ααlrefβ+1-2ααΔp(r,t),$$ {p}_D({{x}}_{{D}},{t}_D)=\frac{2\pi {k}_{\alpha,\beta }h}{{q\mu }}\frac{{\stackrel{\tilde }{\eta }}^{\frac{1-\alpha }{\alpha }}}{{\mathcal{l}}_{\mathrm{ref}}^{\frac{\beta +1-2\alpha }{\alpha }}}\Delta p(r,t), $$(23)

tD=η̃lrefβ+1tα,$$ {t}_D=\frac{\stackrel{\tilde }{\eta }}{{\mathcal{l}}_{\mathrm{ref}}^{\beta +1}}{t}^{\alpha }, $$(24)

and

xD=xlref.$$ {{x}}_{{D}}=\frac{{x}}{{\mathcal{l}}_{\mathrm{ref}}}. $$(25)

In the following, as all measurements are expressed with respect to the distance from a line-source well, henceforth, we replace xD by rD, the distance from the center of the well; we take lref to be the well radius, rw.

The analog of equation (22) for the homogeneous case α = 1 and β = 1 is

tddt[pD(xD,tD)]=12exp(-xD24tD).$$ t\frac{\mathrm{d}}{\mathrm{d}t}\left[{p}_D\left({{x}}_{{D}},{t}_D\right)\right]=\frac{1}{2}\mathrm{exp}\left(-\frac{{{x}}_{{D}}^{\mathbf{2}}}{4{t}_D}\right). $$(26)

To our knowledge, the logarithmic pressure derivative was introduced by Chow (1952). But it was not adopted presumably because of the lack of precision in pressure measuring devices that measured pressures suitable for computing pressure derivatives. Procedures such as the maximum slope method of Garcia-Rivera and Raghavan (1979) are an attempt to circumvent this limitation. The introduction of high precision gauges in the late 1970s was a prime factor in the widespread adoption of the pressure derivative technique; see Bourdet et al. (1983). For completeness, we note that other methods to arrive at equation (22) are available in Luchko (2016, 2017) as well as Boyadjiev and Luchko (2017).

For the special situation 2α = β + 1, we have the closed form expression (Luchko, 2016)

tddt[pD(rD,tD)]=12βu1-ααEα,α(-1u),$$ t\frac{\mathrm{d}}{\mathrm{d}t}\left[{p}_D\left({r}_D,{t}_D\right)\right]=\frac{1}{{2}^{\beta }}{u}^{\frac{1-\alpha }{\alpha }}{E}_{\alpha,\alpha }\left(-\frac{1}{u}\right), $$(27)

where

u=4αtDrDβ+1.$$ u=\frac{{4}^{\alpha }{t}_D}{{r}_D^{\beta +1}}. $$(28)

The right-hand side of equation (28) is the grouping we would expect to use based on our knowledge of the Theis (1935) solution and the subdiffusive solutions in Raghavan and Chen (2019); in fact, Gehlhausen (2009) presumes this to be the case a priori in all situations which do not appear to be the case, as we shall see. This expression is important in many ways. First, if α and β were both 1 it reduces to the classsical expression, the expression that follows from the Exponential Integral solution; see equation (26). Second, it suggests that ideas which we intuitively expect from classical diffusion and even subdiffusion should hold in the combined situation for this condition. Third, it provides a method for correlating and understanding transient responses in a general way; basically one may show all results on a single chart if α were used as a parameter of interest much like the Theis (1935) solution.

Most important, however, it is not clear whether solutions may be correlated in terms of tD/rDβ+1$ {t}_D/{r}_D^{\beta +1}$ in general. This matter would have to be determined solely through computations, although we do get a glimpse of what is to be expected through considering the work of Deng and Schilling (2019) who have presented the long-time approximation corresponding to the above expression in equation (22) given by

tddt[pD(rD,tD)] or p'D(rD,tD)=12βΓ(1-β2)Γ(β+12)Γ(1-α)tD1-ααrD1-β; n=2>β+1,$$ t\frac{\mathrm{d}}{\mathrm{d}t}\left[{p}_D\left({r}_D,{t}_D\right)\right]\enspace \mathrm{or}\enspace {{p\prime}_D({r}_D,{t}_D)=\frac{1}{{2}^{\beta }}\frac{\mathrm{\Gamma }\left(\frac{1-\beta }{2}\right)}{\mathrm{\Gamma }\left(\frac{\beta +1}{2}\right)\mathrm{\Gamma }\left(1-\alpha \right)}\frac{{t}_D^{\frac{1-\alpha }{\alpha }}}{{r}_D^{1-\beta }};\hspace{1em}n=2>\beta +1, $$(29)

where n is number of dimensions. Again, interestingly, many observations are evident: this expression suggests that the time dependence of the above expression is essentially similar to that of Raghavan and Chen (2019) for purely subdiffusive flow where p'D(rD,tD→∞)∼tD1-α/αln tD$ {{p}^{\prime}_D\left({r}_D,{t}_D\to \mathrm{\infty }\right)\sim {t}_D^{1-\alpha /\alpha }\mathrm{ln}{\enspace {t}}_D$ (with the influence of the logarithmic term being small) as p′D(rD,tD)$ {p\mathrm{\prime}_D({r}_D,{t}_D)$ is given by

p′D(rD,tD)=121Γ(2-α)tD1-αα{(1-α)[ln4tDe2γrD2-αψ(2-α)]+α},$$ \begin{array}{c}{p\mathrm{\prime}_D\left({r}_D,{t}_D\right)=\frac{1}{2}\frac{1}{\mathrm{\Gamma }\left(2-\alpha \right)}{t}_D^{\frac{1-\alpha }{\alpha }}\\ \left\{(1-\alpha )\left[\mathrm{ln}\frac{4{t}_D}{{e}^{2\gamma }{r}_D^2}-{\alpha \psi }(2-\alpha )\right]+\alpha \right\},\end{array} $$(30)

and where ψ (·) is the Digamma function. Most interesting and surprising, the superdiffusion component plays no role in the exponent of time, t, which is not the case in 1D flow as discussed in Chen and Raghavan (2015); furthermore, the expression makes clear that in the general case p′D(rD,tD)≠f(tD/rDβ+1)$ {p\mathrm{\prime}_D\left({r}_D,{t}_D\right)\ne f({t}_D/{r}_D^{\beta +1})$ unless 2α were equal to β + 1. This issue could become important when we consider more complex geometries for wellbores such as fractured or horizontal wells which are not considered here. Other points of note in this regard are explored below. Note that for α = 1, equation (30) does yield the semi-logarithmic approximation of the Theis (1935) solution, namely the Cooper and Jacob (1946) approximation (also, see Eq. (32)).

In concluding this section, we note that the results in Estrada-Rodriguez et al. (2018) may also be useful and they have derived Green’s functions.

5 Computational results

The principal goal of the computations considered below is to establish the viability of the proposals made to ensure that they are suitable for applications. Second, we establish the validity of the observations made concerning the structure of the solutions; these are particularly revealing in that they illustrate the nature of downhole transients to arrive at reservoir descriptions. Third, we demonstrate some unusual features of superdiffusion which are surprising at least to us. All solutions discussed here were computed through equation (22). Asymptotic behaviors are demonstrated through equation (29) or through other expressions derived in this work. For convenience and also consistency, computations are displayed in terms of tD/rDβ+1$ {t}_D/{r}_D^{\beta +1}$. Pressure responses are computed by straight forward integration of d/dt[pD(rD,tD)].$ \mathrm{d}/\mathrm{d}t\left[{p}_D\left({r}_D,{t}_D\right)\right].$

Figures 1 and 2 present results wherein 2α = β + 1; as indicated by equation (22) we expect solutions to be correlatable in terms of tD/rDβ+1$ {t}_D/{r}_D^{\beta +1}$ which is the case. Figure 1 considers responses in terms of the logarithmic derivative, pD′(rD,tD)$ {p}_D^\mathrm{\prime}\left({r}_D,{t}_D\right)$, for four values of rD [1 (unbroken line), 25 (circles), 50 (squares) and 100 (triangles)]. The results in Figure 1 are for α = 0.9 and β = 0.8. In actuality, we have compared and correlated solutions for rD as large as 250. The dashed line reflects the expression on the right-hand side of equation (29); at long enough times p′D(rD,tD)∝tD(1-α)/α$ {p\mathrm{\prime}_D({r}_D,{t}_D)\propto {t}_D^{(1-\alpha )/\alpha }$. For completeness the derivative response for the Theis (1935) solution is also shown as a dashed line; see equation (26). Not unexpectedly, responses for the classical case are quite distinct from the heterogeneous solutions. Based on this result, we expect the dimensionless pressure, pD(rD,tD), to also correlate as a function of tD/rDβ+1$ {t}_D/{r}_D^{\beta +1}$; see Figure 2. Here, we present both pD(rD,tD) and p′D(rD,tD)$ {p\mathrm{\prime}_D({r}_D,{t}_D)$ and two combinations of α and β (α = 0.9, β = 0.8 and α = 0.8, β = 0.6). The correlation in terms of the circles, triangles, etc. is compatible with the results shown in Figure 1 for α = 0.9, β = 0.8. The bottom two curves (shown as dashes) are the corresponding p′D(rD,tD)$ {p\mathrm{\prime}_D({r}_D,{t}_D)$ responses; at long enough times they, as already noted, will follow the trend predicted by equation (29); the circles denote the start of the straight line, equation (29). Pressure responses at long enough times, pD(rD,tD→∞)$ {p}_D({r}_D,{t}_D\to \mathrm{\infty })$, do also follow the trend tD1-α/α$ {t}_D^{1-\alpha /\alpha }$ provided that α≠1$ \alpha \ne 1$; the long-time expression derived through the Laplace transformation corresponds to:

pD(rD,tD→∞)=12βΓ(1-β2)Γ(β+12)Γ(2-α)tD1-ααrD1-β; n=2>β+1; α≠1,$$ {p}_D\left({r}_D,{t}_D\to \mathrm{\infty }\right)=\frac{1}{{2}^{\beta }}\frac{\mathrm{\Gamma }\left(\frac{1-\beta }{2}\right)}{\mathrm{\Gamma }\left(\frac{\beta +1}{2}\right)\mathrm{\Gamma }\left(2-\alpha \right)}\frac{{t}_D^{\frac{1-\alpha }{\alpha }}}{{r}_D^{1-\beta }};\hspace{1em}n=2>\beta +1;\hspace{1em}\alpha \ne 1, $$(31)

and, incidentally, for purely subdiffusive flows, β = 1, we have

pD(rD,tD)=121Γ(2-α)tD1-αα[ln4tDe2γrD2-αψ(2-α)].$$ {p}_D({r}_D,{t}_D)=\frac{1}{2}\frac{1}{\mathrm{\Gamma }(2-\alpha )}{t}_D^{\frac{1-\alpha }{\alpha }}\left[\mathrm{ln}\frac{4{t}_D}{{e}^{2\gamma }{r}_D^2}-{\alpha \psi }(2-\alpha )\right]. $$(32)

thumbnail Fig. 1

Derivative responses reflecting anomalous diffusion (super- and subdiffusion) under the constraint 2α = β + 1, equation (22). Three values of rD besides the wellbore shown as an unbroken line are considered to demonstrate that pressure responses are a function of tD/rDβ+1$ {t}_D/{r}_D^{\beta +1}$. The dashed line reflects equation (29). For each value of rD, p′D(rD,tD→∞)∼tD(1-α)/α$ {p\mathrm{\prime}_D\left({r}_D,{t}_D\to \mathrm{\infty }\right)\sim {t}_D^{(1-\alpha )/\alpha }$ in agreement with equation (29).

thumbnail Fig. 2

Pressure and derivative distributions reflecting anomalous diffusion (super- and subdiffusion) under the constraint 2α = β + 1, equation (22). Results shown apply for all values of rD. Two sets of α and β values are considered; for one of the cases we specifically show that results may be correlated in terms of tD/rDβ+1$ {t}_D/{r}_D^{\beta +1}$. For each value of rD, both pD (rD, tD → ∞) and p′D(rD,tD→∞)∼tD(1-α)/α$ {p\mathrm{\prime}_D({r}_D,{t}_D\to \mathrm{\infty })\sim {t}_D^{(1-\alpha )/\alpha }$. The topmost dashed line (α = 0.8, β = 0.6) reflects the trend noted by equation (31).

The straight line commensurate to equation (31), however, unlike the derivative responses which begin much earlier, is only evident at times much longer than those considered in Figure 2; and this is often the case. The topmost dashed line shows the pressure trend reflected by equation (31) for α = 0.8, β = 0.6; it is clear that the line will merge with the pressure response.

Figure 3 considers responses where pure superdiffusion is the governing mechanism; that is α = 1. Responses considered are at the wellbore and are quite distinct from those considered so far or those shown in the literature for subdiffusive flow. The curve β = 1 is the line-source solution and pD′(1,tD→∞)→1/2$ {p}_D^\mathrm{\prime}(1,{t}_D\to \mathrm{\infty })\to 1/2$ as expected. For the other cases although equation (31) does not strictly apply if α = 1, by considering α = 1 − ε we expect pD (rD, tD → ∞) ≈ a constant for β < 1; that is, we expect responses to be akin to a constant pressure boundary. This is the case as both p′D(1,tD,β<1)$ {p\mathrm{\prime}_D(1,{t}_D,\beta < 1)$ curves in Figure 3 display a downward trend. As β decreases further the trend to a steady-state-type response will be even more rapid. Recall that all this is occurring in a system that is infinite in its areal extent and where steady flow should not be the norm. In summary, for the situation under discussion, pure superdiffusion in 2D leads to, as alluded earlier a steady type behavior at long enough times. It does not appear to be possible to make reference to any specific characteristic that would by itself permit us to readily identify superdiffusion in the pressure trace which would lead to the determination of β in a manner similar through equation (29). Note that in equation (32) the logarithmic nature of the pressure response is preserved for α = 1 whereas in the case of equation (31) no conclusion may be drawn. Presumably because of the lack of an expression similar to equation (31) neither Cloot and Botha (2006), nor Gehlhausen (2009) recognize the possibility of steady flow. The above observations also apply to rD > 1. The slopes of the downward sections of the p′D(rD,tD)$ {p\mathrm{\prime}_D({r}_D,{t}_D)$ derivative curves are independent of rD and depend only on β. As no analytical options are available at this time, one option is to determine β empirically through the exponents after determining the slope of the downward section; see Figure 4 which may be represented by

m=0.5006β3-1.470β2+1.992β-1.021; β<1, α=1.$$ m=0.5006{\beta }^3-1.470{\beta }^2+1.992\beta -1.021;\hspace{1em}\beta < 1,\hspace{1em}\alpha =1. $$(33)

thumbnail Fig. 3

Pressure and derivative responses at the wellbore, rD = 1, reflecting the influence of pure anomalous, superdiffusion; that is, responses obtained under the presumption that α = 1. The downward trend in the derivative curves is to be expected as the pressure response at the well approaches a steady response as the predicted by equations (22) and (29).

thumbnail Fig. 4

Exponents of the downward sections of the derivative curves for superdiffusion to determine β.

Simply put, in essence, pure superdiffusion leads to negative slopes at long times.

Because of the unusual aspects seen in Figure 3, Figure 5 explores the characteristic behaviors near the classical-subdiffusive limit by considering the role of subdiffusion in the range 0.9≤α<1$ 0.9\le \alpha < 1$ and we consider p′D$ {p\mathrm{\prime}_D$ for four values of α at the wellbore in Figure 5 (0:9000; 0:9250; 0:9900 and 0:9999). Pressure responses are shown for only α of 0.9000. Straight lines corresponding to the asymptotic solution, namely equation (31), are evident for α≤0.9250$ \alpha \le 0.9250$. For larger values no straight line is evident for the time ranges considered here.

thumbnail Fig. 5

Pressure and derivative responses at the wellbore, rD = 1, reflecting influence of subdifussion near the classical limit of α of 1 under the presumption that both super- and subdiffusion effects govern transients. Downward trend in the derivative curves for α close to the limit of classical diffusion suggests that the well pressure approaches a steady value.

Returning to the results in Figure 3, we call attention to the observation not explicitly detailed earlier that in linear systems the pressure does depend on the long-range exponent, β, for α = 1, for Chen and Raghavan (2015) show that

pD (xD=0,tD)=π2tDβ+1-αα(β+1)Γ(2-αβ+1).$$ {p}_D\enspace ({x}_D=0,{t}_D)=\frac{\pi }{2}\frac{{t}_D^{\frac{\beta +1-\alpha }{\alpha (\beta +1)}}}{\mathrm{\Gamma }\left(2-\frac{\alpha }{\beta +1}\right)}. $$(34)

Furthermore, one important insight that should be garnered on comparing equation (34) with equation (27) is that for 1D flow the response for the neutral case, 2α = β + 1, equation (34) yields the homogeneous solution; see Magin et al. (2013) and Chen and Raghavan (2015) and this is not the case in 2D. Incidentally, the works by Chu and coworkers tellingly illustrate advantages that accrue through the use of equation (34).

Figure 6 considers the general case of combined diffusion for the situation 2α≠β+1$ 2\alpha \ne \beta +1$; nevertheless pressure distributions are plotted in terms of tD/rDβ+1$ {t}_D/{r}_D^{\beta +1}$ for ease of comparison. Responses are typical of the myriad solutions considered by us. It is clear that in general solutions depend on both α and β; the general shapes are no different from those considered earlier. Indeed at long enough times the exponent of t is independent of β, and equation (29) does remain available. The lines are not shown here as the onset of the straight lines happens to be much later than the times considered here. Long-time responses, however, are correlatable in terms of tD/rD1-β$ {t}_D/{r}_D^{1-\beta }$; see equation (27).

thumbnail Fig. 6

Pressure derivative distributions reflecting anomalous diffusion (super- and subdiffusion) for the general case. The constraint, 2α = β + 1, does not apply. Results are presented in terms of tD/rDβ+1$ {t}_D/{r}_D^{\beta +1}$ principally for convenience.

5 Discussion and concluding remarks

The principal goal of this work is to incorporate the influences of macro-and microscale heterogeneities at the Theis scale around a line-source well. Anomalous processes of fractional calculus are used to evaluate pressure responses governing space-time fractional diffusion. Two dimensional super- and subdiffusion solutions are, to our knowledge, yet to be explored. One immediate advantage of this work is that it provides an understanding of long term behaviors of the flow of fluids produced through complex wellbores (fractured wells, horizontal wells), for pseudoradial flow will ultimately prevail in all situations. For example, if a well produces through a finite-conductivity fracture, then exponents of (β+1-α)/[2(β+1)],(β+1-α)/(β+1)$ (\beta +1-\alpha )/[2(\beta +1)],(\beta +1-\alpha )/(\beta +1)$ and 1-α corresponding to the bilinear, linear and pseudoradial flow-periods, respectively, should be evident on the derivative trace whenever all three flow regimes are evident; see Raghavan and Chen (2020). Responses to a line-source well are presented both analytically and graphically. Salient features of the results obtained are as follows: (i) characteristic features of super- and subdiffusion and pure subdiffusion at long enough times are identical, (ii) derivative curves for pure superdiffusion yield negative slopes, and (iii) pressure distributions fall under two classifications, those that meet the restriction 2α = β + 1 in which case pD=f[tD/rDβ+1]$ {p}_D=f[{t}_D/{r}_D^{\beta +1}]$, the grouping expected from the Theis solution, and the general case wherein 2α≠β+1$ 2\alpha \ne \beta +1$; regardless, at long enough times pD∼tD1-α/α$ {p}_D\sim {t}_D^{1-\alpha /\alpha }$ provided that α,β≠1$ \alpha,\beta \ne 1$. Computational results presented here may be used to address many applications in the earth sciences: movement of groundwater, production of hydrocarbon and geothermal resources, sequestration etc. Studies of this type provide a richer understanding of transient diffusion in heterogeneous porous solids and may be extended to other similar applications.

Acknowledgments

We thank Dr. Wei Chun Chu for reviewing the typescript and his comments.

References

All Figures

thumbnail Fig. 1

Derivative responses reflecting anomalous diffusion (super- and subdiffusion) under the constraint 2α = β + 1, equation (22). Three values of rD besides the wellbore shown as an unbroken line are considered to demonstrate that pressure responses are a function of tD/rDβ+1$ {t}_D/{r}_D^{\beta +1}$. The dashed line reflects equation (29). For each value of rD, p′D(rD,tD→∞)∼tD(1-α)/α$ {p\mathrm{\prime}_D\left({r}_D,{t}_D\to \mathrm{\infty }\right)\sim {t}_D^{(1-\alpha )/\alpha }$ in agreement with equation (29).
In the text
thumbnail Fig. 2

Pressure and derivative distributions reflecting anomalous diffusion (super- and subdiffusion) under the constraint 2α = β + 1, equation (22). Results shown apply for all values of rD. Two sets of α and β values are considered; for one of the cases we specifically show that results may be correlated in terms of tD/rDβ+1$ {t}_D/{r}_D^{\beta +1}$. For each value of rD, both pD (rD, tD → ∞) and p′D(rD,tD→∞)∼tD(1-α)/α$ {p\mathrm{\prime}_D({r}_D,{t}_D\to \mathrm{\infty })\sim {t}_D^{(1-\alpha )/\alpha }$. The topmost dashed line (α = 0.8, β = 0.6) reflects the trend noted by equation (31).
In the text
thumbnail Fig. 3

Pressure and derivative responses at the wellbore, rD = 1, reflecting the influence of pure anomalous, superdiffusion; that is, responses obtained under the presumption that α = 1. The downward trend in the derivative curves is to be expected as the pressure response at the well approaches a steady response as the predicted by equations (22) and (29).
In the text
thumbnail Fig. 4

Exponents of the downward sections of the derivative curves for superdiffusion to determine β.
In the text
thumbnail Fig. 5

Pressure and derivative responses at the wellbore, rD = 1, reflecting influence of subdifussion near the classical limit of α of 1 under the presumption that both super- and subdiffusion effects govern transients. Downward trend in the derivative curves for α close to the limit of classical diffusion suggests that the well pressure approaches a steady value.
In the text
thumbnail Fig. 6

Pressure derivative distributions reflecting anomalous diffusion (super- and subdiffusion) for the general case. The constraint, 2α = β + 1, does not apply. Results are presented in terms of tD/rDβ+1$ {t}_D/{r}_D^{\beta +1}$ principally for convenience.
In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.