Introduction

Water is essential for life and at the heart of civilisations. The Earth is drying and possibly warming forever1, and the world water scarcity is expected to reach the worst level ever by 2025 when 1800 million people are projected to live in countries and regions with “absolute water scarcity” (with < 500 m3 per year per capita), and two-thirds of the world population could be under “stress conditions” (with 500 to 1,000 m3 per year per capita)2. Facing these situations, options are urgently needed to harvest rainwater during increasingly variable rainfall to improve water security for human survival, for peace and the environment in general, and for farmland to maintain food production in particular3. Infiltration of water into soils is one of the most important processes on the natural and urban surfaces, which can serve the purpose of increasing water storage in the subsurface and improving water security.

Infiltration into soils and subsequent recharge of aquifers are essential processes to increase subsurface water storage. Models of infiltration are also a major component in flood forecasting and management, and irrigation scheduling etc. While water shortages and droughts frequently occur in many regions around the world, only 20% of land surface runoff infiltrates the soil and recharges aquifers4 and 80% becomes runoff or floods with the latter damaging to the land surface, infrastructure, and human life. Under the present rainfall and surface conditions, groundwater recharge is decreasing in 40% of world basins5. To change this situation and sustain the environment, increasing infiltration into soils and subsequent recharge of aquifers on both rural natural surfaces and urban areas are essential6,7,8. While the surfaces of rural soils are more natural and their underlying aquifers are more variable, the urban surfaces are, in addition to parks, gardens and other areas in the general category of green infrastructure, can be designed to be more efficient for increased infiltration to reduce floods and improve infiltration for more artificial recharge of groundwater4,9.

In order to assess infiltration and subsequent recharge, a reliable and robust equation of water infiltration into soils is needed. A model of water infiltration into soils is an essential component in flood forecasting, assessing soil and land pollution potentials, scheduling irrigation, and in simulations for various water-related sectors etc. For its importance, major progress was made more than a century ago with the equations named after Green and Ampt10, Horton11 and Philip12 being dominant and treated as classic models. However, significant problems with these classic equations have been singled out by Knight13, which include unrealistic assumptions of instantaneous ponding and sharp front etc. which compromise the physics of water flow by highly simplistic concepts in those models. Forty years on since Knight’s comment13, the problems with “elusive theoretical quantity with arbitrary physical meaning” in the classic equations remain14.

In order to improve the mathematical representation of the physics of flow into soils, Philip15,16 demonstrated that a variety of functions for the moisture-dependent diffusivity in the diffusion equation can be used to replace the constant diffusivity. But Philip’s efforts were only applied to diffusion and limited situations17 without a potential for widespread applications in hydrology and soil science due to the complicated mathematics involved in the solutions.

In the recent two decades, efforts are actively directed to the fractional partial differential equation (fPDE)17,18,19,20,21 to improve our understanding of flow dynamics during water movement in the environment in general and infiltration in particular. The advantage of these new equations of infiltration based on fPDEs is their ability to account for fluctuations in flow velocities during water movement in different pores of the soil in terms of temporal and spatial nonlocalities22, which is lacking in the classic models. The verification of the new equations based on fPDEs against the classic equations supports and encourages the application of the new type of models23.

Built on the recent advances in the fPDEs and the promising performance of the new models23 as well as a versatile numerical method named the homotopy perturbation (HPM)24, this paper presents a set of equations of infiltration. These new equations offer significant findings which have the advantages of (1) eliminating the profound unrealistic assumptions13,25, (2) incorporating random parameters to eliminate the negligence of natural variability of soil properties, and (3) nonlinear power function parameters to account for variabilities in these parameters. The nonlinear fPDE is solved with the powerful HPM and verified using field data, and significant conclusions have been suggested.

Governing equations, equations of infiltration, and their random functional parameters and conditions

The fundamental equation governing the water movement in soils considered here is an fPDE, and it is investigated subject to a special form of the gamma function as the initial condition (IC) and mixed boundary conditions (BCs). The fPDE for water movement in soils with uniform pores takes the form17,18,19,

$$\tau^{\beta - 1} \frac{{\partial^{\beta } \theta }}{{\partial t^{\beta } }} = \frac{\partial }{\partial z}\left( {D(\theta )\frac{\partial \theta }{{\partial z}}} \right) - \frac{\partial K(\theta )}{{\partial z}}$$
(1)

where \(\theta\) is the reduced or relative water content, \(\theta = (\Theta - \Theta_{i} )/(\Theta_{s} - \Theta_{i} )\) with \(\Theta\) the volumetric water content, \(\Theta_{s}\) and \(\Theta_{i}\) being the values of \(\Theta\) at saturation and the initial water content, respectively, \(D(\theta )\) and \(K(\theta )\) are the diffusivity and hydraulic conductivity, respectively, \(z\) is the depth of the soil, \(t\) is time, \(\beta\) is the order of Caputo fractional derivatives with \(0 < \beta \le 1\) for slow- or sub-diffusion type of flow in soil pores, and the term \(\tau^{\beta - 1}\) functions as a dimensional correction to ensure the final result conforms to the conventional dimensions.

