Fractional Vertical Inﬁltration

: The inﬁltration phenomena has been studied by several authors for decades, and numerical and approximate results have been shown through the asymptotic solution in short and long times. In particular, it is worth highlighting the works of Philip and Parlange, who used time and volumetric content as independent variables and space as a dependent variable, and found the solution as a power series in t 1/2 that is valid for short times. However, several studies show that these models are not applicable to anomalous ﬂows, in which case the application of fractional calculus is needed. In this work, a fractional time derivative of a Caputo type is applied to model anomalous inﬁltration phenomena. Fractional horizontal inﬁltration phenomena are studied, and the fractional Boltzmann transform is deﬁned. To study fractional vertical inﬁltration phenomena, the asymptotic behavior is described for short and long times considering an arbitrary diffusivity and hydraulic conductivity. Finally, considering a constant ﬂux-dependent relation and a relation between diffusivity and hydraulic conductivity, a fractional cumulative inﬁltration model applicable to various types of soil is built; its solution is expressed as a power series in t ν /2 , where ν ∈ ( 0,2 ) is the order of the fractional derivative. The results show the effect of superdiffusive and subdiffusive ﬂows in different types of soil.


Introduction
More than a century ago, Green and Ampt [1] laid the bases for describing the infiltration process; to this day, there are still open problems in the area, where existing models in the literature cannot explain laboratory results [2]. As one of the pioneers in the area, Philip [3] describes the infiltration phenomenon as a Fokker-Planck type equation. In a semi-infinite domain, homogeneous, its form is where θ is the volumetric water content, D(θ) is the hydraulic diffusivity, K(θ) is the hydraulic conductivity, θ 0 is the initial volumetric water content, θ s is the volumetric water content at saturation, t is the time, and z is the spatial variable taken positive downward.
In particular, when considering horizontal infiltration: with initial conditions by introducing the Boltzmann transform [4]. On the other hand, when considering vertical infiltration, ∂θ ∂t Philip [5] applied the implicit function theorem and rewrote Equation (6) as where q = q(θ, t) is the Darcy flow; Philip proved that, for short times, the solution of Equation (8) is sought in the form z(θ, t) = φ(θ)t 1/2 + χ(θ)t + Ψ(θ)t 3/2 + ω(θ)t 2 + · · · (9) where each term function of θ in the series is defined from integro-differential equations whose expression considering fractional calculus will be shown in Section 4. In particular, the first term in the series (9), φ(θ), matches the solution of the horizontal infiltration problem, Equation (4). However, Philip did not work with the radius of convergence of the series, a problem that remains open to date, but suggests that the range of useful convergence is t grav , where it is obtained from dimensional analysis that gravity effects are as great as capillarity effects, namely, where K s = K(θ s ) is the hydraulic conductivity at saturation, K 0 = K(θ 0 ) is the initial hydraulic conductivity, and S is the sorptivity defined by where φ(θ) is the Boltzmann transform defined by the differential Equation (4). Likewise, Philip [5] proves that, for long times, the solution of Equation (8) has the form where ζ(θ) is a function of volumetric water content and can be expressed in terms of diffusivity and hydraulic conductivity as where θ a = θ s − ∆θ and ∆θ is a definite small positive quality; the calculus of ζ, Equation (13), should be done carefully since, in general, ζ has singularities at θ = θ s , θ 0 .
Philip showed that, when water is applied to uniform soils, it is sufficient to retain only two terms of the series to obtain an adequate description of the vertical infiltration. Talsma and Parlange [6] proved that sorptivity S and conductivity at saturation K s were the soil parameters necessary for the two terms of the Philip series; in addition, the onedimensional infiltration equation is more sensitive to soil heterogeneity than its equivalents in two or three dimensions.
Later, Parlange et al. [7] showed that adding a third parameter to the series is sufficient to describe cumulative infiltration accurately, where this extra parameter is representative of the type of soil in which water infiltration takes place.
Haverkamp et al. [8] recognize that the infiltration process can be described either by solving the Richards equation or by considering cumulative infiltration. They chose the second option. Since previous models used highly time-dependent parameters [9], they developed an equation for infiltration subject to head boundary conditions considering time-independent parameters, which takes into account the possibility of an infinite diffusivity near saturation.
Fuentes et al. [10] analyzed the constraints on different fitting parameters used in a water retention equation and hydraulic conductivity using the infiltration equations to describe several soils types. Later, Fuentes et al. [11] calculated the analytical solution for vertical infiltration considering a constant diffusivity and a hydraulic conductivity resulting from a linear combination of linear and quadratic functions of volumetric water content. They showed that their model can also be applied to soils that do not necessarily meet the hypotheses for his deduction. On the other hand, using a numerical approximation, Saucedo et al. [12] verified the contact time hypothesis to describe water transfer in melgas irrigation by coupling the Saint-Venant and Richards equations.
Over the years, the Philip power series has been applied in different contexts, such as falling head ponded infiltration [13], variable head ponded infiltration on flat surfaces [14], and variable head ponded infiltration on sloping surfaces [15], to name a few. In particular, in the last work, they showed that sorptivity is independent of the slope angle and that cumulative infiltration normal to the slope decreases with increasing slope angles.
Recently, cumulative infiltration was used to describe the effect of pressure head and soil bulk density on moistube irrigation [16]. It was concluded that infiltration index α is negatively correlated with pressure head but positively with soil bulk density, while the infiltration coefficient K was positively correlated with the pressure head but negatively with the soil bulk density.
Likewise, water infiltration has been used to calculate infiltrated depth evolution and humidity profiles by coupling the Saint-Venant and Richards equations with furrow boundary conditions [17]. Considering this model, the relationship between the optimal irrigation flow and the length of the border for different types of soil was found [18].
These examples show the wide range of applications that the Philip power series has for water infiltration in different types of soil; in this sense, it goes without saying how useful infiltration models are in developing an effective design and evaluating surface irrigation systems, whether they are physical (such as Philip's [5] and Richards' [19] models), semi-empirical (such as Holtan's [20] model) or empirical models (such as Sihag's [21] model) [22].
When studying non-heterogeneous soils, the transport phenomenon in soils is recognized as being fundamentally stochastic [23,24]. Moreover, to apply equations, e.g., Richards or Philip power series, it is necessary to assume that water moves in a Brownian motion [25]; however, several studies have shown that this assumption is not valid and that anomalous flows can be created. The reason behind these anomalous flows is explained by several authors; they argue that constraints can be imposed by the structure of the porous media or by solute-surface interactions, resulting in Lévy flights due to the presence of highly conductive fractures, channels, or macropores or because the infiltration front undergoes jerky movements [26][27][28]. Thus, the application of fractional derivatives was considered to explain the anomalous behavior.
Pachepsky et al. [29] considered water transport in horizontal soil columns and generalized the Richards equation to include a fractional time derivative and solved it by generalizing the Boltzmann transform. Sun et al. [30] considered the same problem; however, noting that the fractional derivative fails with the chain rule, they considered a fractal derivative and showed that their results exhibited an anomalous Boltzmann transform, attributed to the fractal nature of the heterogeneous media.
Generalizing the Fokker-Planck type equation, Su [31] applied the fractional derivative to present a fractional cumulated infiltration into swelling soils where convection dominates, applying the most linear hydraulic conductivity. Subsequently [32], he generalized these results to include mobile and immobile zones with and without gravity effects. Later [33], he introduced mass-time fractional derivatives for swelling-shrinking soils and space-time fractional derivatives for non-swelling soils, applying for mobile zones with or without immobile zones by using power functions for diffusivity and hydraulic conductivity. Physically, the spatial fractional derivative means that the concentration change at the point of observation depends on upstream concentrations, while the temporal fractional derivative implies that the concentration change at the point of observation depends on the prior concentration loading [34].
In the present work, a model for infiltration is developed using a fractional calculation as a tool to incorporate the anomalous behavior in the flow of water. Anomalous horizontal infiltration will first be studied to later extend these results to study anomalous vertical infiltration, where asymptotic solutions will be given for short and long times. Finally, a relationship between conductivity and diffusivity will be assumed to yield a valid anomalous infiltration model for all times.
The work is developed in the following way: Section 2 describes the mathematical tools from fractional calculus that will be implemented; Section 3 develops the equation for the fractional horizontal infiltration considering x = x(θ, t) as the dependent variable and θ, t as independent variables; Section 4 extends the results of the previous section by studying vertical infiltration phenomena for both arbitrary diffusivity and hydraulic conductivity through the asymptotic solution in short times and long times; Section 5 considers, as additional hypotheses, a constant flux-dependant relationship and a relation between diffusivity and hydraulic conductivity to build a model for fractional cumulative infiltration, where its solution is developed for short times expressed as a power series of t ν/2 , where ν is the order of the time fractional derivative with 0 < ν < 2; finally, Section 6 summarizes the conclusions reached in this work.

