Analytical Solutions of the Diffusion–Wave Equation of Groundwater Flow with Distributed-Order of Atangana–Baleanu Fractional Derivative

Abstract

A generalized mathematical model of the radial groundwater flow to or from a well is studied using the time-fractional derivative with Mittag-Lefler kernel. Two temporal orders of fractional derivatives which characterize small and large pores are considered in the fractional diffusion–wave equation. New analytical solutions to the distributed-order fractional diffusion–wave equation are determined using the Laplace and Dirichlet-Weber integral transforms. The influence of the fractional parameters on the radial groundwater flow is analyzed by numerical calculations and graphical illustrations are obtained with the software Mathcad.

1. Introduction

Groundwater, in nearly all geological formations, is visible and in movement. Movement is very sluggish in some materials such as rock or clay, as opposed to the groundwater pouring in sand and gravel patterns. However, highly eroded rock may display discharges of groundwater identical to or in surplus of sand and gravel aquifers [1]. In almost every part of the globe, aquifers are essential mineral resources for irrigated agriculture, commercial and drinkable water supply. The study of aquifers is germane to chemical engineers, agricultural engineers, civil engineers, forestry engineers, ecologists, geologists, geographers, geotechnical engineers, irrigation engineers, mining engineers, hydraulic engineers, hydro-geologists, hydrologists, petroleum engineers, soil scientists, and others [2].

In groundwater research, the central issue faced so far is the actual configuration of the geological structure in which water flows into the aquifer under analysis. However, there are many fractured rock aquifers where groundwater flow does not conform to traditional geometries [3], such as the Karoo aquifers in South Africa, described by the existence of a few parallel bedding fractures that act as the primary water channels in aquifers [4]. For the Karoo aquifer, Botha et al. [4] established a three-dimensional prototype and demonstrated that the dominant flow field in these aquifers is vertical and linear and not, as generally thought, horizontal and radial. However, subsequent investigations [5] indicate that the flow is often affected by the geometry of bedding parallel fractures. It was concluded that the most appropriate modification is to generalize the flow dimension to non-integer values after considering a number of potential variations on the models, while maintaining the premises of radial flow and homogeneity [6]. Analogously, in order to investigate a radially symmetric model for the groundwater flow, Cloot and Botha suggested the idea of a non-integer fractional derivative by changing the conventional first order derivative of a piezometric head with a complementary fractional derivative [7]. To investigate groundwater flow, Cloot and Botha [7] have considered the Riemann–Liouville fractional derivatives ([8], chapter 1, Equation (1.25)) to formulate the fractional diffusion–wave equations for the radial solute transport in aquifers, and this formulation is similar to the equation suggested by Metzler et al. [9].

During aquifer analyses in fractured networks, Cello et al. [10] analyzed the impact of flow dimensions and irregular diffusion of groundwater. Quite recently, Atangana and Bildik [11] reused a time-fractional partial differential equation (PDE) proposed by Metzler et al. [9] and three variants of fractional partial differential equations (fPDEs) in a radial coordinate that were considered by Atangana and Vermeulen [12], Atangana [13]. More recently, the study of groundwater flows in aquifers using fPDEs for groundwater transport in unsaturated soils has been addressed [14,15,16,17,18].

In several scientific fields, such as geoscience, demography, physics, mathematics, biology, and bioengineering, the simulation of complex systems is made with fractional order differential/integral operators [19]. Many fractional–differential operators have been studied by researchers. We recall here some of them: namely, the fractional Caputo–Fabrizio operator [20], the integral/derivative Riemann–Liouville operator [21], fractional Caputo derivative [22], the fractional Yang–Srivastava–Machado derivative [23], the Atangana–Baleanu fractional derivative [24,25], etc. The space–fractional diffusion with power–law super diffusivity formulated with the Riemann–Liouville fractional derivative was analyzed by Hristov [26]. The time–fractional Caputo–Fabrizio derivative was implemented by Ahmad et al. [27] to examine the advective diffusion process with low-range memory and concentrated source.