An extension to the fPDE in Eq. (1) is the distributed-order fPDE in which \(b_{1} \frac{{\partial^{{\beta_{1} }} \theta }}{{\partial t^{{\beta_{1} }} }} + b_{2} \frac{{\partial^{{\beta_{2} }} \theta }}{{\partial t^{{\beta_{2} }} }}\) replaces \(\frac{{\partial^{\beta } \theta }}{{\partial t^{\beta } }}\) in Eq. (1) to account for water movement in large and small pores or in mobile and immobile zones of the soil. It should be noted that the time-fractional order, \(\beta\), (or distributed orders) in a PDE is the parameter of waiting time distributions in the continuous-time random walk (CTRW) model which explains the movement of water parcels18. When this parameter appears in the equation of infiltration, it models the movement of water in soil pores in terms of the CTRW theory.

Soils to be analysed here can be non-swelling and swelling soils with Eq. (1) for uniform non-swelling soils where the physical coordinate in the vertical direction is denoted by \(z\), while for swelling soils the material coordinate, denoted by \(m\), is used which includes a factor \(\left( {\gamma_{n} \alpha - 1} \right)\) in the diffusivity and hydraulic conductivity to account for swelling properties of the soil.

The diffusivity and hydraulic conductivity in Eq. (1) for water flow in non-swelling soils take the forms of \(D(\theta ) = D_{0} \theta^{c}\) and \(K(\theta ) = K_{0} \theta^{n}\), respectively, where \(D_{0}\), \(c\), \(K_{0}\) and \(n\) are soil and flow parameters. For swelling soils, both the diffusivity and hydraulic conductivity take similar forms except for a factor \(\left( {\gamma_{n} \alpha - 1} \right)\) included in the fPDE to account for swelling properties (see details in the Methods).

Solutions of the above fPDEs are derived using the versatile iterative numerical method HPM. This iterative numerical method can efficiently solve an fPDE with either an IC, BC, or any combinations of both an IC and BCs. In this article, the fPDEs for non-swelling and swelling soils in uniform pores, and the distributed-order fPDEs for flow in soils with large and small pores are solved for a special form of the gamma function, and mixed BCs.

The special gamma function as the IC takes the form of

$$\theta (z,t) = \theta_{0} \left( {1 + az} \right)e^{ - bz}$$
(2)

where \(\theta_{0}\) is the moisture ratio on the soil surface, and the two parameters, \(a\) and \(b\), determine the shapes of the moisture profile. Equation (2) for \(a = 0\) becomes the exponential function, \(\theta (z,t) = \theta_{0} e^{ - bz}\). The distributions of the soil moisture profiles described by the gamma function are depicted in Fig. 1, which appear under different situations.

Fig. 1
figure 1

The simulated initial moisture distributions using Eq. (2). (A) \(\theta_{0} = 0.40\), \(a = 0.6\) and different values for \(b\); (B) \(\theta_{0} = 0.40\), \(b = 0.35\) and different values for \(a\) with \(a = 0\) representing the exponential distribution.

Full size image

The mixed BCs take the following forms

$$\theta (0,t) = f_{1} (t),\frac{\partial \theta (0,t)}{{\partial z}} = f_{2} (t)$$
(3)

where \(f_{1} (t)\) and \(f_{2} (t)\) are random variables or functions to be analysed in detail with \(f_{1} (t)\) as the moisture content or ratio on the surface, and \(f_{2} (t)\) as the gradient of the moisture across the soil profile. Both \(f_{1} (t)\) and/or \(f_{2} (t)\) can be random or deterministic functions, variables or constants.

The equations of cumulative infiltration and infiltration rate

The equations of cumulative infiltration, \(I(t)\), derived using HPM subject to the IC and mixed BCs detailed in the Methods section all take the same form given by

$$I(t) = B + At + Ft^{\beta }$$
(4)

With different expressions for \(B\) and \(F\) under different ICs and BCs, where \(A\) is the final infiltration rate, \(B\) is a measure of the initial depth of water in the surface soil and in the soil profile, and \(F\) is an integrated transport parameter. The infiltration rate is given by differentiating Eq. (4) with respect to time,

$$i(t) = A + \beta Ft^{\beta - 1}$$
(5)

The different forms of the equations of cumulative infiltration and infiltration rates are summarised in Table 1.

Table 1 Generic equation of cumulative infiltration \(I(t) = B + At + Ft^{\beta }\) with different forms of \(B\) and \(F\) derived with different ICs and mixed BCs.
Full size table

While the equations of infiltration derived from different ICs and BCs can be used, the equations derived with the gamma function IC are the most generic form. The rigorously derived Eq. (4) also supports an empirical model which is similar to Eq. (4) in structure and is regarded as “A general formula exhibiting observed features…” of field infiltration known as the (three-term) modified Kostiakov equation (MKE)26 or Kostiakov-Lewis equation (KLE)27 which, in fact, was analysed earlier by Mezencev28.

The initial condition in Eq. (2) used for deriving Eq. (4) has physical significance: (1) the case with \(a < 0\) corresponds to a decreasing soil moisture profile from the surface initially, and (2) the case with \(a > 0\) in the initial condition means that the soil moisture is low on the surface initially and increases downwards for certain depths, and then decreases again as depicted in Fig. 1B.

The boundary conditions used here are also physically important which result in the same form of the equation of cumulative infiltration. As such the users can apply the general equation with different choice of explanations for the soil and flow conditions.

