Exact and Approximate Solutions of Fractional Partial Differential Equations for Water Movement in Soils

This paper presents solutions of the fractional partial differential equation (fPDE) for analysing water movement in soils. The fPDE explains processes equivalent to the concept of symmetrical fractional derivatives (SFDs) which have two components: the forward fractional derivative (FFD) and backward fractional derivative (BFD) of water movement in soils with the BFD representing the micro-scale backwater effect in porous media. The distributed-order time-space fPDE represents water movement in both swelling and non-swelling soils with mobile and immobile zones with the backwater effect operating at two time scales in large and small pores. The concept of flux-concentration relation is now updated to account for the relative fractional flux of water movement in soils.


Introduction
Water movement on the surface of the soil, into and out of the subsurface terrain, is a process which interfaces the atmosphere and the subsurface through physical infiltration and evaporation as well as transpiration by plants.Understanding of the movements of water into and out of the soils in different saturated states and through aquifers requires knowledge which connects soil science, groundwater hydrology, fluid mechanics and other related topics.
In this paper the water movement in soils of finite depths is investigated by presenting analytical solutions of the fractional partial differential equation (fPDE) of water movement.The fPDE discussed in this paper is one special form of the fPDEs presented earlier [1] following the extensive analysis and parameterisation of the flow processes in soils.The fPDEs incorporate both the material coordinate for water movement in swelling-shrinking soils and the Cartesian coordinates in non-swelling soils.For the derivation of the fPDE for water flow in soils, the connection has been established between the mass-time and space-time fPDEs and the continuous-time random walk (CTRW) theory.We also clarify that the fPDEs presented earlier based on the CTRW concept [2,3] incorporate the forward fractional derivative (FFD) and backward fractional derivative (BFD) introduced by Bochner [4] and demonstrated by Saichev and Zaslavsky [5] and Zaslavsky [6].
This paper further investigates water movement in non-swelling soils by presenting solutions with known boundary conditions for a finite soil profile.The solutions are also applicable to a material coordinate for swelling soils where the finite mass is used instead of the finite depth.These analyses are based on models formulated for water flow in soils with or without an immobile zone, and consider the following major issues in soil water movement: (1) Swelling and non-swelling properties of soils; Hydrology 2017, 4, 8; doi:10.3390/hydrology4010008www.mdpi.com/journal/hydrology(2) Mobile and immobile zones in soils are optional.The mobile-immobile concept, which was first proposed by Barenblatt et al. in 1960 [7], is taken into account in both the mass-time and space-time fPDEs; (3) The connection between the mass-time fPDEs and CTRW theory for soils is further discussed, and (4) Space-time scaling laws [6] are briefly discussed for the fPDEs in the context of mass-time and space-time fractional derivatives.
The CTRW is a further development from the concept of the random walk which was first mathematically illustrated by Crofton [8] in terms of random flights as an alternative term, and later used by Pearson [9].This concept was successfully employed to model solute movement in porous media by Saffman [10].The CTRW theory for particle transport is understood in the framework of the classical renewal theory [11].The mechanics of flow in porous media can be ideally described by the CTRW theory, and the phenomena described by the CTRW concept is anomalous which is a generalisation of the transport processes with the classic diffusion as a special case [3].
The mass-time fPDE is an extension of the space-time fPDE, which is a consequence of the long-time limit or asymptotic result of the two probability density functions of the CTRW model with power-law waiting times and power-law jumps [2,3,6,12].These distributed-order fPDEs with convection included form the mass-time fractional advection-diffusion equation (fADE) and the space-time fADE.This connection offers another option for deriving the fPDE-based models without resorting to the traditional mass balance method for their derivation.
The symmetrical fractional derivatives (SFDs) include both the backward fractional derivative (BFD) and the forward fractional derivatives (FFD) in space.The BFD is capable of accounting for the backwater effect on particle and water parcel motion at a micro-scale, compared to the well-known large-scale backwater effect in hydraulics.
The multi-term fractional derivatives in time used in the fPDE by Hadid and Luchko [13] and Jiang et al. [14] are an ideal way to model flow and particle motion in fractal media with an unlimited level of micro-structures, which is a key characteristic of fractal media.In this paper a two-term Mittag-Leffler function is used in the solution to account for the mobile-immobile zone as an approximation to the multinomial Mittag-Leffler function.The examples of mobile-immobile models include those for solute transport in saturated media by Schumer et al. [15], water flow in unsaturated soils by Su [1,16] and saturated aquifers [17].
By combining the BFD and the widely-used two-term mobile-immobile model, the fPDE is expected to embrace more processes and provide informative insights into processes of water flow in soils.