Fractional Calculus
At present, there is a great diversity of fractional derivatives such as Riemann-Liouville, Caputo, Grünwald, Marchaud, and Riez, among many others. In this section, the basis of fractional calculus that will be applied for the development of this work will be described. For more information regarding these properties, see [35][36][37].
where Γ(α) is the gamma function applied to α defined by The fractional integral is the basis from which various forms of fractional derivative are defined: where n ∈ N with n − 1 < α ≤ n.
Definition 3. The Caputo fractional derivative C D α a+ y of order α ∈ R is defined by where n ∈ N with n − 1 < α ≤ n and y (n) is the n-th derivative.
The Riemann-Liouville fractional derivative and the Caputo fractional derivative are two of the most widely used definitions in fractional calculus; the first because it can be implemented numerically with ease due to the Grünwald-Letnikov algorithm, although it does not allow for the use of the usual initial conditions, while the second allows the classical initial conditions to be implemented when solving fractional differential equations.
In addition, both definitions of the fractional derivative are related by it follows from Equation (18) that both definitions of the fractional derivative are equivalent if y (k) (a) = 0 for k = 0, 1, . . . , n − 1.
In particular, if α, β ∈ R, it can be verified directly from definitions 2 and 3 that Likewise, the Caputo fractional derivative and the Riemann-Liouville fractional integral act among themselves, generalizing the fundamental theorem of calculus, where, for sufficiently good conditions of y(t), it is satisfied that