For infiltration into swelling soils, as it is summarised in Table 1 that the equation of cumulative infiltration is similar to Eq. (4) except that the transport parameter \(F\) is replaced by \(F_{s}\) with a factor \((\gamma_{n} \alpha - 1)\) added to account for swelling, i.e., \(F_{s} = F(\gamma_{n} \alpha - 1)\).

It is important to point out the differences between Eq. (4) and other reported empirical equations of cumulative infiltration: firstly Eq. (4) is rigorously derived from a nonlinear random fPDE with the diffusivity and hydraulic conductivity being random power functions; secondly an arbitrary moisture profile is used which is more realistic under natural conditions without unrealistic assumptions; and thirdly the random parameters accommodate natural variabilities of the parameters.

Fitting the data of cumulative infiltration measured at different initial moisture contents29 to Eq. (4) yields the parameters that are listed in Table 2. It should be noted that surface conditions are different for the same soil which result in different values for the parameters in Eq. (4). Given the very small values of \(B\) in the soils analysed in this paper as a special case, Eq. (4) can be approximated by

$$I(t) = At + Ft^{\beta }$$
(6)

In Table 2, the difference in Eqs. (4) and (6) are compared with published data measured in the field.

Effects of initial soil moisture contents on infiltration

The same data of Talsma29 was fitted to Eq. (6) which yields the parameter values that are also listed in Table 2. Two more data sets30,31 without information on the initial moisture contents are also analysed, respectively, with Eq. (4) and Eq. (6) and the results are listed in Table 2 for comparisons. 

Table 2 Comparisons of parameters in Eqs. (4) and (6).
Full size table

It can be seen in Table 2 that the larger the initial moisture content (or the background flow), the larger is the value of \(A\), which supports the concept of background flow. The three curves and corresponding different levels of initial moisture data are shown in Fig. 2, which, together with the data in Table 2, indicate that the parameter \(B\) is very small for this soil and can be neglected for simplicity in this case. These equations are derived from the random concept and can be also used as deterministic equations as an option.

Fig. 2
figure 2

The effect of the initial moisture content on infiltration. The data from Talsma29 were fitted to Eq. (6).

Full size image

Effects of soil types on infiltration

In addition to the three-term MKE26 identical to Eq. (4), there is another empirical equation of infiltration called (two-term) MKE identical to Eq. (6) in practice32,33, and there are extensive reports on the variability of the parameters in the two-term MKE. Here, two important issues are investigated using Eq. (6) with data originally reported for the two-term MKE: the effects of soil types and land surface management on infiltration.

The types of soil can have a significant impact on infiltration physics which determines the rate and quantity of infiltration. With the reported data on the infiltration rates measured in the field for different types of soils32, the parameters in the two-term MKE are derived by curve-fitting and converted to the equivalent parameters in Eq. (6) here which are listed in Table 3 for comparison.

Table 3 Variability in parameters of Eq. (6) for infiltration into different types of soils with data converted from a published report33.
Full size table

Effects of land surface management on infiltration

Land surface management can also have a significant impact on infiltration physics. Tillage and land surface covers are two important forms of land management practice, which are examined here using the data and parameters derived for the two-term MKE33 and are translated into the parameters in Eq. (6) for re-interpretation and comparison in Table 4. The data listed in Table 4 indicate that surface covers and tillage have different consequences for infiltration which result in different levels of rainfall recharge in the subsurface.

Table 4 Response of infiltration physics to land surface treatments measured by the parameters in the equation of cumulative infiltration.
Full size table

The results in Table 4 indicate that the land management practice with tillage and mulching is the most effective practice for a higher infiltration rate followed by the mulched no-tillage land surface. These findings are very important for land management guidelines in a drying and warming climate to promote rainfall infiltration into soils for better water security, and particularly important for the farming sector around the world for food security. An optimal land surface management is also significant for its low-cost operation for the agricultural sector, for water conservation and climate adaptation without causing too much disturbance on the land and the environment.

Effects of irrigation water quality and management practice on infiltration

In irrigated agriculture in different regions of the world, many factors affect infiltration into soils, which include irrigation water quality that is characterised by variable chemical concentrations27, and different land management options that alter the surface conditions which include tillage, mulching, burning and residue management etc.27,34,35,36. The highly variable values of the three parameters in Eq. (6) are visible in terms of the KLE27.

Conclusion

A very general equation of cumulative infiltration has been presented for assessing water infiltration into soils to recharge the subsurface water storage. This equation of cumulative infiltration in Eq. (4) is shown to generalise a set of equations of cumulative infiltration derived using different initial conditions (ICs) and mixed boundary conditions (BCs). This generalisation implies that the same physical process of infiltration into soils can be represented by one equation and interpreted with different mathematical terminologies whether they are ICs or BCs. These equations have been derived from the very active theory of random fractional partial differential equations without unrealistic assumptions which compromise the physics of flow into soils. The procedures presented here demonstrate that the governing equation with the general nonlinear diffusivity and hydraulic conductivity can be solved with the aid of the homotopy perturbation method with either an initial condition or mixed boundary conditions without compromising the physics and mathematics which describe the movement of water in soils. A summary of this set of equations is given so that readers will have a choice for what is suitable for their situations. An important demonstration of this finding is the evaluation of the effects of land surface covers and tillage to improve infiltration into the subsurface33, which is a very important practice in water management. When the practices of tillage, no tillage and occasional tillage37 are optimised, the reductions in surface runoff and soil erosion on land will accompany the increase in water infiltration and improvement in other environmental benefits. The quantity of water infiltrated into the same soil can be increased by altering land management practices, and the implications are particularly important for managing large areas of farmland in the drying climate without high costs associated with major engineering works which cause major environmental concerns.