The CTRW Theory and Its Connection with the fPDE for Water Movement in Soils with Wandering Processes
In a comprehensive development [1], we have introduced a set of fPDEs for water flow in both swelling and non-swelling soils.For swelling soils, with the material diffusivity given by and the hydraulic conductivity by where D 0 , b, K 0 and k are constants, the equation governing the movement of water in swelling soils is given as b 1 where q is the Darcy's flux, λ = µ + 1 and m is the material coordinate defined by [18,19] where m is the conventional space coordinate in a vertical direction, and e is the void ratio given by e = θ/(1 − θ) with θ the moisture ratio defined as θ = θ l /θ s , and θ l and θ s being the volume fractions of liquid and solid, respectively; b 1 and b 2 are the relative porosities in immobile and mobile zones, respectively, i.e., b 1 = ϕ im ϕ and b 2 = ϕ m ϕ with ϕ im , ϕ m and ϕ being the porosities in the immobile and mobile zones, and total porosity, respectively; β 1 and β 2 are the orders of fractional derivatives for immobile and mobile zones, respectively; γ n is the particle specific gravity, and α is the gradient (or slope) of the shrinkage curve, which is a ratio on the graph of the specific volume, v, versus water content or moisture ratio, θ.
The equation for water movement in one dimension in non-swelling soils is of the form where the power functions similar to Equations ( 1) and ( 2) are used for the diffusivity and hydraulic conductivity for non-swelling soils, respectively.Multidimensional fPDEs have also been presented for water flow in non-swelling soils [1].Equations ( 3) and ( 5) were derived based on the CTRW theory, and the notations used for their derivation are identical to those in Gorenflo and Mainardi [2,20] and Gorenflo et al. [3].Here we should emphasise that the SFDs in [2,3,20] are identical terminologies in Saichev and Zaslavsky [5] who derived fPDEs by defining the concept of "wandering process" with the BFD and FFD illustrated by Bochner [4].See Appendix A in this paper for more details.
For water movement in two directions in swelling soils, the BFD and FFD are denoted using the sign |m| and Equation (3) can be more generally written as which is the mass-time distributed-order fADE for water movement in swelling soils with the material coordinate being dependent on the moisture ratio.In other words, the term ∂ λ θ ∂|m| λ in Equation ( 6) includes both the forward and backward fractional components if the particle or water parcel motion is regarded as a wandering process [5].
Equation ( 6) is based on the CTRW theory with distributed-order time fractional derivatives resulting from the two time-scaling property and convection due to a shift jump size distribution in the CTRW theory [21].The connection between CTRW and the distributed-order fPDE has been established [22], which is also applicable to Equation (6) and is an extension of the CTRW theory and fPDE to swelling soils.
The parameters and the material coordinate in Equation ( 6) characterise the different flow patterns in swelling soils with mobile and immobile zones which have a variable material diffusivity, D m (θ) and hydraulic conductivity, K m (θ) as functions of the moisture ratio.
Similar to Equation ( 6) for swelling soils, Equation ( 5) for non-swelling soils can be more generally written as where z is the usual physical coordinate.The large-scale backwater effect is well known in hydraulics.With the BFD concept in space, the backwater effect is now extended to flow at a micro-scale, and the distributed-order time fPDE for water movement with the forward and backward fractional derivatives accounts for the anomalous water movement with the backwater effect at two time-scales in soils.
Similar to Equations ( 8) and (A6) and (A7) in the Appendix A, the signs of ∂ λ θ ∂|z| λ and D λ 0 θ for non-swelling soils are identical, then ∂ λ θ ∂|z| λ = D λ 0 θ is used throughout the text as the synonyms of the SFDs.
Compared to the classic advection-diffusion equation (ADE), the space-and time-fractional models are shown to better represent the first and second moments of solute transport at early times and the tails of tracer plumes [21].As water is the carrier of the solutes in soils, the fractional model is then expected to perform better than its classic counterpart for water flow in soils, particularly at the early and final stages of infiltration processes at the soil surface.