The fractional derivative with exponential and non-singular kernel defined by Caputo and Fabrizio is largely used in the modeling of dissipative phenomena by fractional derivatives. The Caputo–Fabrizio integral operator with exponential memory (so called fractional derivative with a non-singular kernel) is oriented to more complex relaxation processes than the ones modeled by the classical fractional derivatives with power–law (fractional derivatives with singular kernel) memories. Applications of the Caputo–Fabrizio integral operator with memory in the modeling of the dissipative phenomena with exponentially decaying relaxation functions have been discussed by Hrisotv [28].

It is known that the most popular function fitting the time–domain relaxation patterns is the Kohlraush stretched-exponential function [29] $\phi (k,t)=\exp (-{(at)}^{{\alpha }_k})$ where ${\alpha }_k, 0\le {\alpha }_k\le 1$ is the stretching exponent and $a$ is an inverse of the characteristic material time constant. Barbera-Santos et al. [30] showed that the Kohlraush relaxation function is a kernel derived from experiments.

It is easy to notice that the Caputo–Fabrizio derivative operator has been constructed with a simple exponential memory; namely, the Kohlraush stretched exponent is set to one. In a similar way, a generalization of the Caputo–Fabrizio derivative was constructed by considering a kernel that generalizes the exponential function, namely the Mittag–Leffler function with a single parameter. The new operator is called the Atangana–Baleanu fractional derivative.

The choosing of the memory kernel and the recovery of the relaxation function from the experimental data and fitting it with a well-defined function is a task strongly dependent on the physics of the process of interest. The memory kernels in the fractional differential operators used to model linear viscoelasticity correspond to the responses of the material, namely to the relaxation moduli and compliances [28]. It is important to note that Mittag–Leffler functions are solutions of some fractional differential equations; therefore, they are non-local functions. In consequence, the time–fractional Atangana–Baleanu derivative is a fractional derivative with non-singular and non-local kernel. It is expected that the non-locality of the Mittag–Leffler kernel allows better description of the memory within the structure and media with different scales. Also, it is known that Mittag–Leffler functions are suitable for interesting generalizations. Thus, generalizations of these functions with three or four parameters were studied. Considered to be the kernel of fractional derivatives, these functions become useful in describing some complex physical processes. For example, the three-parameter Mittag–Leffler functions studied by the Indian mathematician Prabhakar play an important role in describing anomalous dielectric relaxation, stochastic processes, or renewal processes [31].

Long et al. [32] studied four fractional viscoelastic models, namely the fractional Maxwell model, fractional Kelvin–Voigt model, fractional Zener model and fractional Poynting–Thomson model using the Riemann–Liouvile and Caputo derivatives, as well as the Caputo–Fabrizio and Atangana–Baleanu derivatives. For each fractional viscoelastic model, the stress relaxation modulus, creep compliance and dynamic modulus have been derived and compared. Their results show that the fractional Maxwell model and fractional Zener model with the Mittag–Leffler function kernel do not provide accurate descriptions of the stress relaxation modulus at the shortest time and the storage modulus at the highest frequency. However, these results are not a verdict against the derivatives with non-singular and non-local kernel that could describe different complex processes.

It is important to emphasize the different assessments of researchers in the field of fractional calculus and their applications, regarding the criteria according to which such an operator is really a fractional operator or is the equivalent of an operator of integer order. Most of the discussions were clarified by Creson and Szafranska, [33], respectively, Tarasov and Tarasova [34] formulated criteria for this selection.

In [35], the authors have investigated the distributed-order fractional diffusion–wave equation based on the Caputo fractional derivative in time, and determined solutions for the fluid flow to or from a well. Two temporal orders of fractional derivatives to describe large and small pores, respectively, are incorporated in the distributed-order fractional diffusion–wave equation.

In this paper we investigate a generalized mathematical model of the radial groundwater flow to or from a well using the time–fractional derivative with Mittag–Lefler kernel defined by Atangana and Baleanu. Following Su et al. [35], we consider the fractional diffusion–wave equation with two temporal orders of fractional derivatives to characterize the diffusion process.