Methods

This section is the summary of Methods for the solutions of two flow equations using the HPM. Two governing equations are used here: one fPDE is for uniform pores in soils, and the second fPDE for flow in soils with large and small pores or mobile-immobile zones.

The equation of flow in soils with uniform pores

One of the fPDEs presented by the author17,18,19 for water movement in rigid soils analysed in this paper is of the form,

$$\tau^{\beta - 1} \frac{{\partial^{\beta } \theta }}{{\partial t^{\beta } }} = \frac{\partial }{\partial z}\left( {D(\theta )\frac{\partial \theta }{{\partial z}}} \right) - \frac{\partial K(\theta )}{{\partial z}}$$
(7)

where \(D(\theta )\) is the diffusivity, \(K(\theta )\) is the hydraulic conductivity, and both of which are related by the relationship37

$$D(\theta ) = K(\theta )\frac{d\psi }{{d\theta }}$$
(8)

with \(\psi\) being the water potential. For water movement in swelling soils, Eq. (7) is modified as19,20

$$\tau^{\beta - 1} \frac{{\partial^{\beta } \theta }}{{\partial t^{\beta } }} = \frac{\partial }{\partial m}\left[ {D_{m} (\theta )\frac{\partial \theta }{{\partial m}}} \right] - \left( {\gamma_{n} \alpha - 1} \right)\frac{{dK_{m} (\theta )}}{d\theta }\frac{\partial \theta }{{\partial m}}$$
(9)

where \(m\) is the material coordinate, \(\gamma_{n}\) is the particle specific gravity, \(\alpha\) 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, \(\theta\), and \(D_{m}\left(\theta\right)\) is the material diffusivity given by Smiles and Raats38

$$D_{m} (\theta ) = \frac{{K_{m} (\theta )}}{1 + \theta }\frac{d\Phi }{{d\theta }}$$
(10)

with \(\Phi\) being the unloaded matrix potential, and \(K_{m} \left( \theta \right)\) being the material-based unsaturated material hydraulic conductivity. With the saturated material hydraulic conductivity given as \(k_{m} = K\left( \theta \right)\theta_{s}\), where \(\theta_{s}\) is the fraction of solids and \(K\left( \theta \right)\) the conventional unsaturated hydraulic conductivity, the material-based hydraulic conductivity takes the form of

$$K_{m} (\theta ) = K(\theta )(\gamma_{n} \alpha - 1)$$
(11)

The equation of flow in soils with large and small pores

When the pores in soils are not uniform, but can be grouped into large and small pores, a distributed-order fPDE18,39 can be used,

$$\tau^{\beta - 1} \left( {b_{1} \frac{{\partial^{{\beta_{1} }} \theta }}{{\partial t^{{\beta_{1} }} }} + b_{2} \frac{{\partial^{{\beta_{2} }} \theta }}{{\partial t^{{\beta_{2} }} }}} \right) = \frac{\partial }{\partial z}\left( {D(\theta )\frac{\partial \theta }{{\partial z}}} \right) - \frac{dK(\theta )}{{dz}}\frac{\partial \theta }{{\partial z}}$$
(12)

for non-swelling soils with \(0 < \beta_{1} < \beta_{2} \le 1.0\). In this paper, both fPDEs for uniform and non-uniform pores are analysed. For swelling soils, Eq. (11) shows that the difference between the material-based hydraulic conductivity, \(K_{m} \left( \theta \right)\), and its conventional counterpart, \(K(\theta )\) is only a factor \((\gamma_{n} \alpha - 1)\).

Essentials for the homotopy perturbation method for solutions of the governing equations

The HPM is an iterative numerical method developed by He24 which has extensive applications in applied mathematics and numerical analysis for different types of ordinary differential equations (ODEs), partial differential equations (PDEs) and their fractional counterparts such as the fPDEs. In water-related fields, a compendium of HPM applications has been published40, and the essentials of this method are outlined here only for the convenience of readers.

Homotopy is a concept in topology which is central to the HPM and related methods. Topology is a field in mathematics which deals with geometric properties of a mathematical object which remains unaffected under continuous deformation40.

Homotopy defines a connection between two mathematical objects by specifying that the two objects are homotopic if one can be continuously deformed into the other. This concept is the key to the approximate solution of a differential equation when the original problem is expanded in a mathematical expression known as the homotopy.

Take two functions, \(f(t)\) and \(g(t)\), of a dimension \(t\) (either for time or space) for example, a homotopy is defined as

$$H(t;q) = (1 - q)f(t) + qg(t),$$
(13)

which is defined in the range of zero and one, namely, \(q \in [0,1][0,1]\), where \(q\) is the perturbation parameter. The homotopy \(H(t;q)\) between \(f(t)\) and \(g(t)\) is itself a continuous function, and Eq. (13) shows that \(H(t;q) = f(t)\) for \(q = 0\) and \(H(t;q) = g(t)\) for \(q = 1\).

For ODEs, PDEs, fractional ODEs (fODEs) and fPDEs, the HPM has been extensively applied to develop solutions in a series form24. For an analytical function, \(f(r)\), with a general differential operator, \(D\), define