Fractional Horizontal Infiltration
In this section, the fractional horizontal infiltration equation is developed and solved by replacing the time derivative by the Caputo fractional derivative and neglecting the term associated with the gravitational field.
Note that applying the implicit function theorem to Equation (2), to express θ and t as independent variables and x = x(θ, t) as a dependent variable, it follows that where q = q(θ, t) is the Darcy flow. Therefore, the fractional differential equation that describes the anomalous horizontal infiltration is where ν ∈ (0, 2) is the order of the fractional derivative, and τ c is a constant introduced to maintain dimensional balance with units of time and initial and boundary conditions shown in Equation (3). It can be verified that the solution to anomalous horizontal infiltration, Equation (22), is where φ ν can be called the fractional Boltzmann transform and is defined by the following ordinary differential equation: In addition, fractional cumulative infiltration for fractional horizontal infiltration is defined by where S ν can be called the ν-sorptivity. Note that, when ν → 1, Equations (24) become Equation (4), and Equation (26) becomes Further, Equation (24) and (4) are equivalent if therefore,

Fractional Vertical Infiltration
The fractional differential equation that describes the complete anomalous infiltration is obtained analogously to the previous section, namely, expressing it in closed form, we obtain with initial and boundary conditions as in (7). This equation has been analytically solved for ν = 1 considering particular behaviors for diffusivity and hydraulic conductivity; however, considering arbitrary behavior for D and K, the asymptotic behavior for infiltration can be deduced.
Finally, from Equation (42), the asymptotic behavior at long times for fractional cumulative infiltration is

