Fractional governing equations of transient groundwater flow in unconfined aquifers with multi-fractional dimensions in fractional time

In this study, a dimensionally consistent governing equation of transient unconfined groundwater flow in fractional time and multi-fractional space is developed. First, a fractional continuity equation for transient unconfined groundwater flow is developed in fractional time and space. For the equation of groundwater motion within a multi-fractional multidimensional unconfined aquifer, a previously developed dimensionally consistent equation for water flux in unsaturated/saturated porous media is combined with the Dupuit approximation to obtain an equation for groundwater motion in multi-fractional space in unconfined aquifers. Combining the fractional continuity and groundwater motion equations, the fractional governing equation of transient unconfined aquifer flow is then obtained. Finally, two numerical applications to unconfined aquifer groundwater flow are presented to show the skills of the proposed fractional governing equation. As shown in one of the numerical applications, the newly developed governing equation can produce heavy-tailed recession behavior in unconfined aquifer discharges.

1 Introduction

Nearly 70 years ago in his hydrologic studies of the Aswan High Dam, Hurst (1951) discovered that the flow time series of the Nile River demonstrated fluctuations whose rescaled range may not be proportional to the square root of the observation duration but may be proportional to the duration raised to a power H (the so-called Hurst coefficient) that is larger than 0.5 but less than 1. This finding, now called the “Hurst phenomenon” implies that in such river flows the integral scale (the integral of the flow autocorrelation function with respect to the time lag, over the range from zero to infinity) may not exist, putting the process outside the Brownian domain of finite-memory processes where the integral scale is finite. Since the Hurst phenomenon amounts to the clustering of wet years with wet years and dry years with dry years, the so-called “Joseph effect” from the Bible (Mandelbrot, 1977), it has important consequences on the planning and operation of water storage systems over long periods (Koutsoyiannis, 2005). The Hurst phenomenon in hydrologic flow processes was later demonstrated convincingly by various researchers, including Eltahir (1996), Radziejewski and Kundzewicz (1997), Montanari et al. (1997), and Vogel et. al. (1998) among others. In order to model the Hurst phenomenon in river flows, the fractional Gaussian noise (FGN), where the rescaled range for the time series of a flow process in a time interval [0,t] is proportional to t^H for $\mathrm{0.5}<H<\mathrm{1}$, was introduced by Mandelbrot and Wallis (1969). The FGN model was later extended by Koutsoyiannis (2002) in order to satisfactorily model a range of timescales, including the conventional Brownian finite-memory flow processes. Aside from the FGN models, physically based models of the Hurst phenomenon were also developed by various authors, including Klemes (1974), Beran (1994), and Koutsoyiannis (2003). However, a physically based model that explains the Hurst phenomenon explicitly in terms of the hydrologic process mechanisms is still missing. Yevjevich (1963, 1971) provided a plausible physical explanation for the Markovian structure of the annual river flows within a river basin by linking the annual evolution of the water storage in the basin to the exponential recession in baseflow of the basin runoff. Meanwhile, baseflow in the basin runoff is mainly due to unconfined aquifer flow to the neighboring stream network of the basin. As shall be shown in a numerical example later in this paper, the conventional unconfined groundwater flow equation with integer powers does result in the hydraulic head of and the discharge from the aquifer to decay exponentially, which would result in the Markovian finite-memory behavior of the river outflow from the basin. Such an exponentially decaying baseflow, while it can be explained by the mechanics of the conventional unconfined groundwater flow governing equation with integer powers, may not produce the heavy-tailed recession behavior necessary for the long-range dependence in river flows, the basic characteristic of the Hurst phenomenon, reported in annual river flow series in the abovementioned studies. The conventional integer-power governing equations of the unconfined groundwater flow, having finite memory, are fundamentally in the Brownian domain and may not model the heavy-tailed baseflow recession behavior that would be necessary to model the Hurst phenomenon in annual river flows. What is needed is a new structure for the governing equation of unconfined groundwater flow that can reproduce heavy-tailed behavior with time in the hydraulic head and aquifer discharge recession; this would then lead to heavy-tailed recession behavior in the baseflow of the river basin. Furthermore, various researchers also reported long-range dependence in groundwater level fluctuations (e.g., Li and Zhang, 2007; Yu et al., 2016; Tu et al., 2017; and the references therein). One possible way to reproduce heavy-tailed recession behavior in the hydraulic head and discharge of an unconfined aquifer is by means of a new governing equation of unconfined groundwater flow with fractional powers. Such behavior in an anisotropic confined groundwater aquifer with time and space fractional operators in its governing equation was recently demonstrated (Kavvas et al., 2017a; Tu et al., 2017). Accordingly, the study reported herein will follow a similar approach to develop a new governing equation for unconfined groundwater aquifers.

Reporting that conventional geometries cannot characterize groundwater flow in many fractured rock aquifers (Black et al., 1986) and the observed drawdown tends to be underestimated in early times and overestimated at later times by the conventional radial groundwater flow model (Van Tonder et al., 2001), Cloot and Botha (2006) developed a fractional governing equation for radial groundwater flow in integer time and fractional space in a uniform homogeneous aquifer. They used the Riemann–Liouville (RL) fractional derivative form (please see Podlubny, 1998, pp. 62–77, for a comprehensive explanation of the RL fractional derivative) in their model formulation. Atangana and Bildik (2013), Atangana (2014), and Atangana and Vermeulen (2014) then reformulated the fractional radial groundwater flow model of Cloot and Botha (2006) by using the Caputo differentiation framework (to be detailed in the next section) and reported better performance. Compared to the Riemann–Liouville derivative approach, the Caputo framework has a fundamental advantage of being able to accommodate physically interpretable real-life initial and boundary conditions (Podlubny, 1998). In simple terms, a differential equation that is based on the RL fractional derivative requires the limit values of the RL fractional derivative for its initial and boundary values, which have no known physical interpretation (Podlubny, 1998, p. 78). Meanwhile, “Caputo derivatives take on the same form as for integer-order differential equations, i.e., contain the limit values of integer-order derivatives” (Podlubny, 1998, p. 79), incorporating the real-world initial and boundary conditions into the solution to a fractional governing equation. Atangana and Baleanu (2014) presented a new radial groundwater flow model in fractional time based on a new fractional derivative definition, “conformable derivative” (Khalil et al., 2014). Most recently, Su (2017) proposed a time–space fractional Boussinesq equation, and he claimed this fractional equation is a general groundwater flow equation and can be applied to groundwater flow in both confined and unconfined aquifers. However, all of the aforementioned studies only presented the formulated fractional governing groundwater flow equations, and no detailed derivations of these governing equations from the fundamental conservation principles were provided.

Wheatcraft and Meerschaert (2008) derived the groundwater flow continuity equation in the fractional form by using the fractional Taylor series approximation. They further removed the linearity, or piecewise linearity, restriction for the flux and the infinitesimal control volume restriction. When developing the fractional continuity equation, the groundwater flow process was considered in fractional space but in integer time by Wheatcraft and Meerschaert (2008). They further assumed the same fractional power in every direction of the fractional porous media space. Furthermore, only the mass conservation was considered in their derivation, not the fractional water flux equation. Mehdinejadiani et al. (2013) expanded the approach of Wheatcraft and Meerschaert (2008) to the derivation of a governing equation of groundwater flow in an unconfined aquifer in fractional space but in integer time. In their derivation, they used the conventional Darcy formulation for the water flux with an integer spatial derivative, while utilizing fractional spatial derivatives in their continuity equation.