$$D(u) = f(r)$$
(14)

and the operator \(D\) can be split into two parts, a nonlinear component denoted by \(N\), and a linear component by \(L\). Then a nonlinear differential equation can be written as

$$N(u) + L(u) - f(r) = 0$$
(15)

A homotopic of the general differential equation in Eq. (14) is written as24,41,

$$H(v,p) = (1 - p)\left[ {L(v) - L(u{}_{0})} \right] + p\left[ {D(v) - f(r)} \right] = 0,p \in [0,1]$$
(16)

where \(u{}_{0}\) is the initial approximation of Eq. (14). Note that a series decomposition solution is defined in term of \(v\) as approximation to \(u\). Equation (16) defines

$$H(v,0) = L(v) - L(u{}_{0}) = 0$$
(17)

and

$$H(v,1) = D(v) - f(r) = 0$$
(18)

with \(p\) being the perturbation parameter, the solution of Eq. (16) is assumed to be a series solution of the form,

$$\begin{aligned} \nu \left( {x,p} \right) & = \nu_{0} \left( x \right) + \nu_{1} \left( x \right)p + \nu_{2} \left( x \right)p^{2} + \ldots + \\ & = \sum\limits_{k = 0}^{\infty } {\nu_{k} \left( x \right)p^{k} } \\ \end{aligned}$$
(19)

which, for \(p = 1\), gives the desired solution of Eq. (16) as

$$v(x,p) = v_{0} (x) + v_{1} (x) + v_{2} (x) + ... + = \sum\limits_{k = 0}^{\infty } {v_{k} (x)}$$
(20)

As an alternative to traditional methods43, HPM has been applied to solve various types of nonlinear problems with satisfactory accuracy in nonlinear and complex fluid flow problems compared with the traditional approximation methods and classic numerical methods44.

Equations for infiltration derived using the homotopy perturbation method

Non-swelling soils with uniform pores

Based on the earlier development18,44, we further investigate the random boundary conditions and parameters in fPDEs for water flow in soils. The following nonlinear diffusivity, \(D(\theta )\), and the hydraulic conductivity, \(K(\theta )\), are used,

$$D(\theta ) = D_{0} \theta^{c} ,K(\theta ) = K_{0} \theta^{n}$$
(21)

The special form of the gamma function as the initial condition is given here as

$$\theta (z,t) = \theta_{0} \left( {1 + az} \right)e^{ - bz} \,{\text{at}}\,t = 0\,and\,z > 0$$
(22)

where \(\theta_{0}\) is the water content or water ratio on the surface at \(t = 0\), \(a\) and \(b\) are constants which determine the shape of the initial moisture distribution across the soil profile. In particular, the values of \(a\) ranges from negative to positive values, and a positive \(a\) implies the increase in the moisture from the surface of the soil while a negative \(a\) means the decrease of the soil moisture from the surface.

Simulations indicate that the shape parameter \(a\) is in the range of \(- 2 < a < 2\) which is suitable for many natural shapes of the initial moisture profile. In the special case of \(a = 0\), Eq. (22) becomes the exponential function, \(\theta (z,t) = \theta_{0} e^{ - bz}\) at \(t = 0\) and \(z > 0\). Both forms of the initial conditions indicate that the initial moisture content at the surface is \(\theta (0,0) = \theta_{0}\) at \(z = 0\) and \(t = 0\).

The solutions of Eq. (7) with the two coefficients as functions in Eq. (21) subject to the condition in Eq. (22) are, respectively, given as an infinite series of the approximate solution, \(u\),

$$\theta (z,t) = \sum\limits_{i = 0}^{\infty } {u_{i} }$$
(23)

The procedures for the solution with the initial condition Eq. (22) derived using the HPM can be found in the literature41,45,46,47,48,49,50,51. The general form of the HPM procedure for solutions of Eq. (7) subject to Eq. (22) with the diffusivity and hydraulic conductivity in Eq. (21) is written as