Using the Laplace and Dirichlet–Weber transform, the analytical solutions of the distributed-order fractional diffusion–wave equation are determined. The obtained solutions are new in the literatures and generate the solutions corresponding to ordinary process for fractional parameters equal to one. The influence of the fractional parameters on the radial groundwater flow is analyzed by numerical calculations and graphical illustrations are obtained with the software Mathcad.

2. The Formulation of the Generalized Problem

The flow to or from a well is described by the following diffusion equation [1,36]

$$\frac{\partial H}{\partial t}=\frac{1}{r}\frac{\partial }{\partial r}\left(Dr\frac{\partial H}{\partial r}\right),$$

(1)

which can be written equivalently as

$$C_1\frac{\partial H}{\partial t}=\left(\frac{{\partial }^{2}H}{\partial r^2}+\frac{1}{r}\frac{\partial H}{\partial r}\right),$$

(2)

where $H$ being the normalized groundwater depth and is given by

$$H=\frac{h-h_0}{h_m-h_0},$$

(3)

with $h$ is the groundwater height above the datum, $h_0$ is the original height of groundwater above the datum before the test, $h_m$ is the extreme groundwater height following the injection or pumping, $D$ is the aquifer diffusivity, and $r$ is the radial distance from the well center. $C_1=\frac{1}{D}$ or $D=\frac{1}{C_1}$ for both unconfined and restricted aquifers described by Bras [36] as:

(1)

for an unconfined aquifer, the diffusivity $D$ is

$$D=\frac{1}{C_1}=\frac{T}{n_e},$$

(4)

(2)

for a confined aquifer, the diffusivity $D$ is

$$D=\frac{1}{C_1}=\frac{T}{S},$$

(5)

where $T=K{h}_m$ is called the transmissivity with $K$ being the hydraulic conductivity, $n_e$ being the effective porosity in an unconfined aquifer and $S$ being the specific storage coefficient [37,38].

Equation (4) can also be written

$$D=K_e{h}_m$$

(6)

with $K_e$ being called the effective hydraulic conductivity given as

$$K_e=\frac{K}{n_e}$$

(7)

where $T=K{h}_m$, which occurs in Equation (4) suggests that Equation (1) is a linearization of Boussinesq’s radial coordinate equation [39]. Equation (6) results when the height perturbation of injected water or pumped drawdown is small enough to assume that $H\approx h_m$.

To take into account the impact of pores of various sizes on saturated flow as a complement of the distributed-order fractional partial differential equation for water flow in an unsaturated medium [18], we incorporate distributed-order in Equation (1) to construct the distributed-order fractional partial differential equation of diffusion.

We rewrite Equation (1) with the Atangana–Balean fractional derivative in time with the temporal distributed orders in the following form,

$$a^{ABC}{\mathcal{D}}_t^{\alpha }H(r,t)+b^{ABC}{\mathcal{D}}_t^{\beta }H(r,t)=\frac{D_f}{r}\frac{\partial }{\partial r}\left(r\frac{\partial H(r,t)}{\partial r}\right)$$

(8)

where: $a$ and $b$ are, respectively, the ratio of immobile and mobile section porosities with the total porosity and satisfy the relation $a+b=1$; $0<\alpha \le 1$ and $0<\beta \le 1$ are, respectively, the fractional derivative parameters for the mobile and the immobile sections; and $D_f$ is the fractional diffusivity, which also takes two forms for unconfined and confined aquifers as their traditional complements in Equations (4) and (5). As this can be seen in Equations (4) and (5) for two significant groundwater kinds, i.e., unconfined and confined aquifers, Df has specific meaning, but the major component of Equation (8) remains mathematically the same. An alternative approach to describe the differentiated drainage from the two main forms of pores is the use of the distributed orders in Equation (8). Of course, further orders to account for main types of pores can be put in temporal–fractional derivatives, but the determination of their values will become a technical issue.

The Atangana–Baleanu fractional derivative in the Caputo sense is defined as [24,25,40]