Olsen et al. (2016) pointed out that the derivations in Wheatcraft and Meerschaert (2008) and Mehdinejadiani et al. (2013) utilized the fractional Taylor series, as formulated by Odibat and Shawagfeh (2007), which utilized local Caputo derivatives. In order to expand the local Caputo derivatives in the abovementioned studies, Olsen et al. (2016) utilized the fractional mean value theorem from Diethelm (2012) to develop a continuity equation of groundwater flow with left and right fractional nonlocal Caputo derivatives in fractional space but in integer time. Olsen et al. (2016) did not address the water flux formulation in fractional space and, hence, did not develop a complete governing equation of groundwater flow. They also did not address the multi-fractional spatial derivatives in order to address anisotropy within an aquifer. Around that time, Kavvas et al. (2017a) utilized the mean value formulation from Usero (2008), Odibat and Shawagfeh (2007), and Li et al. (2009) to derive a complete governing equation of transient groundwater flow in an anisotropic confined aquifer with fractional time and multi-fractional space derivatives which addressed not only the continuity but also the water flux (motion) in fractional time–space and the effect of a sink/source term. By employing the abovementioned fractional mean value formulations, Kavvas et al. (2017a) developed the governing equation of confined groundwater flow in fractional time–space in nonlocal form.

As mentioned above, unconfined groundwater flow is the fundamental component of the watershed runoff baseflow since it is the fundamental contributor to the network streamflow within a watershed during dry periods. As such, the behavior of unconfined groundwater flow is key to the physically based understanding of the long memory in watershed runoff. Meanwhile, as will be seen in the following derivation of its governing equation, unconfined aquifer groundwater flow is uniquely different from the confined aquifer groundwater flow. The fundamental difference between the two aquifer flows is that while the flow in a confined aquifer is linear and compressible, the flow in an unconfined aquifer is nonlinear and incompressible due to the unconfined aquifer being phreatic, its top surface boundary being open to the atmosphere. Accordingly, hydrologists have developed unique governing equations of unconfined aquifer groundwater flow (Bear, 1979; Freeze and Cherry, 1979). Starting with the next section, first the continuity equation of transient unconfined groundwater flow within an anisotropic heterogeneous aquifer under a time–space varying sink/source will be developed in fractional time and fractional space. Then, this fractional continuity equation will be combined with a fractional groundwater motion equation to obtain a transient groundwater flow equation in fractional time and multi-fractional space within an anisotropic heterogeneous unconfined aquifer.

Analogous to the traditional governing groundwater flow equations, as outlined by Freeze and Cherry (1979) and Bear (1979), the fractional unconfined groundwater flow equations must have specific features (Kavvas et al., 2017a):

• i.

  In order for the governing equation to be prognostic, the form of the equation must be known completely from the outset.

• ii.

  The fractional governing equations must be dimensionally consistent and be purely differential equations, containing only differential operators without difference operators.

• iii.

  As the fractional derivative powers go to integer values, the fractional unconfined groundwater flow equations must converge to the corresponding conventional integer-order governing equations.

Within this framework, the governing equations of unconfined groundwater flow in fractional time and fractional space will be developed in the following.

2 Derivation of the continuity equation for transient unconfined groundwater flow in a heterogeneous anisotropic multi-fractional medium in fractional time

The β-order Caputo fractional derivative $D_a^{k\mathit{\beta }}f\left(x\right)$ of a function f(x) may be defined as (Odibat and Shawagfeh, 2007; Podlubny, 1998; Usero, 2008, and Li et al., 2009)

$$\begin{array}{}\text{(1)}& \begin{array}{rl}& D_a^{\mathit{\beta }}f\left(x\right)={\displaystyle \frac{\mathrm{1}}{\mathrm{\Gamma }(\mathrm{1}-\mathit{\beta })}}\underset{a}{\overset{x}{\int }}{\displaystyle \frac{f^{\prime }\left(\mathit{\xi }\right)}{(x-\mathit{\xi }{)}^{\mathit{\beta }}}}\mathrm{d}\mathit{\xi },\\ & \mathrm{0}<\mathit{\beta }<\mathrm{1},x\ge a,\end{array}\end{array}$$

where ξ represents a dummy variable in the equation.

It was shown in Kavvas et al. (2017b) that one can obtain a ${\mathit{\beta }}_{x_i}$-order approximation (i=1, 2) to a function f(x_i) around x_i−Δx_i by the following:

$$\begin{array}{}\text{(2)}& \begin{array}{rl}& f\left(x_i\right)=f\left(x_i-\mathrm{\Delta }{x}_i\right)+{\displaystyle \frac{{\left(\mathrm{\Delta }{x}_i\right)}^{{\mathit{\beta }}_{x_i}}}{\mathrm{\Gamma }\left({\mathit{\beta }}_{x_i}+\mathrm{1}\right)}}{D}_{x_i-\mathrm{\Delta }{x}_i}^{{\mathit{\beta }}_{x_i}}f\left(x_i\right);\\ & i=\mathrm{1},\mathrm{2}.\end{array}\end{array}$$

In Eq. (2), an analytical relationship between Δx_i and ${\left(\mathrm{\Delta }{x}_i\right)}^{{\mathit{\beta }}_{x_i}}$ (i=1, 2) that will be universally applicable throughout the modeling domain is possible when the lower limit of the above Caputo derivative in Eq. (2) is taken as zero (that is, Δx_i=x_i) for f(x_i)=x_i (Kavvas et al., 2017b).

Figure 1 The mass flux through the control volume of an unconfined aquifer.

Under the Dupuit approximation of horizontal flow streamlines (for a very small water table gradient) (Bear, 1979), the net mass flux through the control volume of an unconfined aquifer with a flat bottom confining layer, as depicted in Fig. 1, which also has a sink/source mass flux, ρq_vΔx₁Δx₂, can be formulated as

$$\begin{array}{}\text{(3)}& \begin{array}{rl}& \left[\mathit{\rho }{Q}_{x_{\mathrm{1}}}\left(x_{\mathrm{1}},x_{\mathrm{2}};t\right)-\mathit{\rho }{Q}_{x_{\mathrm{1}}}\left(x_{\mathrm{1}}-\mathrm{\Delta }{x}_{\mathrm{1}},x_{\mathrm{2}};t\right)\right]\mathrm{\Delta }{x}_{\mathrm{2}}\\ & +\left[\mathit{\rho }{Q}_{x_{\mathrm{2}}}\left(x_{\mathrm{1}},x_{\mathrm{2}};t\right)-\mathit{\rho }{Q}_{x_{\mathrm{2}}}\left(x_{\mathrm{1}},x_{\mathrm{2}}-\mathrm{\Delta }{x}_{\mathrm{2}};t\right)\right]\mathrm{\Delta }{x}_{\mathrm{1}}\\ & -\mathit{\rho }{q}_v\mathrm{\Delta }{x}_{\mathrm{1}}\mathrm{\Delta }{x}_{\mathrm{2}},\end{array}\end{array}$$

where $Q_{x_i}$ is the discharge across a vertical plane of unit width in the ith direction (i=1, 2), ρ is the fluid density, and q_v is the source/sink (recharge/leakage) per unit horizontal area. Then, combining Eq. (2) with Eq. (3), with Δx_i=x_i (i=1, 2), and expressing the resulting Caputo derivative, $D_{\mathrm{0}}^{{\mathit{\beta }}_{x_i}}f\left(x_i\right)$, as $\frac{{\partial }^{{\mathit{\beta }}_{x_i}}f\left(x_i\right)}{{\left(\partial x_i\right)}^{{\mathit{\beta }}_{x_i}}}$ (i=1, 2) for convenience, yields the net mass flux through the control volume in Fig. 1, to the orders of ${\left(\mathrm{\Delta }{x}_{\mathrm{1}}\right)}^{{\mathit{\beta }}_{x_{\mathrm{1}}}}$ and ${\left(\mathrm{\Delta }{x}_{\mathrm{2}}\right)}^{{\mathit{\beta }}_{x_{\mathrm{2}}}}$, expressed as