Fractional Parlange Solution
In the previous sections, the fractional equation that describes the behavior of anomalous infiltration was deduced by neglecting the term associated with the gravitational field, generalizing the Boltzmann transform, as well as the anomalous cumulative infiltration and its relationship with the expressions corresponding to the non-anomalous case; otherwise, for both an arbitrary diffusivity and hydraulic conductivity, the fractional equation was deduced to describe the full anomalous vertical infiltration, and the asymptotic behavior was found for fractional cumulative infiltration in short and long times. Next, considering an arbitrary but related hydraulic diffusivity and conductivity, as well as the approximation of sorptivity given by Parlange [38], the explicit form of fractional equation for anomalous vertical infiltration will be deduced, and its solution will be given in terms of a power series in t ν/2 with ν ∈ (0, 2).
Starting from the flow-concentration relationship applied to fractional flow, Equations (39) and (40) are with 0 ≤ F(θ, t) ≤ 1. Substituting Darcy's law and making a first integration, it follows that Considering the previous expression, the fractional cumulative infiltration is obtained by substituting in Equation (42), namely, Considering this, in long times, the flow concentration relationship behaves as follows: assuming that this behavior is valid for all times. We consider the following relationship between diffusivity and hydraulic conductivity [7]: where β is a parameter associated with the soil type that satisfies 0 < β < 1. Equation (52) interpolates between the Green and Ampt model and the Talsma and Parlange model, for which β = 0 and 1, respectively; moreover, β can be calculated from the experimentally measured hydrodynamic characteristics data using the relation [39] In fact, Parlange et al. [7] recommend β ≈ 0.85 to represent the two soil types sand and clay. Hereafter, β will be considered as a fixed but arbitrary constant parameter.
By substituting Equations (51) and (52) into Equation (50), we have where S is the sorptivity defined by the following approximation [38]: Considering the dimensionless variables the fractional cumulative infiltration, Equation (54), is expressed through the following dimensionless fractional differential equation: with initial condition I Dν (t D = 0) = 0. Note that, for ν = 1, I Dν = I D and Equation (56) can be integrated to obtain the following result: as in [7], where, for β → 0, the last equation is reduced to I D = t D + ln(1 + I D ), as in [1]; for β → 1, the last equation is reduced to I D = t D + 1 − exp(−I D ), as in [6]. The asymptotic behavior of Equation (57) is, for short times, for long times, it is where, for the singular value β = 0, the last expression becomes I D = t D + ln(t D ). Note that, for ν ∈ (0, 1), the integral equation equivalent to Equation (56) is Equation (60) expresses the general solution to fractional cumulative infiltration before different types of soil considering an anomalous flow. However, since the integral term has a singularity for t D → 0, its solution requires numerical methods; nevertheless, an approximate solution in short times can yield useful information.
As in Equation (44), consider the approximation of I Dν (t D ) through the dimensionless power series where S D ν,n for n = 1, 2, . . . are obtained by matching terms. Thus, by substituting the approximation of Equation (61) into the right hand side of Equation (56), we have where, the first four terms are By introducing the approximation shown in Equation (62) into the fractional differential equation, Equation (56), we have where the following coefficients relationship between the series is found: Since b D ν,n depends on S D ν,k for k = 1, 2, . . . , n and n = 1, 2, . . . , it is necessary to solve non-linear equations to find the coefficients of the series I Dν (t D ). The expressions for the first four coefficients of I Dν are shown below: 1 12 It can be seen that the short-time asymptotic behavior shown in Equation (58) is obtained from Equations (67) and (68), making ν = 1. Figure 1 shows the first four terms of the series I Dν (t D ), S D ν,n t nν/2 D for n = 1, 2, 3, 4, where S D ν,n is as in Equations (67)-(70). In particular, it can be observed that the first two terms of the series provide the main contribution to I Dν , which is why it is understood that the approximation to two terms used in the literature is considered sufficient to approximate the cumulative infiltration; however, it is from the fourth term that the contribution provided by S D ν,n is less than 5% of the cumulative infiltration value.
(a) ν = 0.50 (b) ν = 1.00 (c) ν = 1.50 Figure 1. First four terms from Series (61), for the solution in short times to the fractional cumulative infiltration given by the fractional differential Equation (56), considering subdiffusive flow ν = 0.5, classic flow ν = 1.0, and superdiffusive flow ν = 1.5 for a soil type represented by β = 0.85, which is valid for both sands and clays and τ cD = 1. Figure 2 shows the behavior of anomalous infiltration for different values of ν, where ν = 1 is the classic case solved in [7,40]. The phenomena described by the anomalous infiltration with ν < 1 are known as subdifusive phenomena, which indicate that the flow modifies the porous media such that particles blocking fractures or channels are moved, where, for a lower value of ν, there is greater blockage to the flow, showing a lower value to fractional cumulative infiltration. The phenomena described with ν > 1 is known as superdiffusive phenomena, where the flow modifies the porous media such that the particles create new fractures or channels that allow greater mobility to the fluid, where, for a higher value of ν, there is a greater flow and therefore greater fractional cumulative infiltration.  (61) for fractional cumulative infiltration considering subdiffusive flows ν < 1, classical flows ν = 1, and superdiffusive flows ν > 1, for a soil type represented by β = 0.85 that is valid for sands and clays and τ c = 1. Figure 3 shows the effect of β on the fractional accumulated infiltration I Dν (t D ) for several values of ν. As previously described, β is associated with hydrodynamic characteristics of the soil type, where β = 0 corresponds to a linear soil type solved in [1] for ν = 1. β = 1 corresponds to the solution given by Talsma and Parlange solved in [6] for ν = 1. Further, Figure 3 shows how the effect of ν is greater than β in determining the value of the fractional cumulative infiltration I Dν ; that is, the anomalous soil behavior has a greater impact on the cumulative infiltration compared to the soil type. On the other hand, it can be observed that, for a fixed value of ν, the extreme values of β limit the behavior for cumulative infiltration.