$$\left\{ \begin{gathered} \tau^{\beta - 1} \frac{{\partial^{\beta } \theta }}{{\partial t^{\beta } }} - p\left[ {\frac{\partial }{\partial z}\left( {D_{0} \theta^{c} \frac{{\partial^{2} \theta }}{{\partial z^{2} }}} \right) - nK_{0} \theta^{n - 1} \frac{\partial \theta }{{\partial z}}} \right] = 0 \hfill \\ \theta (z,t) = \theta_{0} \left( {1 + az} \right)e^{ - bz} ,t = 0 \hfill \\ \end{gathered} \right.$$
(24)

where the perturbation parameter varies as \(0 \le p \le 1\). For \(p = 0\), Eq. (24) becomes41,45,46

$$\left\{ \begin{gathered} p^{0} :\tau^{\beta - 1} \frac{{\partial^{\beta } u_{0} }}{{\partial t^{\beta } }} = 0 \hfill \\ u_{0} = \theta_{0} \left( {1 + az} \right)e^{ - bz} ,t = 0 \hfill \\ \end{gathered} \right.$$
(25)

which yields the first component of the solution

$$u_{0} = \theta_{0} \left( {1 + az} \right)e^{ - bz}$$
(26)

For \(p = 1\), Eq. (24) with a homogeneous condition is written as

$$\left\{ \begin{gathered} p^{1} :\tau^{\beta - 1} \frac{{\partial^{\beta } u_{1} }}{{\partial t^{\beta } }} - \left[ {\frac{\partial }{\partial z}\left( {D_{0} u_{0}^{c} \frac{{\partial u_{0} }}{\partial z}} \right) - nK_{0} u_{0}^{n - 1} \frac{{\partial u_{0} }}{\partial z}} \right] = 0 \hfill \\ u_{1} = 0,t = 0 \hfill \\ \end{gathered} \right.$$
(27)

where \(u_{0}\) is given by Eq. (26). Equation (27) is a nonlinear fPDE, and its solution is given as follows by fractional integration52

$$u_{1} = \left\{ \begin{gathered} D_{0} \left[ {\theta_{0} \left( {a - b - abz} \right)} \right]^{2 + c} e^{ - (2 + c)bz} \hfill \\ + cD_{0} \left[ {\theta_{0} \left( {a - b - abz} \right)} \right]^{1 + c} e^{ - (1 + c)bz} \hfill \\ - nK_{0} \left[ {\theta_{0} \left( {a - b - abz} \right)} \right]^{n} e^{ - nbz} \hfill \\ \end{gathered} \right\}\frac{{\tau^{1 - \beta } }}{\Gamma (1 + \beta )}t^{\beta }$$
(28)

and the two-order approximate solution is the sum of the order zero \(u_{0}\) in Eq. (26), and \(u_{1}\) in (28), i.e., \(\theta = \theta_{0} \left( {1 + az} \right)e^{ - bz} + u_{1}\), which is used in the following integral37

$$I(t) = At + \int\limits_{0}^{\infty } {\theta dz}$$
(29)

to derive the equation of cumulative infiltration. Equation (29) can be completed to yield a very complicated form, and to simplify the presentation, it is reasoned that the term \(abz\) can be neglected as the values of \(a\) and \(b\) are very small (see Fig. 1), particular for small depths associated with infiltration, then, the integration in Eq. (29) yields

$$I(t) = B + At + Ft^{\beta }$$
(30)

with

$$B = \frac{1}{b}\left( {\theta_{0} + \frac{a}{b}} \right)$$
(31)

and

$$F = \left\{ \begin{gathered} \frac{{D_{0} \left[ {\theta_{0} \left( {a - b} \right)} \right]^{2 + c} }}{(2 + c)} + \frac{{cD_{0} }}{(1 + c)}\left[ {\theta_{0} \left( {a - b} \right)} \right]^{1 + c} \hfill \\ - K_{0} \left[ {\theta_{0} \left( {a - b} \right)} \right]^{n} \hfill \\ \end{gathered} \right\}\frac{{\tau^{1 - \beta } }}{b\Gamma (1 + \beta )}$$
(32)

where A is the final infiltration rate, \(B\) accounts for the initial quantity of the moisture contents and is a constant for a saturated soil, and \(F\) is an integrated parameter. As \(\theta_{0}\) is continuous and random by definition, \(B\) is also continuous and random, and its derivative with respect to time is \(dB(\theta_{0} )/dt = 0\). The infiltration rate is given by differentiating Eq. (30) with respect to time,

$$i(t) = A + \beta Ft^{\beta - 1}$$
(33)

When the exponential function as the initial condition is used, the solution of Eq. (24) is slightly different, and the equation of cumulative infiltration is identical in form to Eq. (30) except for the following parameters,

$$B = \frac{{\theta_{0} }}{b}$$
(34)

and

$$F = \left[ {bD_{0} \theta_{0}^{c + 1} \left( {\frac{{\theta_{0} }}{2 + c} + \frac{c}{1 + c}} \right) + K_{0} \theta_{0}^{n} } \right]\frac{{\tau^{1 - \beta } }}{\Gamma (\beta + 1)}$$
(35)

which is not fully identical to Eq. (32) for \(a = 0\) because of the approximation in Eq. (35) by neglecting the term \(abz\). The infiltration rate in this case is similar to Eq. (33) but with \(B\) in Eq. (34) and \(F\) in Eq. (35) instead.

Swelling soils with uniform pores

For swelling soils, Eq. (9) is the governing equation for water movement. Comparing the formulation for non-swelling soils and its solutions as well as the equations of infiltration, the solution of Eq. (9) derived using the HPM for swelling soils is given as

$$\theta = \theta_{0} \left( {1 + am} \right)e^{ - bm} + u_{1}$$
(36)

The equation of cumulative infiltration based on Eqs. (29) and (36) is

$$I(t) = B + At + F_{s} t^{\beta }$$
(37)

with

$$B = \frac{1}{b}\left( {\theta_{0} + \frac{a}{b}} \right)$$
(38)

and

$$F_{s} = \left\{ \begin{gathered} \frac{{D_{0} \left[ {\theta_{0} \left( {a - b} \right)} \right]^{2 + c} }}{(2 + c)} + \frac{{cD_{0} }}{(1 + c)}\left[ {\theta_{0} \left( {a - b} \right)} \right]^{1 + c} \hfill \\ - K_{0} \left[ {\theta_{0} \left( {a - b} \right)} \right]^{n} \hfill \\ \end{gathered} \right\}\frac{{\tau^{1 - \beta } (\gamma_{n} \alpha - 1)}}{b\Gamma (1 + \beta )}$$
(39)

The rate of infiltration into swelling soils is given by \(i(t) = A + \beta Ft^{\beta - 1}\) with \(F\) given by Eq. (39).

Non-swelling soils with large and small pores

In this case, Eq. (12) is used, and the homotopy above is extended to this case,

$$\left\{ \begin{gathered} p^{1} :\tau^{\beta - 1} \left( {b_{1} \frac{{\partial^{{\beta_{1} }} u_{1} }}{{\partial t^{{\beta_{1} }} }} + b_{2} \frac{{\partial^{{\beta_{2} }} u_{1} }}{{\partial t^{{\beta_{2} }} }}} \right) - \left[ {\frac{\partial }{\partial z}\left( {D_{0} u_{0}^{c} \frac{{\partial u_{0} }}{\partial z}} \right) - nK_{0} u_{0}^{n - 1} \frac{{\partial u_{0} }}{\partial z}} \right] = 0 \hfill \\ u_{1} = 0,u_{0} = \theta_{0} \left( {1 + az} \right)e^{ - bz} ,t = 0 \hfill \\ \end{gathered} \right.$$
(40)

with \(u_{0} = \theta_{0} \left( {1 + az} \right)e^{ - bz}\) as the first component of the solution. The second term on the right-hand-side of Eq. (40) is the water supply rate on the surface, and can be written as

$$q = - \frac{\partial }{\partial z}\left[ {D_{0} u_{0}^{c} \frac{{\partial u_{0} }}{\partial z} - K_{0} u_{0}^{n} } \right]$$
(41)

which, in the content of infiltration, is equivalent to the infiltration rate, i.e., \(q = i(t)\). By applying Laplace transform to Eq. (40), rearranging the term \(\left( {b_{1} s^{{\beta_{1} }} + b_{2} s^{{\beta_{2} }} } \right)\), and then applying the inverse Laplace transform gives

$$u_{1} = \tau^{1 - \beta } e_{\alpha ,\beta } \left( {t;\lambda } \right)i(t)$$
(42)

where \(\lambda = \frac{{b_{2} }}{{b_{1} }}\), and \(e_{\alpha ,\beta } \left( {t;\lambda } \right)\) is the generalised Prabhakar function (GPF)53,54 given by

$$e_{\alpha ,\beta } \left( {t;\lambda } \right) = t^{{\beta_{2} - 1}} E_{{\beta_{2} - \beta_{1} ,\beta_{2} }} \left( { - \lambda t^{\alpha } } \right)$$
(43)

with \(E_{{\beta_{2} - \beta_{1} ,\beta_{2} }} \left( { - \lambda t^{\alpha } } \right)\) being the Mittag–Leffler function53, \(\alpha = \beta_{2} - \beta_{1}\) and \(\beta = \beta_{2}\). The two-term approximation of the solution for the water content or water ratio in the soil is then given by

$$\theta (z,t) = u = \theta_{0} \left( {1 + az} \right)e^{ - bz} + \tau^{1 - \beta } e_{\alpha ,\beta } \left( {t;\lambda } \right)i(t)$$
(44)

which, for infiltration that applies to the surface at \(z = 0\), is rearranged to give

$$i(t) = \frac{{\tau^{\beta - 1} \left[ {\theta (0,t) - \theta_{0} } \right]}}{{e_{\alpha ,\beta } \left( {t;\lambda } \right)}}$$
(45)

where \(\theta (0,t)\) is the water content or water ratio on the surface as a function of time compared to the initial water content on the surface, \(\theta_{0}\). Equation (45) relates the infiltration rate to the gradient of the moisture content or ratio on the surface. The GPF in Eq. (43) and Eq. (44) has asymptotic solutions under different conditions as time approaches large or \(t \to + \infty\). For \(0 < \beta_{1} < \beta_{2} \le 1.0\)54, the asymptotic result is

$$e_{\alpha ,\beta } \left( {t;\frac{{b_{2} }}{{b_{1} }}} \right) \sim \frac{1}{{\Gamma (\beta_{1} )}}\left( {\frac{{b_{2} }}{{b_{1} }}t} \right)^{{\beta_{1} - 1}} \,0 < \beta_{1} \le 1.0$$
(46)

which shows that as the time approaches very large, Eq. (45) can be written as

$$i(t) \sim \tau^{\beta - 1} \left[ {\theta (0,t) - \theta_{0} } \right]\Gamma (\beta_{1} )\left( {\frac{{b_{2} }}{{b_{1} }}t} \right)^{{1 - \beta_{1} }}$$
(47)

which means that as time becomes very large, infiltration into the soil is determined by the ratio of porosities, \(\lambda\), and small pores represented by the order of fractional derivatives for small pores, \(\beta_{1}\). Here the symbol \(\sim\) means “approaches” as time becomes very large or asymptotically, and for the purpose of approximate expressions, it can be replaced by the approximately equal sign \(\approx\) in practice provided that the large time condition is stated.

If analysis and applications are desired to relate infiltration to more soil parameters, two options can be taken based on the infiltration integral55 with Eq. (45) incorporated: one is the Parlange-Smith equation (PSE)55,56, and the second is the Green and Ampt equation (GAE)10.

With the PSE, cumulative infiltration, \(I(t)\), is given by combining Eq. (45) with the PSE55 to yield,

$$I(t) = At + \left( {\theta_{s} - \theta_{0} } \right)G\ln \left( {\frac{{i(t) - K_{i} }}{{i(t) - K_{s} }}} \right)$$
(48)

where \(i(t)\) is given by Eq. (45), and

$$G = \frac{{D_{0} \left( {\theta_{s}^{c + 1} - \theta_{0}^{c + 1} } \right)}}{{(c + 1)\left( {K_{s} - K_{i} } \right)}}$$
(49)

with \(\theta_{s}\) and \(K_{s}\) being the saturated water content (or water ratio) and saturated hydraulic conductivity of the soil, respectively. At ponding time when the surface water supply rate such as the rainfall intensity equals the infiltration rate, i.e., \(i(t) = i_{p}\), the cumulative infiltration in Eq. (48) is expressed as \(I_{p} (t)\).

With the GAE, the equation of cumulative infiltration55 is combined with Eq. (45) and with \(A = K_{s}\) to yield

$$I(t) = At + \frac{S}{i(t) - A}$$
(50)

with

$$S = \frac{{D_{0} \left( {\theta_{s} - \theta_{0} } \right)\left( {\theta_{s}^{c + 1} - \theta_{0}^{c + 1} } \right)}}{(c + 1)}$$
(51)

The assignment of \(A = K_{s}\) make the discussion in Smith55 more general because \(A\) covers more cases including \(K_{s}\), which has been discussed by Youngs57 and Philip58. Equation (50) can be rearranged to give the infiltration rate, \(i(t) = A + \frac{S}{I(t) - At}\), which means that the infiltration rate \(i(t)\) approaches the final infiltration rate, \(A\), as time becomes very large, i.e., \(i(t) \to A\) as \(t \to \infty\). Note that \(i(t)\) in Eqs. (48) and (50) can be replaced by Eq. (45).

Solutions of the fPDE and equation of cumulative infiltration subject to mixed boundary conditions

The solutions of Eq. (7) can be derived subject to the following mixed boundary conditions49,59,

$$\theta (0,t) = f_{1} (t),\frac{\partial \theta (0,t)}{{\partial z}} = f_{2} (t)$$
(52)

with both \(f_{1} (t)\) and \(f_{2} (t)\) as random variables and/or deterministic functions or constants.

In the case of constant \(f_{1} (t) = \theta_{0}\) and \(f_{2} (t) = q_{0}\) in Eq. (52), the first component of the governing equation using the HPM is given as

$$u_{0} = \theta_{0} + zq_{0}$$
(53)

and the next order component of the solution is given by solving the fractional differential equation

$$u_{1} = \left[ {cq_{0} D_{0} \left( {\theta_{0} + zq_{0} } \right)^{c - 1} - nK_{0} \left( {\theta_{0} + zq_{0} } \right)^{n - 1} } \right]\frac{{q_{0} \tau^{1 - \beta } t^{\beta } }}{\Gamma (\beta + 1)}$$
(54)

The equation of cumulative infiltration can be shown to be

$$I(t) = B_{b} + At\, + F_{b} t^{\beta }$$
(55)

where

$$B_{c} = \frac{{\theta_{0}^{2} }}{{q_{0} }}$$
(56)

is the initial depth of water content in the soil, and

$$F_{b} = \frac{{\tau^{1 - \beta } }}{\Gamma (\beta + 1)}\left( {K_{0} \theta_{0}^{n} - D_{0} q_{0} \theta_{0}^{c} } \right)$$
(57)

Swelling soils with large and small pores

As with the uniform soils, the key difference is the factor \(\left( {\gamma_{n} \alpha - 1} \right)\) associated with the hydraulic conductivity and conductivity for swelling soils. As discussed in the main text, the distributed-order fPDE \(b_{1} \frac{{\partial^{{\beta_{1} }} \theta }}{{\partial t^{{\beta_{1} }} }} + b_{2} \frac{{\partial^{{\beta_{2} }} \theta }}{{\partial t^{{\beta_{2} }} }}\) replaces \(\frac{{\partial^{\beta } \theta }}{{\partial t^{\beta } }}\) in Eq. (9) to account for water movement in large and small pores or in mobile and immobile zones of the soil. Then the formulation above from Eq. (52) to Eq. (57) for non-swelling soils can be extended to swelling soils by including \(\left( {\gamma_{n} \alpha - 1} \right)\) to where \(D_{0}\) and \(K_{0}\) appear so that \(F_{m} = \left( {\gamma_{n} \alpha - 1} \right)F_{b}\) for swelling soils.

The procedures with the homotopy perturbation method leading to the above solutions are based on the mean square approach to random variables which are similar to the mean square Laplace transform60 in that the random variables \(A\), \(B\) and \(F\) discussed here are, in fact, synonyms of the expressions \(A = E[A_{n} ]\) for simplicity, where \(E\) symbolises expectation (or mean), and each random variable is an expectation of many observations, i.e., \(E[A_{n} ]\to _{n \to \infty } E[A] \to A\), \(E[F_{n} ]\to _{n \to \infty } E[F] \to F\), and \(E[B_{n} ]\to _{n \to \infty } E[B] \to B\). In hydrological practice such as during parameter estimation etc., the parameters \(A\), \(B\) and \(F\) can also be treated as deterministic quantities so that no sophistication in the concept is needed.