$$\begin{array}{}\text{(4)}& \begin{array}{rl}& {\displaystyle \frac{\mathrm{1}}{\mathrm{\Gamma }\left({\mathit{\beta }}_{x_{\mathrm{1}}}+\mathrm{1}\right)}}{\left({\displaystyle \frac{\partial }{\partial x_{\mathrm{1}}}}\right)}^{{\mathit{\beta }}_{x_{\mathrm{1}}}}\left(\mathit{\rho }{Q}_{x_{\mathrm{1}}}\left(x_{\mathrm{1}},x_{\mathrm{2}};t\right)\right){\left(\mathrm{\Delta }{x}_{\mathrm{1}}\right)}^{{\mathit{\beta }}_{x_{\mathrm{1}}}}\mathrm{\Delta }{x}_{\mathrm{2}}\\ & +{\displaystyle \frac{\mathrm{1}}{\mathrm{\Gamma }\left({\mathit{\beta }}_{x_{\mathrm{2}}}+\mathrm{1}\right)}}{\left({\displaystyle \frac{\partial }{\partial x_{\mathrm{2}}}}\right)}^{{\mathit{\beta }}_{x_{\mathrm{2}}}}\left(\mathit{\rho }{Q}_{x_{\mathrm{2}}}\left(x_{\mathrm{1}},x_{\mathrm{2}};t\right)\right)\mathrm{\Delta }{x}_{\mathrm{1}}{\left(\mathrm{\Delta }{x}_{\mathrm{2}}\right)}^{{\mathit{\beta }}_{x_{\mathrm{2}}}}\\ & -\mathit{\rho }{q}_v\mathrm{\Delta }{x}_{\mathrm{1}}\mathrm{\Delta }{x}_{\mathrm{2}},\end{array}\end{array}$$

where different powers for fractional space derivatives are utilized in different directions due to the anisotropy in the flow medium.

Kavvas et al. (2017b) have shown that, to ${\mathit{\beta }}_{x_i}$-order fractional increments in space in the ith direction (i=1, 2),

$$\begin{array}{}\text{(5)}& {\left(\mathrm{\Delta }{x}_i\right)}^{{\mathit{\beta }}_{x_i}}={\displaystyle \frac{\mathrm{\Gamma }\left({\mathit{\beta }}_{x_i}+\mathrm{1}\right)\mathrm{\Gamma }(\mathrm{2}-{\mathit{\beta }}_{x_i})}{x_i^{\mathrm{1}-{\mathit{\beta }}_{x_i}}}}\mathrm{\Delta }{x}_i,i=\mathrm{1},\mathrm{2}.\end{array}$$

Combining Eqs. (5) and (4) yields, for the net mass outflow through the control volume in Fig. 1 (to the orders of ${\left(\mathrm{\Delta }{x}_i\right)}^{{\mathit{\beta }}_{x_i}}$, i=1, 2),

$$\begin{array}{}\text{(6)}& \begin{array}{rl}& {\displaystyle \frac{\mathrm{\Gamma }\left(\mathrm{2}-{\mathit{\beta }}_{x_{\mathrm{1}}}\right)}{x_{\mathrm{1}}^{\mathrm{1}-{\mathit{\beta }}_{x_{\mathrm{1}}}}}}{\left({\displaystyle \frac{\partial }{\partial x_{\mathrm{1}}}}\right)}^{{\mathit{\beta }}_{x_{\mathrm{1}}}}\left(\mathit{\rho }{Q}_{x_{\mathrm{1}}}\left(\stackrel{\mathrm{‾}}{x};t\right)\right)\mathrm{\Delta }{x}_{\mathrm{1}}\mathrm{\Delta }{x}_{\mathrm{2}}\\ & +{\displaystyle \frac{\mathrm{\Gamma }(\mathrm{2}-{\mathit{\beta }}_{x_{\mathrm{2}}})}{x_{\mathrm{2}}^{\mathrm{1}-{\mathit{\beta }}_{x_{\mathrm{2}}}}}}{\left({\displaystyle \frac{\partial }{\partial x_{\mathrm{2}}}}\right)}^{{\mathit{\beta }}_{x_{\mathrm{2}}}}\left(\mathit{\rho }{Q}_{x_{\mathrm{2}}}\left(\stackrel{\mathrm{‾}}{x};t\right)\right)\mathrm{\Delta }{x}_{\mathrm{1}}\mathrm{\Delta }{x}_{\mathrm{2}}\\ & -\mathit{\rho }{q}_v\mathrm{\Delta }{x}_{\mathrm{1}}\mathrm{\Delta }{x}_{\mathrm{2}},\stackrel{\mathrm{‾}}{x}=\left(x_{\mathrm{1}},x_{\mathrm{2}}\right).\end{array}\end{array}$$

Denoting the water volume within the control volume in Fig. 1 by V_w, the specific yield (effective porosity) S_y of a phreatic aquifer (Bear and Verruijt, 1987) is expressed as

$$\begin{array}{}\text{(7)}& S_y={\displaystyle \frac{\mathrm{\Delta }{V}_{\mathrm{w}}}{\mathrm{\Delta }h}}{\displaystyle \frac{\mathrm{1}}{\mathrm{\Delta }{x}_{\mathrm{1}}\mathrm{\Delta }{x}_{\mathrm{2}}}},\end{array}$$

where ΔV_w is the change in water volume in the control volume per change Δh in the hydraulic head (the elevation of the phreatic surface (water table) above the flat bottom of the aquifer). The time rate of change in mass within the control volume in Fig. 1 may be written as (Bear and Verruijt, 1987)

$$\begin{array}{}\text{(8)}& {\displaystyle \frac{S_y\left(\mathit{\rho }h\left(\stackrel{\mathrm{‾}}{x};t\right)-\mathit{\rho }h\left(\stackrel{\mathrm{‾}}{x};t-\mathrm{\Delta }t\right)\right)}{\mathrm{\Delta }t}}\mathrm{\Delta }{x}_{\mathrm{1}}\mathrm{\Delta }{x}_{\mathrm{2}},\end{array}$$

which can then be expressed in terms of the approximation (Eq. 2) with respect to the time dimension as

$$\begin{array}{}\text{(9)}& {\displaystyle \frac{S_y}{\mathrm{\Delta }t}}\left[{\displaystyle \frac{\mathrm{\Delta }{t}^{\mathit{\alpha }}}{\mathrm{\Gamma }(\mathit{\alpha }+\mathrm{1})}}{\left({\displaystyle \frac{\partial }{\partial t}}\right)}^{\mathit{\alpha }}\left(\mathit{\rho }h\right)\right]\mathrm{\Delta }{x}_{\mathrm{1}}\mathrm{\Delta }{x}_{\mathrm{2}}.\end{array}$$

To α-order fractional increments in time (Kavvas et al., 2017b)

$$\begin{array}{}\text{(10)}& {\left(\mathrm{\Delta }t\right)}^{\mathit{\alpha }}={\displaystyle \frac{\mathrm{\Gamma }\left(\mathit{\alpha }+\mathrm{1}\right)\mathrm{\Gamma }(\mathrm{2}-\mathit{\alpha })}{t^{\mathrm{1}-\mathit{\alpha }}}}\mathrm{\Delta }t.\end{array}$$