Conclusions
A new model is proposed to capture the phenomenon of anomalous infiltration through the use of fractional calculus considering time t and volumetric water content θ as dependent variables and space z as an independent variable. We began the study of fractional horizontal infiltration phenomena by neglecting the term associated with the gravitational field, finding as a solution the fractional Boltzmann transform φ v and the corresponding accumulated infiltration. For the study of fractional vertical infiltration phenomena, the asymptotic solution is found in short and long times for both an arbitrary diffusivity and hydraulic conductivity; the short-time approximation was expressed as a power series in t ν/2 , where ν is the order of fractional derivative, and the first terms are expressed as integro-differential equations. Finally, a constant flow-concentration relationship and a relation between diffusivity and hydraulic conductivity were considered to construct a fractional accumulated infiltration model that is valid for all times and for different soil types. The model was solved through a series that was valid for short times. It is shown that the first two terms of the series provide the main contribution to the anomalous infiltration. It is shown that, for subdifusive flows, there is less accumulated infiltration than in the classic case; for superdifusive flows, the accumulated infiltration is greater than the classic case. The effect of considering different types of soils is also shown, concluding that, for a fixed value of ν, the extreme values of β limit the behavior of cumulative infiltration, and the anomalous soil behavior has a greater impact on cumulative infiltration compared to the soil type.

Funding:
The corresponding author of this work is grateful to CONACYT for the scholarship grant, scholarship number 548429.