Solutions of the Distributed-Order Fractional Partial Differential Equations Incorporating Forward and Backward Motion of Water Flow in Soils
Here we present solutions of Equation ( 8), the space-time fPDE for water flow in non-swelling soils, which are also applicable to Equation ( 6) with the only difference being m replacing z.The solutions presented here are modifications of those given by Jiang et al. [14] with only two terms of time fractional derivatives retained to account for the two-zone model (the mobile-immobile or large-small porosity model) for flow in porous media.The solutions are subject to the following initial condition (IC) and boundary conditions (BCs), where L is the depth of the soil profile.The moisture ratios on the surface and the depth z = L are, respectively, θ(0, t) and θ(L, t).The complete solution of Equation ( 8) subject to the above conditions is Equation (B12) in Appendix B, and the approximate solution by retaining the first term of the full solution is as follows, where is from Equation (B6) for j = 0, 1 and i = 0, 1, and is the temporal component of the moisture ratio, a 1 and a 2 are from Equations (B9) and (B10), respectively.Equation ( 13) is identical to Equation (B14) through the identity, 1 − cos 2πz L = 2 sin 2 πz L , in trigonometry [23].
The approximate solution in Equation ( 13) can be used to describe the movement of water in non-swelling soils in terms of the moisture ratio, θ(z, t) and in the flux-concentration relation to be discussed below.
At z = 0, Equation ( 13) becomes θ(0, t), and at z = L, it becomes θ(z, t) = θ(L, t).Clearly, for t = 0 Equation ( 13) yields the initial moisture ratio profile as To illustrate the effects of the two orders of fractional derivatives on F(t) in Equation ( 14), Figures 1  and 2 are generated to examine the effects of each order of the fractional derivatives on the moisture distribution patterns.
Hydrology 2017, 4, 8 5 of 14 The approximate solution in Equation ( 13) can be used to describe the movement of water in non-swelling soils in terms of the moisture ratio, ( , ) z t θ and in the flux-concentration relation to be discussed below.At 0 z = , Equation ( 13) becomes (0, ) t θ , and at L z = , it becomes ( , ) ( , ) . Clearly, for 0 t = Equation ( 13) yields the initial moisture ratio profile as ( ) To illustrate the effects of the two orders of fractional derivatives on ( ) in Equation ( 14), Figures 1 and 2 are generated to examine the effects of each order of the fractional derivatives on the moisture distribution patterns.The patterns of moisture distributions in Figure 1 show that for a fixed small pore structure represented by 1 β (Figure 1A), the flow of soil water is faster as 2 β increases, and that for a fixed large pore structure represented by 2 β , the flow slows down as the value of 1 β representing small pores increases.These numerical values and flow patterns are intuitive.It is clearly seen in Figure 2 that is a monotonically decreasing function of time and the increase in both values of the two orders of fractional derivatives enhances the steepness of the curves implying the faster movement of soil water in the profile.It is clear that compared to the patterns in Figure 1 the increase in the flow rate by increasing 2 β overtakes the role of 1 β leading to the overall Figure 2. The facilitating role of β 2 overtaking β 1 in increasing the flow rate when both orders increase.The illustration of Equation ( 14) with pairs of β 2 and β 1 (the joint parameter K 1 = 0.01).
The patterns of moisture distributions in Figure 1 show that for a fixed small pore structure represented by β 1 (Figure 1A), the flow of soil water is faster as β 2 increases, and that for a fixed large pore structure represented by β 2 , the flow slows down as the value of β 1 representing small pores increases.These numerical values and flow patterns are intuitive.
It is clearly seen in Figure 2 that F(t) is a monotonically decreasing function of time and the increase in both values of the two orders of fractional derivatives enhances the steepness of the curves implying the faster movement of soil water in the profile.It is clear that compared to the patterns in Figure 1 the increase in the flow rate by increasing β 2 overtakes the role of β 1 leading to the overall increase in the flow rate.