Substituting Eq. (10) into Eq. (9), one can obtain the time rate of change of mass in the control volume, as shown in Fig. 1,

$$\begin{array}{}\text{(11)}& S_y{\displaystyle \frac{\mathrm{\Gamma }(\mathrm{2}-\mathit{\alpha })}{t^{\mathrm{1}-\mathit{\alpha }}}}{\left({\displaystyle \frac{\partial }{\partial t}}\right)}^{\mathit{\alpha }}\left(\mathit{\rho }h\right)\mathrm{\Delta }{x}_{\mathrm{1}}\mathrm{\Delta }{x}_{\mathrm{2}}.\end{array}$$

As the time rate of change of mass within the control volume, as shown in Fig. 1, must be inversely proportional to the net mass flux passing through the control volume, one may combine Eqs. (6) and (11) to obtain

$$\begin{array}{}\text{(12)}& {\displaystyle }\begin{array}{rl}& [{\displaystyle \frac{\mathrm{\Gamma }\left(\mathrm{2}-{\mathit{\beta }}_{x_{\mathrm{1}}}\right)}{x_{\mathrm{1}}^{\mathrm{1}-{\mathit{\beta }}_{x_{\mathrm{1}}}}}}{\left({\displaystyle \frac{\partial }{\partial x_{\mathrm{1}}}}\right)}^{{\mathit{\beta }}_{x_{\mathrm{1}}}}\left(\mathit{\rho }{Q}_{x_{\mathrm{1}}}\left(\stackrel{\mathrm{‾}}{x};t\right)\right)\\ & +{\displaystyle \frac{\mathrm{\Gamma }\left(\mathrm{2}-{\mathit{\beta }}_{x_{\mathrm{2}}}\right)}{x_{\mathrm{2}}^{\mathrm{1}-{\mathit{\beta }}_{x_{\mathrm{2}}}}}}{\left({\displaystyle \frac{\partial }{\partial x_{\mathrm{2}}}}\right)}^{{\mathit{\beta }}_{x_{\mathrm{2}}}}\left(\mathit{\rho }{Q}_{x_{\mathrm{2}}}\left(\stackrel{\mathrm{‾}}{x};t\right)\right)\\ & -\mathit{\rho }{q}_v]\mathrm{\Delta }{x}_{\mathrm{1}}\mathrm{\Delta }{x}_{\mathrm{2}}=\\ & -S_y{\displaystyle \frac{\mathrm{\Gamma }(\mathrm{2}-\mathit{\alpha })}{t^{\mathrm{1}-\mathit{\alpha }}}}{\left({\displaystyle \frac{\partial }{\partial t}}\right)}^{\mathit{\alpha }}\left(\mathit{\rho }h\right)\mathrm{\Delta }{x}_{\mathrm{1}}\mathrm{\Delta }{x}_{\mathrm{2}};\end{array}\text{(13)}& {\displaystyle }\begin{array}{rl}& {\displaystyle \frac{\mathrm{\Gamma }\left(\mathrm{2}-{\mathit{\beta }}_{x_{\mathrm{1}}}\right)}{x_{\mathrm{1}}^{\mathrm{1}-{\mathit{\beta }}_{x_{\mathrm{1}}}}}}{\left({\displaystyle \frac{\partial }{\partial x_{\mathrm{1}}}}\right)}^{{\mathit{\beta }}_{x_{\mathrm{1}}}}\left(\mathit{\rho }{Q}_{x_{\mathrm{1}}}\left(\stackrel{\mathrm{‾}}{x};t\right)\right)\\ & +{\displaystyle \frac{\mathrm{\Gamma }\left(\mathrm{2}-{\mathit{\beta }}_{x_{\mathrm{2}}}\right)}{x_{\mathrm{2}}^{\mathrm{1}-{\mathit{\beta }}_{x_{\mathrm{2}}}}}}{\left({\displaystyle \frac{\partial }{\partial x_{\mathrm{2}}}}\right)}^{{\mathit{\beta }}_{x_{\mathrm{2}}}}\left(\mathit{\rho }{Q}_{x_{\mathrm{2}}}\left(\stackrel{\mathrm{‾}}{x};t\right)\right)-\mathit{\rho }{q}_v=\\ & -S_y{\displaystyle \frac{\mathrm{\Gamma }(\mathrm{2}-\mathit{\alpha })}{t^{\mathrm{1}-\mathit{\alpha }}}}{\left({\displaystyle \frac{\partial }{\partial t}}\right)}^{\mathit{\alpha }}\left(\mathit{\rho }h\right)\end{array}\end{array}$$

for $\mathrm{0}<\mathit{\alpha },{\mathit{\beta }}_{x_{\mathrm{1}}},{\mathit{\beta }}_{x_{\mathrm{2}}}<\mathrm{1}$, $\stackrel{\mathrm{‾}}{x}=\left(x_{\mathrm{1}},x_{\mathrm{2}},\right)$.

Within the framework of fluid incompressibility in the unconfined aquifer, Eq. (13) reduces further to

$$\begin{array}{}\text{(14)}& \begin{array}{rl}& {\displaystyle \frac{\mathrm{\Gamma }\left(\mathrm{2}-{\mathit{\beta }}_{x_{\mathrm{1}}}\right)}{\mathrm{\Gamma }(\mathrm{2}-\mathit{\alpha })}}{\displaystyle \frac{t^{\mathrm{1}-\mathit{\alpha }}}{x_{\mathrm{1}}^{\mathrm{1}-{\mathit{\beta }}_{x_{\mathrm{1}}}}}}{\left({\displaystyle \frac{\partial }{\partial x_{\mathrm{1}}}}\right)}^{{\mathit{\beta }}_{x_{\mathrm{1}}}}\left(Q_{x_{\mathrm{1}}}\left(\stackrel{\mathrm{‾}}{x};t\right)\right)\\ & +{\displaystyle \frac{\mathrm{\Gamma }\left(\mathrm{2}-{\mathit{\beta }}_{x_{\mathrm{2}}}\right)}{\mathrm{\Gamma }(\mathrm{2}-\mathit{\alpha })}}{\displaystyle \frac{t^{\mathrm{1}-\mathit{\alpha }}}{x_{\mathrm{2}}^{\mathrm{1}-{\mathit{\beta }}_{x_{\mathrm{2}}}}}}{\left({\displaystyle \frac{\partial }{\partial x_{\mathrm{2}}}}\right)}^{{\mathit{\beta }}_{x_{\mathrm{2}}}}\left(Q_{x_{\mathrm{2}}}\left(\stackrel{\mathrm{‾}}{x};t\right)\right)\\ & -{\displaystyle \frac{t^{\mathrm{1}-\mathit{\alpha }}}{\mathrm{\Gamma }(\mathrm{2}-\mathit{\alpha })}}{q}_v=-S_y{\displaystyle \frac{{\partial }^{\mathit{\alpha }}h}{{\left(\partial t\right)}^{\mathit{\alpha }}}},\end{array}\end{array}$$

for $\mathrm{0}<\mathit{\alpha },{\mathit{\beta }}_{x_{\mathrm{1}}},{\mathit{\beta }}_{x_{\mathrm{2}}}<\mathrm{1}$, $\stackrel{\mathrm{‾}}{x}=\left(x_{\mathrm{1}},x_{\mathrm{2}},\right)$, as the time–space fractional continuity equation of transient groundwater flow in an anisotropic unconfined aquifer with multi-fractional dimensions and in fractional time.

Performing a dimensional analysis of Eq. (14) yields