$${}^{ABC}{\mathcal{D}}_t^{\beta }g(r,t)=\{\begin{array}{l}\frac{N(\beta )}{1-\beta }{\displaystyle {\int }_0^t\dot{g}(r,\tau )}{E}_{\beta }(\frac{-\beta }{1-\beta }{(t-\tau )}^{\beta })d\tau , 0<\beta <1,\\ \frac{\partial g(r,t)}{\partial t}, \beta =1.\end{array}$$

(9)

where $N(\beta )$ is the normalization function such that $N(0)=N(1)=1$, $\dot{g}(r,\tau )={\frac{\partial g(r,t)}{\partial t}|}_{t=\tau }$ and $E_{\beta }(·)$ is the Mettag–Leffler function [41]. In this paper we consider $N(\beta )=1$.

If $g(r,0)=0$, the Laplace transform of the Atangana–Baleanu fractional derivative ${}^{ABC}{\mathcal{D}}_t^{\beta }g(r,t)$ is given by

$$L{\{}^{ABC}{\mathcal{D}}_t^{\beta }g(r,t);p\}=\frac{p^{\beta }}{(1-\beta )p^{\beta }+\beta }\cdot \overline{g}(r,p),$$

(10)

where $\overline{g}(r,p)={\displaystyle \underset{0}{\overset{\infty }{\int }}g(r,t)e^{-pt}dt}$ denotes the Laplace transform of $g(r,t)$.

3. Solution to the Problem

In this section, the analytical solution of the fractional differential Equation (8) will be determined along with the initial and boundary conditions

$$H(r,0)=0,\hspace{1em}r\in [R,\infty ),$$

(11)

$$\frac{\partial H(r,t)}{\partial r}{|}_{r=R}=-\frac{Q}{2\pi R{T}_f}$$

(12)

$$\sqrt{r}\frac{\partial H(r,t)}{\partial r}\to 0, \sqrt{r}\frac{{\partial }^{2}H(r,t)}{\partial r^2}\to 0 as r\to \infty .$$

(13)

Applying the Laplace transform to Equations (8), (12) and (13), and by using the initial condition (11), we get the transformed problem

$$\frac{1}{r}\frac{\partial }{\partial r}(r\frac{\partial \overline{H}(r,s)}{\partial r})=\frac{1}{D_f}[\frac{a{s}^{\alpha }}{(1-\alpha )s^{\alpha }+\alpha }+\frac{b{s}^{\beta }}{(1-\beta )s^{\beta }+\beta }]\overline{H}(r,s),\hspace{1em}\alpha ,\beta \in (0,1],$$

(14)

$$\frac{\partial \overline{H}(r,s)}{\partial r}{|}_{r=R}=-\frac{Q}{2\pi R{T}_f}\frac{1}{s}$$

(15)

$$\sqrt{r}\frac{\partial \overline{H}(r,s)}{\partial r}\to 0,\hspace{1em}\sqrt{r}\frac{{\partial }^2\overline{H}(r,s)}{\partial r^2}\to 0, \mathrm{as} r\to \infty$$

(16)

In order to find the solution to Equations (14)–(16) we use the Dirichlet–Weber transform defined as [42]

$$\widehat{\overline{H}}(\xi ,s)={\displaystyle {\int }_R^{\infty }r}\overline{H}(r,s)[J_{\nu }(r\xi )Y_{\nu }(R\xi )-J_{\nu }(R\xi )Y_{\nu }(r\xi )]dr,$$

(17)

where $J_{\nu }(\cdot )$, $Y_{\nu }(\cdot )$ are Bessel functions of the first kind with order $\nu$.

The inverse formula of the integral transform (17) is [42]

$$\overline{H}(r,s)={\displaystyle {\int }_0^{\infty }\frac{J_{\nu }(r\xi )Y_{\nu }(R\xi )-J_{\nu }(R\xi )Y_{\nu }(r\xi )}{J_{\nu }^2(R\xi )+Y_{\nu }^2(R\xi )}}\xi ·\widehat{\overline{H}}(\xi ,s)d\xi .$$

(18)

Using the notation,

$$\overline{{\phi }_0}(s)=\frac{1}{D_f}[\frac{a{s}^{\alpha }}{(1-\alpha )s^{\alpha }+\alpha }+\frac{b{s}^{\beta }}{(1-\beta )s^{\beta }+\beta }]$$

(19)

into Equation (14), we obtain the transformed equation,

$$\frac{{\partial }^2\overline{H}(r,s)}{\partial r^2}+\frac{1}{r}\frac{\partial \overline{H}(r,s)}{\partial r}={\overline{\phi }}_0(s)\overline{H}(r,s)$$

(20)

which can be written in the equivalent form

$$\frac{{\partial }^2}{\partial r^2}\frac{\partial \overline{H}}{\partial r}+\frac{1}{r}\frac{\partial }{\partial r}\frac{\partial \overline{H}}{\partial r}-\frac{1}{r^2}\frac{\partial \overline{H}}{\partial r}={\overline{\phi }}_0(s)\frac{\partial \overline{H}}{\partial r}$$

(21)

Making the change of unknown function $\overline{H}(r,s)$ into $\overline{u}(r,s)=\frac{\partial \overline{H}(r,s)}{\partial r}$, Equation (21) becomes

$$\frac{{\partial }^2\overline{u}}{\partial r^2}+\frac{1}{r}\frac{\partial \overline{u}}{\partial r}-\frac{1}{r^2}\overline{u}={\overline{\phi }}_0(s)\overline{u}$$

(22)

The function $\overline{u}(r,s)$ has to satisfy the boundary condition

$$\overline{u}(R,s)=-\frac{Q}{2\pi R{T}_f}\frac{1}{s}$$

(23)

$$\sqrt{r}\frac{\partial \overline{u}}{\partial r}\to 0,\sqrt{r}\overline{u}(r,s)\to 0,\hspace{1em}\mathrm{as}r\to \infty$$

(24)

Using the asymptotic expansions

$$J_{\alpha }(z)\cong \sqrt{\frac{2}{\pi z}}\cos (z-\frac{(2\alpha +1)\pi }{4});\hspace{1em}{Y}_{\alpha }(z)\cong \sqrt{\frac{2}{\pi z}}\sin (z-\frac{(2\alpha +1)\pi }{4}), z\to \infty$$

and the identity

$$J_0(R\xi )Y_1(R\xi )-J_1(R\xi )Y_0(R\xi )=-\frac{2}{\pi R\xi }$$

it easily turns out that

$$I={\displaystyle {\int }_R^{\infty }r}[\frac{{\partial }^2\overline{u}}{\partial r^2}+\frac{1}{r}\frac{\partial \overline{u}}{\partial r}-\frac{\overline{u}}{r^2}][J_1(r\xi )Y_1(R\xi )-J_1(R\xi )Y_1(r\xi )]dr=-{\xi }^2\widehat{\overline{u}}(\xi ,s)-\frac{2}{\pi }\overline{u}(R,s)$$

(25)

Now, applying the Dirichlet–Weber transform to Equation (22), and using Equation (25) and the boundary condition (23), we obtain the transformed equation

$$\widehat{\overline{u}}(\xi ,s)=\frac{k_0}{s({\xi }^2+{\overline{\phi }}_0(s))}$$

(26)

where $k_0=\frac{Q}{{\pi }^{2}R{T}_f}$. Using the function $f(r)=\frac{R}{r}$ with the Dirichlet–Weber transform $\widehat{f}(\xi )=-\frac{2}{\pi {\xi }^2}$, we write the function of Equation (26) in the suitable form

$$\widehat{\overline{u}}(\xi ,s)=\frac{k_0}{{\xi }^{2}s}-\frac{k_0{\overline{\phi }}_0(s)}{{\xi }^{2}s({\xi }^2+{\overline{\phi }}_0(s))}=\left(\frac{-Q}{2\pi R{T}_f}\right)\left(\frac{-2}{\pi {\xi }^{2}s}\right)-\frac{k_0{\overline{\phi }}_0(s)}{{\xi }^{2}s({\xi }^2+{\overline{\phi }}_0(s))}$$

(27)

We introduce the notation

$$\widehat{\overline{U}}(\xi ,s)=-\frac{k_0{\overline{\phi }}_0(s)}{{\xi }^{2}s({\xi }^2+{\overline{\phi }}_0(s))}.$$

(28)

The expression of function $\widehat{\overline{U}}(\xi ,s)$ is written as

$$\widehat{\overline{U}}(\xi ,s)=\frac{-k_0}{{\xi }^{2}s}\frac{p_0{s}^{\alpha +\beta }+p_1{s}^{\alpha }+p_2{s}^{\beta }}{q_0(\xi )s^{\alpha +\beta }+q_1(\xi )s^{\alpha }+q_2(\xi )s^{\beta }+q_3(\xi )}.$$

(29)

where

$$p_0=a(1-\beta )+b(1-\alpha ), p_1=a\beta , p_2=b\alpha ,$$

(30)

$$\begin{array}{c}{q}_0(\xi )=D_f(1-\alpha )(1-\beta ){\xi }^2+a(1-\beta )+b(1-\alpha )\\ q_1(\xi )=D_f(1-\alpha )\beta {\xi }^2+a\beta \\ q_2(\xi )=D_f(1-\beta )\alpha {\xi }^2+b\alpha \\ q_3(\xi )=D_f\alpha \beta {\xi }^2\end{array}$$

(31)

We consider here three cases.

3.1. The Case $\alpha \le \beta$

$$\widehat{\overline{U}}(\xi ,s)=\frac{-k_0}{{\xi }^2}{\displaystyle \sum _{k=0}^{\infty }{\displaystyle \sum _{m=0}^k\frac{\frac{{(-1)}^k}{m!}\frac{k!}{(k-m)!}\frac{q_1^m(\xi )q_3^{k-m}(\xi )}{q_0^{k+1}(\xi )}}{{(s^{\alpha }+\frac{q_2}{q_0})}^{k+1}}\left[p_0{s}^{\alpha -1+m\alpha -k\beta }+p_1{s}^{\alpha -\beta -1+m\alpha -k\beta }+p_2{s}^{-1+m\alpha -k\beta }\right]}}$$

(32)

Using the generalized G-functions of Lorenzo–Hartley [43], defined as

$$G_{a,b,c}(t,d)={\mathfrak{L}}^{-1}\left\{\frac{s^b}{{(s^a-d)}^c}\right\},\hspace{1em}Re(ac-b)>0,$$

(33)

we obtain

$$\widehat{\overline{U}}(\xi ,t)=\frac{-k_0}{{\xi }^2}{\displaystyle \sum _{k=0}^{\infty }{\displaystyle \sum _{m=0}^k\frac{{(-1)}^k}{m!}\frac{k!}{(k-m)!}\frac{q_1^m(\xi )q_3^{k-m}(\xi )}{q_0^{k+1}(\xi )}\left[\begin{array}{l}{p}_0{G}_{\alpha ,\alpha (m+1)-k\beta -1,k+1}\left(t,-\frac{q_2(\xi )}{q_0(\xi )}\right)+p_1{G}_{\alpha ,\alpha (m+1)-k\beta -1-\beta ,k+1}\left(t,-\frac{q_2(\xi )}{q_0(\xi )}\right)\\ +p_2{G}_{\alpha ,\alpha m-k\beta -1,k+1}\left(t,-\frac{q_2(\xi )}{q_0(\xi )}\right)\end{array}\right]}}$$

(34)

3.2. The Case $\alpha =\beta =1$ (The Ordinary Case)

In this case,

$$\begin{array}{l}{p}_0=0, p_1=a, p_2=b, \\ q_0(\xi )=0, q_1(\xi )=a, q_2(\xi )=b, q_3(\xi )=D_f{\xi }^2\end{array}$$

(35)

$$\widehat{\overline{U}}(\xi ,s)=\frac{-k_0}{{\xi }^2}\frac{1}{s+D_f{\xi }^2},$$

(36)

with the inverse Laplace transform

$$\widehat{U}(\xi ,t)=-\frac{k_0}{{\xi }^2}{e}^{-D_f{\xi }^{2}t}$$

(37)

3.3. The Case $\beta \le \alpha$

Equation (32) can be written in following equivalent form

$$\widehat{\overline{U}}(\xi ,s)=\frac{-k_0}{{\xi }^2}{\displaystyle \sum _{k=0}^{\infty }{\displaystyle \sum _{m=0}^k\frac{\frac{{(-1)}^k}{m!}\frac{k!}{(k-m)!}\frac{q_2^m(\xi )q_3^{k-m}(\xi )}{q_0^{k+1}(\xi )}}{{(s^{\beta }+\frac{q_1}{q_0})}^{k+1}}\left[p_0{s}^{(m+1)\beta -1-k\alpha }+p_1{s}^{m\beta -k\alpha -1}+p_2{s}^{(m+1)\beta -(k+1)\alpha -1}\right]}}$$

(38)

whose inverse Laplace transform is

$$\widehat{U}(\xi ,t)=\frac{-k_0}{{\xi }^2}{\displaystyle \sum _{k=0}^{\infty }{\displaystyle \sum _{m=0}^k\frac{{(-1)}^k}{m!}\frac{k!}{(k-m)!}\frac{q_2^m(\xi )q_3^{k-m}(\xi )}{q_0^{k+1}(\xi )}\left[\begin{array}{l}{p}_0{G}_{\beta ,(m+1)\beta -k\alpha -1,k+1}\left(t,-\frac{q_1(\xi )}{q_0(\xi )}\right)+p_1{G}_{\beta ,\beta m-k\alpha -1,k+1}\left(t,-\frac{q_1(\xi )}{q_0(\xi )}\right)+\\ p_2{G}_{\beta ,\beta (m+1)-\alpha (k+1)-1,k+1}\left(t,-\frac{q_1(\xi )}{q_0(\xi )}\right)\end{array}\right]}}$$

(39)

Applying the inverse Dirichlet–Weber transform to Equation (27), and using the Dirichlet–Weber transform of the function $f(r)=\frac{1}{r},\hspace{1em}r\in [R,\infty )$, namely,

$$\begin{array}{l}\widehat{f}(\xi )={\displaystyle {\int }_R^{\infty }\frac{1}{r}}r[J_1(r\xi )Y_1(R\xi )-J_1(R\xi )Y_1(r\xi )]dr=\frac{Y_1(R\xi )}{\xi }{\displaystyle {\int }_{R\xi }^{\infty }{J}_1}(z)dz-\frac{J_1(R\xi )}{\xi }{\displaystyle {\int }_{R\xi }^{\infty }{Y}_1}(z)dz=\\ \frac{Y_1(R\xi )}{\xi }[-J_0(z)]{|}_{{}_{R\xi }}^{{}^{\infty }}-\frac{J_1(R\xi )}{\xi }[-Y_0(z)]{|}_{{}_{R\xi }}^{{}^{\infty }}=\frac{1}{\xi }[J_0(R\xi )Y_1(R\xi )-J_1(R\xi )Y_0(R\xi )]=-\frac{2}{\pi R{\xi }^2},\end{array}$$

(40)

with the inverse Dirichlet–Weber transform,

$$f(r)=\frac{R}{r}=\{\begin{array}{cc}1,\hfill & r=R\hfill \\ {\displaystyle \underset{0}{\overset{\infty }{\int }}\frac{J_1(r\xi )Y_1(R\xi )-J_1(R\xi )Y_1(r\xi )}{J_1^2(R\xi )+Y_1^2(R\xi )}R\xi \left(-\frac{2}{\pi R{\xi }^2}\right)d\xi ,}& r>R\hfill \end{array}$$

(41)

we obtain

$$u(r,t)=\left(\frac{-Q}{2\pi T_f}\right)\frac{1}{r}+{\displaystyle {\int }_0^{\infty }\frac{J_1(r\xi )Y_1(R\xi )-J_1(R\xi )Y_1(r\xi )}{J_1^2(R\xi )+Y_1^2(R\xi )}}\xi ·\widehat{U}(\xi ,t)d\xi$$

(42)

where $\widehat{U}(\xi ,t)$ is given by Equation (34) for $\alpha ,\beta \in (0,1), \alpha \le \beta ,$ respectively $\widehat{U}(\xi ,t)$ is given by Equation (37) for $\alpha =\beta =1$ and $\widehat{U}(\xi ,t)$ is given by Equation (39) for $\alpha ,\beta \in (0,1), \beta \le \alpha .$

Using the transformation $u(r,t)=\frac{\partial H(r,t)}{\partial r}$ and the following properties of Bessel functions ${\displaystyle \int J_1(z)dz=-J_0(z)+C,} {\displaystyle \int Y_1(z)dz=-Y_0(z)+C,}$ we obtain the expression of the normalized depth of groundwater as,

$$H(r,t)=\left(\frac{-Q}{2\pi T_f}\right)\ln (r)-{\displaystyle {\int }_0^{\infty }\frac{J_0(r\xi )Y_1(R\xi )-J_1(R\xi )Y_0(r\xi )}{J_1^2(R\xi )+Y_1^2(R\xi )}}·\widehat{U}(\xi ,t)d\xi$$

(43)

4. Conclusions and Discussions

The distributed-order fractional diffusion–wave equation for the radial flow to or from a well has been investigated by considering the time–fractional derivative with the Mittag–Leffler function of one parameter. The diffusion–wave equation is featured by two temporal fractional orders α and β to characterize different diffusion types.

It is known that Mittag–Leffler functions are solutions to some fractional differential equations; therefore, they are non-local functions. In consequence, the time–fractional Atangana–Baleanu derivative is a fractional derivative with non-singular and non-local kernel. It is expected that the non-locality of the Mittag–Leffler kernel allows better description of the memory within structure and media with different scales. Also, it is known that Mittag–Leffler functions are suitable for interesting generalizations. The fractional derivatives with Prabhakar kernel are suitable tools to describe complex physical processes, such as the anomalous dielectric relaxation or the renewal processes.

The analytical solutions of the initial-boundary value fractional problem have been determined using the Laplace and Dirichlet–Weber integral transforms. The solutions from the general case are suitable for generating the solutions corresponding to the ordinary case, namely for the case $\alpha =\beta =1$.

Recall that Su et al. [35] studied a similar problem by using Caputo time-fractional derivative. It is important to note that the Caputo fractional derivative cannot be obtained as a particular case of the Atangana–Baleanu derivative, so the results obtained in this paper cannot be identical to those in Su’s article. However, for the case α = β = 1, the results obtained by us are equivalent to those obtained by Su for the same values of the fractional parameters.

To highlight the influence of the fractional parameters α and β on the normalized depth of flow, numerical values of the function H(r,t) have been determined and plotted with the software Mathcad. In the performed analysis, we used the following values of the parameters: $D_f=0.5, Q=-0.5, T_f=20, R=0.5$.

In Figure 1 are the presented profiles of the normalized depth $H(r,t)$, versus r for a = b = 0.5 and for different values of the time t and the fractional parameters α and β. Note that the variation with radial distance of the normalized depth is not significant. The values of the fractional parameters α and β have a significant influence on the values of the function $H(r,t)$. For $\alpha >\beta$, the values of the normalized depth are higher than for $\alpha <\beta$. This is due to the different mode of the diffusion process in areas with different pores.

The influence of the fractional parameters α and β on the normalized depth $H(r,t)$ is analyzed in Figure 2 for a = 0.4 and b = 0.6 for different values of the time t. It is observed in Figure 2 that for a constant value of the fractional parameter α, the function $H(r,t)$ is increasing with the fractional parameter β. A similar behavior is observed for β constant and increasing values of the parameter α. The increasing of values of a fractional parameter means that the porosity of the medium is larger; therefore the diffusion process will be enhanced.

The influence of the parameter a on the normalized depth $H(r,t)$ is presented in Figure 3 for t = 6 and α = 0.885 and different values of the fractional parameter β. As expected, the normalized depth decreases with the parameter a. This result is due to a larger region of the medium with big pores where the flow is faster.