Flux-Concentration Relations at Different Depths
By definition in terms of mass conservation, the flux of moisture is given as Differentiating Equation ( 13) and using it with Equation ( 16) gives where For simplicity denoting F t = F(t), and for very long time, F t → 0 , Equation ( 17) approaches and on the surface, z = 0, the flux given by Equation ( 17) for F t → 0 is q 0 , q = K(θ(0, t)) − ξD 0 (20) The ratio is termed as the flux-concentration relation (FCR) [24,25], which is simply the relative flux at different depths.
For very long time, F t → 0 and from Equations ( 19) and ( 20), the FCR in Equation ( 21) is which mean that the hydraulic conductivity and diffusion coefficient play an important role in determining FCR at a large time.
In practice, different functions have been used for the unsaturated hydraulic conductivity and the diffusion coefficient, and Equation have different shapes depending on the types of the functions for these two parameters.

Conclusions and Discussions
The purpose of this paper is to investigate the movement of water in soils using a newly developed fPDE [1].The results of this investigation in this paper include the (1) Analytical solutions and their approximations are presented for a distributed-order mass-time and space-time fPDE for water movement in soils of finite depths.We limit our analysis to the model with two-term fractional distributed orders in the fPDE to account for the large-small pores (or mobile-immobile zones), which is widely used in soil science and hydrology.The solutions derived for non-swelling soils are identical for solutions of water movement in swelling soils by changing z to m with relevant parameters included.(2) It is shown that the fPDE results from the asymptotic or long-time approximation of the CTRW model with power laws as the two transitional probability distribution functions for the length of jumps and waiting time intervals.The symmetrical fractional derivatives include the backward and forward fractional derivatives with the former representing the wandering process of soil water movement.The backward fractional derivative accounts for the backwater effect at a micro-scale which is a counterpart of the well-known large-scale backwater effect in hydraulics.With these properties the symmetrical fractional derivatives are ideal for describing stochastic movement of water in porous media.The presentations here and earlier [1] include the convection term due to gravity, ∂|z| for non-swelling soils and ∂K(θ) ∂|m| for swelling soils.Due to the fact that the diffusivity, D(θ), and the hydraulic conductivity, K(θ), are related functions [1], the fractional connections between D(θ) and K(θ) are more complex which need to be explored.In the present case, a constant diffusivity and hydraulic conductivity are used and, with functions for D(θ) and K(θ), the fPDE with a fractional term ∂|m| η with 0 < η < 1 needs to be explored.For a general background to fractional calculus, its applications, and the methods for solutions as well as related topics, the reader is referred to other sources [2,3,6,12,15,21,[26][27][28][29].While more details of the fPDE-based models are being explored and refined, there is a need to develop methods for estimating the parameters in the models for unsaturated flow using data from both laboratory and field measurements, and the method suggested in [30] for deriving the conductivity from the fractional boundary-value problem is one option.Detailed methods for estimating these parameters are beyond the scope of this paper.
with ∂ λ θ ∂|z| λ and ∂ η θ ∂|z| η being Riesz space fractional derivatives (RSFD).For material movement due to diffusion and convection in a homogeneous media, the RSFD is defined as [5], The forward motion is represented by the forward fractional derivatives, D + ∂ λ + θ ∂z λ + and V + ∂ η + θ ∂z η + , while the backward fractional derivatives, represent the backward motion of the water in the soil.These terms with different parameters, D + , D − , V + and V − manifest the anisotropic properties of the media.
For a simplified situation when the media is isotropic, where Saichev and Zaslavsky [5] defined the fractional diffusion coefficient, D, which can be extended to the fractional velocity, V, Then the symmetric fractional derivatives (SFD) are defined as Note that the SFD above is defined for homogeneous and anisotropic media.In the above derivations, θ = θ(z, t).
Two important properties of SFDs are that ∂ λ θ ∂|x| λ = 0 is defined for x > 0 and ∂ λ θ ∂|x| λ = 0 for x < 0, and on the contrary, ∂ λ θ ∂|−x| λ = 0 is defined for x < 0 and ∂ λ θ ∂|−x| λ = 0 for x > 0 [5].Umarov and Gorenflo [31] show that the pseudo-differential operator for fractional derivatives are identical to the fractional power of the fractional Laplace operator, For one-dimensional fractional derivatives of the known function θ, Equation (A8) gives For a finite domain [0, L; 0, T] with the homogeneous boundary conditions θ(0, t) = θ(L, t) = 0, Jiang et al. (2012, Equation ( 9)) show that the following equality holds in one dimension where 0 D λ x θ and x D λ L θ are the left-sided and right-sided fractional derivatives, respectively, and Comparing Equations (A9) and (A10) shows that the signs for the SFDs ∂ λ θ ∂|x| λ and D λ 0 θ are identical, then D λ 0 θ is used throughout the text as the synonyms of the symmetric fractional derivatives.In Equations (A4) and (A5), the dimensions of D and V depend on the value of λ and η.Their dimensions in exact solutions should remain the same as in the main model equations.However, due to approximations, their dimensions could change depending on how the approximations are made.For further discussions of the symmetric fractional calculus, the reader is referred to Saichev and Zaslavsky [5] and Umarov and Gorenflo [31].8) Subject to the Constant Boundary Conditions in Equations ( 10) to (12) The solution of Equation ( 8) with b = 0 and k = 1 for water movement in non-swelling soils can be derived by modifying the solution of Jiang et al. [14] as follows,