$$\begin{array}{}\text{(15)}& \begin{array}{rl}{\displaystyle \frac{L}{T^{\mathit{\alpha }}}}& ={\displaystyle \frac{T^{\mathrm{1}-\mathit{\alpha }}}{L^{\mathrm{1}-{\mathit{\beta }}_{x_{\mathrm{1}}}}}}{\displaystyle \frac{\mathrm{1}}{L^{{\mathit{\beta }}_{x_{\mathrm{1}}}}}}{\displaystyle \frac{L^{\mathrm{2}}}{T}}={\displaystyle \frac{T^{\mathrm{1}-\mathit{\alpha }}}{L^{\mathrm{1}-{\mathit{\beta }}_{x_{\mathrm{2}}}}}}{\displaystyle \frac{\mathrm{1}}{L^{{\mathit{\beta }}_{x_{\mathrm{2}}}}}}{\displaystyle \frac{L^{\mathrm{2}}}{T}}\\ & ={\displaystyle \frac{T^{\mathrm{1}-\mathit{\alpha }}}{\mathrm{1}}}{\displaystyle \frac{L}{T}}={\displaystyle \frac{L}{T^{\mathit{\alpha }}}},\end{array}\end{array}$$

where L denotes length and T denotes time. Also, α, ${\mathit{\beta }}_{x_{\mathrm{1}}}$, and ${\mathit{\beta }}_{x_{\mathrm{2}}}$ are respectively the fractional powers in time and in x₁ and x₂ spatial dimensions. In Eq. (15), starting from the left-hand side (LHS), the first term shows the final dimensions of Eq. (14), the second term shows in detail the dimensions of the individual components of the first term on the LHS of Eq. (14), the third term shows in detail the dimensions of the individual components of the second term on the LHS of Eq. (14), the fourth term shows in detail the dimensions of the individual components of the third term on the LHS of Eq. (14), and the fifth term shows in detail the dimensions of the individual components on the right-hand side (RHS) of Eq. (14). Hence, the left-hand and right-hand sides of the continuity Eq. (14) for transient groundwater flow in an unconfined aquifer in multi-fractional space and fractional time are consistent, as shown in Eq. (15).

For $n-\mathrm{1}<\mathit{\alpha },{\mathit{\beta }}_{x_i}<n$, where n is any positive integer, as α and ${\mathit{\beta }}_{x_i}\to n$, the Caputo fractional derivative of a function f(y) to order α or ${\mathit{\beta }}_{x_i}$ (i=1, 2) yields the standard nth derivative of the function f(y) (Podlubny, 1998). Then, when α and ${\mathit{\beta }}_{x_i}\to \mathrm{1}$ (i=1, 2), the continuity Eq. (14) becomes the conventional continuity equation for transient groundwater flow in an unconfined aquifer:

$$\begin{array}{}\text{(16)}& -S_y{\displaystyle \frac{\partial h}{\partial t}}={\displaystyle \frac{\partial }{\partial x_{\mathrm{1}}}}\left(Q_{x_{\mathrm{1}}}\left(\stackrel{\mathrm{‾}}{x};t\right)\right)+{\displaystyle \frac{\partial }{\partial x_{\mathrm{2}}}}\left(Q_{x_{\mathrm{2}}}\left(\stackrel{\mathrm{‾}}{x};t\right)\right)-q_v.\end{array}$$

3 Motion equation (specific discharge equation) in fractional multidimensional unconfined aquifers

Recently, Kavvas et al. (2017a, b) derived a governing equation for water flux (specific discharge), $q_{x_i}$ (i=1, 2, 3), in a saturated or unsaturated porous medium with fractional dimensions in the form

$$\begin{array}{}\text{(17)}& \begin{array}{rl}& q_i\left(\stackrel{\mathrm{‾}}{x},t\right)=-K_{\mathrm{s},x_i}\left(\stackrel{\mathrm{‾}}{x}\right){\displaystyle \frac{\mathrm{\Gamma }\left(\mathrm{2}-{\mathit{\beta }}_{x_i}\right)}{x_i^{\mathrm{1}-{\mathit{\beta }}_{x_i}}}}{\displaystyle \frac{{\partial }^{{\mathit{\beta }}_{x_i}}h}{{\left(\partial x_i\right)}^{{\mathit{\beta }}_{x_i}}}},\\ & i=\mathrm{1},\mathrm{2},\mathrm{3};\stackrel{\mathrm{‾}}{x}=\left(x_{\mathrm{1}},x_{\mathrm{2}},x_{\mathrm{3}}\right),\end{array}\end{array}$$

where $K_{\mathrm{s},x_i}\left(\stackrel{\mathrm{‾}}{x}\right)$ is the saturated hydraulic conductivity in the ith spatial direction (i=1, 2, 3). Meanwhile, under the Dupuit approximation of essentially horizontal unconfined aquifer flow (for a very small water table slope) (Bear, 1979), referring to Fig. 1, the discharge per unit width in the ith direction (i=1, 2) can be expressed as

$$\begin{array}{}\text{(18)}& Q_{x_i}\left(\stackrel{\mathrm{‾}}{x},t\right)=h{q}_i\left(\stackrel{\mathrm{‾}}{x},t\right),i=\mathrm{1},\mathrm{2};\stackrel{\mathrm{‾}}{x}=\left(x_{\mathrm{1}},x_{\mathrm{2}},\right).\end{array}$$

Then combining Eqs. (18) and (17) results in

$$\begin{array}{}\text{(19)}& \begin{array}{rl}& Q_{x_i}\left(\stackrel{\mathrm{‾}}{x},t\right)=-K_{\mathrm{s},x_i}\left(\stackrel{\mathrm{‾}}{x}\right){\displaystyle \frac{\mathrm{\Gamma }\left(\mathrm{2}-{\mathit{\beta }}_{x_i}\right)}{x_i^{\mathrm{1}-{\mathit{\beta }}_{x_i}}}}h{\displaystyle \frac{{\partial }^{{\mathit{\beta }}_{x_i}}h}{{\left(\partial x_i\right)}^{{\mathit{\beta }}_{x_i}}}},\\ & i=\mathrm{1},\mathrm{2};\stackrel{\mathrm{‾}}{x}=\left(x_{\mathrm{1}},x_{\mathrm{2}}\right),\end{array}\end{array}$$

as the governing equation of groundwater motion within an unconfined aquifer with a flat bottom confining layer. In Eq. (19) h is the unconfined aquifer thickness or the phreatic surface elevation above the bottom confining layer.

A dimensional analysis of Eq. (19) yields L²∕T for the units of both the LHS and the RHS of the equation, establishing its dimensional consistency.

Applying the abovementioned result of Podlubny (1998) to the convergence of a fractional derivative to a corresponding integer derivative for ${\mathit{\beta }}_{x_i}\to \mathrm{1}$ (i=1, 2) reduces the fractional motion Eq. (19) for unconfined groundwater flow to the conventional equation (Bear, 1979):

$$\begin{array}{}\text{(20)}& Q_{x_i}\left(\stackrel{\mathrm{‾}}{x},t\right)=-K_{\mathrm{s},x_i}\left(\stackrel{\mathrm{‾}}{x}\right)h{\displaystyle \frac{\partial h\left(\stackrel{\mathrm{‾}}{x},t\right)}{\partial x_i}},i=\mathrm{1},\mathrm{2},\end{array}$$

for the case of integer spatial dimensions. As such, the fractional motion Eq. (19) for unconfined groundwater flow in fractional spatial dimensions is consistent with the conventional motion equation for the integer spatial dimensions.

4 The complete equation for transient unconfined groundwater flow in multi-fractional space and fractional time

Combining the fractional motion Eq. (19) of groundwater flow in an unconfined aquifer with the fractional continuity Eq. (14) of unconfined groundwater flow results in the equation

$$\begin{array}{}\text{(21)}& \begin{array}{rl}& S_y{\displaystyle \frac{{\partial }^{\mathit{\alpha }}h}{{\left(\partial t\right)}^{\mathit{\alpha }}}}={\displaystyle \frac{\mathrm{\Gamma }\left(\mathrm{2}-{\mathit{\beta }}_{x_{\mathrm{1}}}\right)}{x_{\mathrm{1}}^{\mathrm{1}-{\mathit{\beta }}_{x_{\mathrm{1}}}}}}{\left({\displaystyle \frac{\partial }{\partial x_{\mathrm{1}}}}\right)}^{{\mathit{\beta }}_{x_{\mathrm{1}}}}\\ & \left(K_{\mathrm{s},x_{\mathrm{1}}}\left(\stackrel{\mathrm{‾}}{x}\right){\displaystyle \frac{t^{\mathrm{1}-\mathit{\alpha }}}{x_{\mathrm{1}}^{\mathrm{1}-{\mathit{\beta }}_{x_{\mathrm{1}}}}}}{\displaystyle \frac{\mathrm{\Gamma }\left(\mathrm{2}-{\mathit{\beta }}_{x_{\mathrm{1}}}\right)}{\mathrm{\Gamma }(\mathrm{2}-\mathit{\alpha })}}h{\displaystyle \frac{{\partial }^{{\mathit{\beta }}_{x_{\mathrm{1}}}}h}{{\left(\partial x_{\mathrm{1}}\right)}^{{\mathit{\beta }}_{x_{\mathrm{1}}}}}}\right)\\ & +{\displaystyle \frac{\mathrm{\Gamma }\left(\mathrm{2}-{\mathit{\beta }}_{x_{\mathrm{2}}}\right)}{x_{\mathrm{2}}^{\mathrm{1}-{\mathit{\beta }}_{x_{\mathrm{2}}}}}}{\left({\displaystyle \frac{\partial }{\partial x_{\mathrm{2}}}}\right)}^{{\mathit{\beta }}_{x_{\mathrm{2}}}}\\ & \left(K_{\mathrm{s},x_{\mathrm{2}}}\left(\stackrel{\mathrm{‾}}{x}\right){\displaystyle \frac{t^{\mathrm{1}-\mathit{\alpha }}}{x_{\mathrm{2}}^{\mathrm{1}-{\mathit{\beta }}_{x_{\mathrm{2}}}}}}{\displaystyle \frac{\mathrm{\Gamma }\left(\mathrm{2}-{\mathit{\beta }}_{x_{\mathrm{2}}}\right)}{\mathrm{\Gamma }(\mathrm{2}-\mathit{\alpha })}}h{\displaystyle \frac{{\partial }^{{\mathit{\beta }}_{x_{\mathrm{2}}}}h}{{\left(\partial x_{\mathrm{2}}\right)}^{{\mathit{\beta }}_{x_{\mathrm{2}}}}}}\right)\\ & +{\displaystyle \frac{t^{\mathrm{1}-\mathit{\alpha }}}{\mathrm{\Gamma }(\mathrm{2}-\mathit{\alpha })}}{q}_v,\end{array}\end{array}$$

for $\mathrm{0}<\mathit{\alpha },{\mathit{\beta }}_{x_{\mathrm{1}}},{\mathit{\beta }}_{x_{\mathrm{2}}}<\mathrm{1}$, $\stackrel{\mathrm{‾}}{x}=\left(x_{\mathrm{1}},x_{\mathrm{2}},\right)$, as the time–space fractional governing equation of transient unconfined groundwater flow in an anisotropic medium.

Performing a dimensional analysis of Eq. (21) yields

$$\begin{array}{}\text{(22)}& \begin{array}{rl}{\displaystyle \frac{L}{T^{\mathit{\alpha }}}}& ={\displaystyle \frac{\mathrm{1}}{L^{\mathrm{1}-{\mathit{\beta }}_{x_{\mathrm{1}}}}}}{\displaystyle \frac{\mathrm{1}}{L^{{\mathit{\beta }}_{x_{\mathrm{1}}}}}}{\displaystyle \frac{L}{T}}{\displaystyle \frac{T^{\mathrm{1}-\mathit{\alpha }}}{L^{\mathrm{1}-{\mathit{\beta }}_{x_{\mathrm{1}}}}}}L{\displaystyle \frac{L}{L^{{\mathit{\beta }}_{x_{\mathrm{1}}}}}}\\ & ={\displaystyle \frac{\mathrm{1}}{L^{\mathrm{1}-{\mathit{\beta }}_{x_{\mathrm{2}}}}}}{\displaystyle \frac{\mathrm{1}}{L^{{\mathit{\beta }}_{x_{\mathrm{2}}}}}}{\displaystyle \frac{L}{T}}{\displaystyle \frac{T^{\mathrm{1}-\mathit{\alpha }}}{L^{\mathrm{1}-{\mathit{\beta }}_{x_{\mathrm{2}}}}}}{\displaystyle \frac{L^{\mathrm{2}}}{L^{{\mathit{\beta }}_{x_{\mathrm{2}}}}}}={\displaystyle \frac{T^{\mathrm{1}-\mathit{\alpha }}}{\mathrm{1}}}{\displaystyle \frac{L}{T}}={\displaystyle \frac{L}{T^{\mathit{\alpha }}}},\end{array}\end{array}$$

where L denotes length and T denotes time. Hence, the left-hand and right-hand sides of the governing Eq. (21) for transient groundwater flow in an unconfined aquifer in multi-fractional space and fractional time are consistent.

Specializing the above-discussed result of Podlubny (1998) to n=1, for α and ${\mathit{\beta }}_{x_i}\to \mathrm{1}$ (i=1, 2), reduces the governing fractional Eq. (21) to the conventional governing equation for transient groundwater flow in an unconfined aquifer (Bear, 1979):

$$\begin{array}{}\text{(23)}& \begin{array}{rl}& S_y{\displaystyle \frac{\partial h}{\partial t}}={\displaystyle \frac{\partial }{\partial x_{\mathrm{1}}}}\left(K_{\mathrm{s},x_{\mathrm{1}}}\left(\stackrel{\mathrm{‾}}{x}\right)h{\displaystyle \frac{\partial h(\stackrel{\mathrm{‾}}{x},t)}{\partial x_{\mathrm{1}}}}\right)\\ & +{\displaystyle \frac{\partial }{\partial x_{\mathrm{2}}}}\left(K_{\mathrm{s},x_{\mathrm{2}}}\left(\stackrel{\mathrm{‾}}{x}\right)h{\displaystyle \frac{\partial h(\stackrel{\mathrm{‾}}{x},t)}{\partial x_{\mathrm{2}}}}\right)+q_v.\end{array}\end{array}$$

Figure 2 The sketch of numerical application 1; water seepage through the body of a dam as an unconfined groundwater flow.

Figure 3 Results for numerical application 1: (a) the normalized groundwater head at $x=L/\mathrm{2}$ through time under different fractional derivative powers; (b) the normalized groundwater discharge per unit width at $x=L/\mathrm{2}$ through time under different fractional derivative powers; t is time and the simulation time T is 120 000 min; (c) the normalized groundwater head at t=40 000 min through length of the aquifer (through the body of the dam) and the analytical solution to the conventional governing equation of unconfined groundwater flow when ${\mathit{\beta }}_x=\mathit{\alpha }=\mathrm{1}$ at the steady state.