Appendix B. The Solution of Equation (
where θ n (z) is given as (Equation (38) in [14]) and where G 2 β (t) is the two-term Mittag-Leffler function (2MLF) given by the following expression (Equation (19) in [32]) with Then, Equation (B3) can be written as The solution of Equation ( 6) for swelling soils is similar to Equation (B1) for non-swelling soils except for m replacing z in the solution and with the joint parameter, K n , in Equation (B5) being The solution in Equation (B1) is not complete in that an integral in Equation (B2) appears uncompleted.In the following section, we assign specific values to the IC and BCs so that the integration in Equation (B2) can be completed and explicit solutions of Equation (B1) expressed in algebraic forms.

Appendix B.1. Particular Solutions with Constant Known Initial Condition and Boundary Conditions
The solutions in Equation (B1) are to be completed by specifying particular values for the IC, θ(z, 0), and BCs, θ(0, 0) and θ(L, 0).
Let us specify that the initial moisture ratio in the soil, θ(z, 0), and their values at the two boundaries, θ(0, 0) and θ(L, 0), are all constant, then we rewrite Equation (B2) as and Equation (B8) can be integrated (see Equation (1), p. 227 in [33]), combined with the upper boundary condition which determines the constant of integration to be zero, to yield The solution in Equation (B1) can be now written as It is clear that Equation (B12) yields the following moisture values at the boundaries θ(z, t) = θ(0, t), at z = 0, t > 0 and θ(z, t) = θ(L, t) at z = L, t > 0.
The initial profile at t = 0 is The series solution in Equation (B12) can be achieved for practical use by retaining only limited terms.One approximation is for n = 1, which yields where with the aid of the properties of trigonometric functions (see [23], p. 165).
The solution in Equation (B12) and its approximation in (B14) are also applicable to swelling soils when m replaces z in the solution with K n in (B3) given by Equation (B7).
In the solutions presented in Appendix B, the 2MLF can be approximated by retaining limited terms.Here we retain only two leading terms in the 2MLF in Equation (B4) with j = 0, 1 and i = 0, 1 resulting in where K 1 is given by Equation (B15).Then, the solution of Equation ( 8) subject to the condition in Equations ( 10)-( 12) is now given by θ(z, t) = θ(0, t) + [θ(L, t) − θ(0, t)]z L + ∞ ∑ n=1 θ n (z)F n (t) sin nπz L (B21) with θ n (z) given by Equation (B19).The solution of Equation ( 6) for swelling soils is similar to Equation (B21) except for m replacing z and Equation (B7) for K n .

Figure 1 .
Figure 1.Moisture distributions affected by the orders of fractional derivatives in Equation (14).(A): Fixed 1 β with variable 2 β indicating the facilitating role of 2 β in large pores; (B): fixed 2 β with variable 1 β indicating the trapping role of 1 β in small pores.In both cases, the joint parameter is 1 0.01 K = .

( 3 )
The flux-concentration relation is shown to include fractional parameters in the fPDE, and a large-time asymptote is given.(4) The temporal component of the solutions are illustrated to examine the effect of the model parameters β 2 and β 1 on flow processes which are shown to explain realistic physical processes.