5 Numerical application

To demonstrate the skills of the proposed fractional governing equation of unconfined aquifer groundwater flow, two numerical applications are performed using the proposed fractional governing equation. The first application (numerical application 1) follows the physical setting of an example from Wang and Anderson (1995), as depicted in Fig. 2. The numerical problem of seepage through a dam under a sudden change in the water surface elevation at the downstream section of the dam is modified based on seepage through a dam, p. 53 and Problem 4.4(a), p. 89 in Wang and Anderson (1995), as shown in Fig. 2. The water seepage through the body of the dam may be interpreted as one-dimensional groundwater flow through an unconfined aquifer. The unconfined flow system locates the top boundary of the saturated zone in an earthen dam and the bottom of the dam rests on impermeable rock. For this example, the unconfined aquifer length L is 100 m. The initial water level in the upstream and downstream sections of the dam and through the body of the dam is 16 m. Then immediately after the initial time, the water level at the downstream section of the dam is suddenly dropped to 11 m and remains at 11 m afterwards. The unconfined aquifer parameters are S=0.2 for the specific yield and K=0.002 m min⁻¹ for the hydraulic conductivity. The analytical solution to this problem at the steady state is

$$\begin{array}{}\text{(24)}& h=\sqrt{{\displaystyle \frac{h_{\mathrm{2}}^{\mathrm{2}}-h_{\mathrm{1}}^{\mathrm{2}}}{L}}x+h_{\mathrm{1}}^{\mathrm{2}}},\end{array}$$

where h is the depth of the unconfined groundwater surface from the bottom layer, L is the aquifer length, x is the distance from the upstream location of the dam body, and h₁and h₂ are as shown in Fig. 2.

In Fig. 3a and b, the normalized groundwater head and normalized groundwater discharge per unit width at location $x=L/\mathrm{2}$ through time under different fractional power values are shown. Meanwhile, Fig. 3c shows the normalized groundwater head at the instance in time t=40 000 min as a function of location throughout the body of the dam and the analytical solution to the standard governing equation of unconfined groundwater flow when ${\mathit{\beta }}_x=\mathit{\alpha }=\mathrm{1}$ at the steady state. The considered fractional derivative powers in space and time are ${\mathit{\beta }}_x=\mathit{\alpha }=\mathrm{0.7}$, 0.8, 0.9, 1.0. As can be seen from Fig. 3a, the hydraulic head recession in time slows down with the decrease of β_x=α from 1. The hydraulic heads in Fig. 3a have heavier tails as orders of time and space fractional derivative powers decrease from 1 towards 0.7. Furthermore, normalized groundwater discharge per unit width in Fig. 3b goes to 1 at a slower rate as fractional derivative powers decrease from 1 towards 0.7. Meanwhile, Fig. 3c shows that the numerical solution to the governing fractional equation at ${\mathit{\beta }}_x=\mathit{\alpha }=\mathrm{1.0}$ and at a very long time after the initial condition perfectly matches the steady-state analytical solution, Eq. (24), to the standard Eq. (23) with the specified initial and boundary conditions.

Figure 4 The sketch numerical application 2; the downstream groundwater head is fixed at 11 m and the initial upstream groundwater head is 16 m.

Figure 5 Results for numerical application 2: (a) the normalized initial groundwater head in the unconfined aquifer and the normalized groundwater head at time t=60 000 min through length of the aquifer under different fractional derivative powers; (b) the normalized groundwater head at $x=L/\mathrm{2}$ through time under different fractional derivative powers; (c) the normalized groundwater discharge per unit width at $x=L/\mathrm{2}$ through time under different fractional derivative powers; t is time and the simulation time T is 60 000 min.

The second application (numerical application 2) deals with a transient unconfined groundwater flow from a hillslope toward a stream (Fig. 4). The upstream boundary plane separates the region of flow from the adjacent hillslope that feeds the adjacent tributary system; therefore $\frac{\partial h}{\partial x}=\mathrm{0}$ (Freeze, 1978) at x=0. The normalized initial groundwater head in the unconfined aquifer and the normalized groundwater head at time t=60 000 min through the length of the aquifer under different fractional derivative powers are shown in Fig. 5a. The normalized groundwater head and normalized groundwater discharge per unit width at x=L/2 through time under different fractional derivative powers are demonstrated in Fig. 5b and c. As can be seen from Fig. 5b and c, the hydraulic head and groundwater discharge recession in time slows down with the decrease of β_x=α from 1. The hydraulic heads and groundwater discharges in Fig. 5b and c have heavier tails as orders of time and space fractional derivative powers decrease from 1 towards 0.7.

6 Discussion

From the standard governing Eq. (23) of unconfined groundwater flow in integer time–space the saturated hydraulic conductivity may be interpreted as a diffusion coefficient for the nonlinear diffusion of groundwater in an unconfined aquifer. The basic difference between confined and unconfined groundwater flow is that the former can be interpreted as a linear diffusion of groundwater while the latter is a nonlinear diffusion of groundwater within an aquifer. Similar to saturated hydraulic conductivities in Eq. (26) in Kavvas et al. (2017a) for the fractional confined aquifer groundwater flow, the saturated hydraulic conductivities in Eq. (21), which governs fractional unconfined aquifer groundwater flow, are modulated by the ratios of fractional time to fractional space, $\frac{t^{\mathrm{1}-\mathit{\alpha }}}{x_i^{\mathrm{1}-{\mathit{\beta }}_{x_i}}}$, i=1, 2. In other words, the confined and unconfined groundwater diffusions in fractional time–space are modulated by the above fractional time–space ratios.

Numerical applications demonstrated that as the powers of the space and time fractional derivatives decrease from 1, the recession rate of the nondimensional groundwater hydraulic head slows down compared to the case by the conventional governing equation (i.e., with integer-order derivatives). This behavior also indicates the modulation of the nonlinear diffusion of the groundwater by the fractional powers of time and space.

As mentioned in the Introduction section, unconfined groundwater flow is the fundamental component of the watershed runoff baseflow since it is the fundamental contributor to the streamflow network within a watershed during dry periods. As such, the behavior of unconfined groundwater flow is key to the physically based understanding of the long memory in watershed runoff. As seen from the numerical applications in Figs. 3 and 5, the powers of the fractional derivatives in time and space can modulate the speed of the groundwater discharge evolution. Hence, they can modulate the memory of the unconfined aquifer flow, which, in turn, can modulate the memory of the watershed baseflow. Meanwhile, the Caputo derivative, as defined in its special form $D_{\mathrm{0}}^{{\mathit{\beta }}_{x_i}}f\left(x_i\right)$ in space in this study, was shown by Kavvas and Ercan (2016) and Ercan and Kavvas (2017) to be a nonlocal quantity where the effect of the boundary conditions on the groundwater flow within the flow domain can have long spatial memories with the decrease in the powers of the spatial fractional derivatives from unity. Similarly, it was shown by Kavvas et al. (2017a) that the Caputo derivative in time, $D_{\mathrm{0}}^{\mathit{\alpha }}f\left(t\right)$, as defined in this study, is nonlocal in time and can carry the effect of initial conditions on the groundwater flow for long time periods as the power in the time fractional derivative decreases from 1. Therefore, the fractional governing equation of unconfined groundwater flow in fractional time and multi-fractional space has the potential to describe the long-memory characteristics of baseflow within a watershed. This important topic shall be explored in the near future.

7 Conclusion

A dimensionally consistent fractional governing equation of transient unconfined aquifer groundwater flow was derived within a fractional differentiation framework. After developing a fractional continuity equation, a previously developed dimensionally consistent equation for water flux in unsaturated/saturated porous media was combined with the Dupuit approximation to obtain an equation for groundwater motion in multi-fractional space in unconfined aquifers. By combining the fractional continuity and motion equations, the governing equation of transient unconfined aquifer groundwater flow in a multi-fractional medium in fractional time was obtained. To demonstrate the skills of the proposed fractional governing equation of unconfined aquifer groundwater flow, two numerical applications were performed. As demonstrated in the numerical application results, the orders of the fractional space and time derivatives modulate the speed of groundwater discharge and groundwater flow evolution, slowing the process with the decrease in the powers of the fractional derivatives from 1. It is also shown that the proposed dimensionally consistent fractional governing equations approach the corresponding conventional equations as the fractional orders of the derivatives go to 1.

Appendix A: Numerical solution for the one-dimensional case

One-dimensional time–space fractional groundwater flow in the unconfined aquifer with no recharge or leakage can be written as

$$\begin{array}{}\text{(A1)}& \begin{array}{rl}& S_y{\displaystyle \frac{{\partial }^{\mathit{\alpha }}h}{{\left(\partial t\right)}^{\mathit{\alpha }}}}={\displaystyle \frac{\mathrm{\Gamma }(\mathrm{2}-\mathit{\beta })}{x^{\mathrm{1}-\mathit{\beta }}}}{\left({\displaystyle \frac{\partial }{\partial x}}\right)}^{\mathit{\beta }}\\ & \left(K_{\mathrm{s}}\left(\stackrel{\mathrm{‾}}{X}\right){\displaystyle \frac{t^{\mathrm{1}-\mathit{\alpha }}}{x^{\mathrm{1}-\mathit{\beta }}}}{\displaystyle \frac{\mathrm{\Gamma }(\mathrm{2}-\mathit{\beta })}{\mathrm{\Gamma }(\mathrm{2}-\mathit{\alpha })}}h{\displaystyle \frac{{\partial }^{\mathit{\beta }}h}{{\left(\partial x\right)}^{\mathit{\beta }}}}\right).\end{array}\end{array}$$

The fractional time and space derivatives are estimated in the same manner as seen in Tu et al. (2018), where the Caputo fractional space and time derivatives in the fractional governing equation are estimated by the numerical algorithm in Odibat (2009) and the algorithm reported by Murio (2008), respectively. The Caputo fractional space derivative ${\left.D_x^{\mathit{\beta }}f\left(x\right)\right|}_{x=L}{\left.D_x^{\mathit{\beta }}g\left(x\right)\right|}_{x=L}$ at the location L for $m-\mathrm{1}<\mathit{\beta }\le m$ (m∈N) for a given space interval [0,L] is estimated as

$$\begin{array}{}\text{(A2)}& \begin{array}{rl}& {\left.D_x^{\mathit{\beta }}g\left(x\right)\right|}_{x=L}\approx {\displaystyle \frac{\mathrm{\Delta }{L}^{m-\mathit{\beta }}}{\mathrm{\Gamma }\left(m+\mathrm{2}-\mathit{\beta }\right)}}\mathit{\left\{}\right[{\left(N-\mathrm{1}\right)}^{m-\mathit{\beta }+\mathrm{1}}\\ & -\left(N-m+\mathit{\beta }-\mathrm{1}\right)\mathrm{\Delta }{L}^{m-\mathit{\beta }}]g^{\left(m\right)}\left(\mathrm{0}\right)+g^{\left(m\right)}\left(L\right)\\ & +\sum _{i=\mathrm{1}}^{N-\mathrm{1}}[{\left(N-i+\mathrm{1}\right)}^{m-\mathit{\beta }+\mathrm{1}}-\mathrm{2}{\left(N-i\right)}^{m-\mathit{\beta }+\mathrm{1}}\\ & +{\left(N-i-\mathrm{1}\right)}^{m-\mathit{\beta }+\mathrm{1}}\left]g^{\left(m\right)}\left(l_i\right)\mathit{\right\}},\end{array}\end{array}$$

where N is the number of equally spaced subintervals on [0,L]; the subinterval length is $\mathrm{\Delta }L=L/N$ and l_i=iΔL, for $i=\mathrm{0},\mathrm{1},\mathrm{2},\mathrm{\dots },N$.

The Caputo fractional time derivative ${\left.D_t^{\mathit{\alpha }}g\left(x,t\right)\right|}_{x=l_i,t=t_n}$ for $\mathrm{0}<\mathit{\alpha }\le \mathrm{1}$ for a given time interval [0,T], which is divided into M equal subintervals with a time window of $\mathrm{\Delta }t=T/M$ by using the nodes t_n=nΔt, $n=\mathrm{0},\mathrm{1},\mathrm{2},\mathrm{\dots },M$, can be approximated as

$$\begin{array}{}\text{(A3)}& \begin{array}{rl}& D_t^{\mathit{\alpha }}{g}_i^n={\displaystyle \frac{\mathrm{\Delta }{t}^{-\mathit{\alpha }}}{\mathrm{\Gamma }\left(\mathrm{2}-\mathit{\alpha }\right)}}\sum _{k=\mathrm{1}}^n\left[k^{\mathrm{1}-\mathit{\alpha }}-{\left(k-\mathrm{1}\right)}^{\mathrm{1}-\mathit{\alpha }}\right]\\ & \left(g_i^{n-k+\mathrm{1}}-g_i^{n-k}\right).\end{array}\end{array}$$

Then the 1-D governing equation in fractional time and space for Cartesian groundwater flow in an unconfined aquifer can be discretized as

• for n=1,

  $$\begin{array}{}\text{(A4)}& h_i^n=h_i^{n-\mathrm{1}}+{\displaystyle \frac{t_n^{\mathrm{1}-\mathit{\alpha }}}{S_y\mathrm{\Delta }{t}^{-\mathit{\alpha }}}}{\displaystyle \frac{\mathrm{\Gamma }\left(\mathrm{2}-\mathit{\beta }\right)}{l_i^{\mathrm{1}-\mathit{\beta }}}}{G}_i^{n-\mathrm{1}}.\end{array}$$

• for n≥2,

  $$\begin{array}{}\text{(A5)}& \begin{array}{rl}& h_i^n=h_i^{n-\mathrm{1}}+{\displaystyle \frac{t_n^{\mathrm{1}-\mathit{\alpha }}}{S_y\mathrm{\Delta }{t}^{-\mathit{\alpha }}}}{\displaystyle \frac{\mathrm{\Gamma }\left(\mathrm{2}-\mathit{\beta }\right)}{l_i^{\mathrm{1}-\mathit{\beta }}}}{G}_i^{n-\mathrm{1}}\\ & -\sum _{k=\mathrm{2}}^n\left[k^{\mathrm{1}-\mathit{\alpha }}-{\left(k-\mathrm{1}\right)}^{\mathrm{1}-\mathit{\alpha }}\right]\left(h_i^{n-k+\mathrm{1}}-h_i^{n-k}\right),\end{array}\end{array}$$

where $G={\left(\frac{\partial }{\partial x}\right)}^{\mathit{\beta }}\left[K_{\mathrm{s}}\left(\stackrel{\mathrm{‾}}{X}\right)\frac{\mathrm{\Gamma }\left(\mathrm{2}-\mathit{\beta }\right)}{x^{\mathrm{1}-\mathit{\beta }}}h\frac{{\partial }^{\mathit{\beta }}h}{\partial x^{\mathit{\beta }}}\right]$, and the space and time fractional derivatives in G are estimated as in Eqs. (A2) and (A3).