METHOD
10.1029/2023WR035782
Key Points:• A tempered stable law withsubordination built a unified model
governing water equilibrium andkinetics in saturated/unsaturated soils• Fractional‐derivative equations fitted
soil‐water data sets, such as non‐Darcian flow, as well as or better thanclassical models• The assumed tempered stable density
distribution for soil‐hydraulic properties led to multi‐rate moisturedynamics

Supporting Information:
Supporting Information may be found in
the online version of this article.

Correspondence to:Y. Zhang and H. Sun,yzhang264@ua.edu;
shg@hhu.edu.cn

Citation:Zhang, Y., van Genuchten, M. Th., Zhou,
D., Zhang, G. J., & Sun, H. (2024). A
unified phenomenological model captureswater equilibrium and kinetic processes insoil. Water Resources Research, 60,
e2023WR035782. https://doi.org/10.1029/
2023WR035782
Received 8 JULY 2023
Accepted 26 FEB 2024
Author Contributions:Conceptualization: Yong Zhang,Martinus Th. van Genuchten
Data curation: Yong Zhang,
Dongbao Zhou, Golden J. Zhang
Formal analysis: Yong Zhang,
Dongbao Zhou, Golden J. Zhang
Funding acquisition: Yong Zhang,
HongGuang Sun
Investigation: Yong Zhang, Golden
J. Zhang
Methodology: Yong Zhang, Golden
J. Zhang
Project administration: Yong Zhang,
HongGuang Sun© 2024. The Authors.
This is an open access article under theterms of the Creative Commons
Attribution‐NonCommercial‐NoDerivs
License, which permits use anddistribution in any medium, provided theoriginal work is properly cited, the use isnon‐commercial and no modifications or
adaptations are made.

ZHANG ET AL.

A Unified Phenomenological Model Captures Water
Equilibrium and Kinetic Processes in SoilYong Zhang1 , Martinus Th. van Genuchten2,3, Dongbao Zhou4, Golden J. Zhang5, andHongGuang Sun6
1

Department of Geological Sciences, University of Alabama, Tuscaloosa, AL, USA, 2Department of Nuclear Engineering,
LASME, Federal University of Rio de Janeiro, UFRJ, Rio de Janeiro, Brazil, 3Department of Earth Sciences, Utrecht
University, Utrecht, The Netherlands, 4School of Architecture and Civil Engineering, Anhui Polytechnic University, Wuhu,
China, 5College of Engineering, Texas A&M University, College Station, TX, USA, 6The National Key Laboratory of
Water Disaster Prevention, College of Mechanics and Materials, Hohai University, Nanjing, China

Abstract Soil water sustains life on Earth, and how to quantify water equilibrium and kinetics in soil
remains a challenge for over a century despite significant efforts. For example, various models were proposed tointerpret non‐Darcian flow in saturated soils, but none of them can capture the full range of non‐Darcian flow.
To unify the different models into one overall framework and improve them if needed, this technical note
proposes a theory based on the tempered stable density (TSD) assumption for the soil‐hydraulic property
distribution, recognizing that the underlying hydrologic processes all occur in the same, albeit very complex andnot measurable at all the relevant scales, soil‐water system. The TSD assumption forms a unified fractional‐
derivative equation (FDE) using subordination. Preliminary applications show that simplified FDEs, withproposed hydrological interpretations and TSD distributed properties, effectively capture core equilibrium andkinetic water processes, spanning non‐Darcian flow, water retention, moisture movement, infiltration, andwetting/drying, in the soil‐water system with various degrees and scales of system heterogeneity. Modelcomparisons and evaluations suggest that the TSD may serve as a unified density for the properties of a broad
range of soil‐water systems, driving multi‐rate mass, momentum, and energy equilibrium/kinetic processes
often oversimplified by classical models as single‐rate processes.

Plain Language Summary This technical note introduced a framework aimed at unifying andenhancing existing theories/models for interpreting various soil‐water processes in soils. A unified physical law,termed the “tempered stable density law,” was assumed for soil‐water properties, leading to a unified model
which well captures the core soil water processes that have been modeled (many times inefficiently) by tens ofdifferent and sometimes competing theories/models over the past century. This unified theory/model has thepotential for expansion to other land‐surface equilibrium and kinetic processes that cover a wide range ofspatiotemporal scales and share a similar dynamic nature as water in soils.

1. Introduction
The soil‐water system sustains terrestrial ecosystem functions (Knapp et al., 2008), and many physical laws andmodels have been proposed to quantify water equilibrium and kinetics for over a century. For example, various
empirical and mechanistic models have been developed and are used routinely to describe moisture dynamics orequilibrium states in unsaturated soils, including unsaturated flow (Richards, 1931), soil water retention (Brooks& Corey, 1964; van Genuchten, 1980), and transient infiltration (Horton, 1933; Lewis, 1937). Some of thesemodels, however, exhibit limitations in real‐world applications (Sun et al., 2013). Modeling saturated flow,
particularly non‐Darcian flow, also remains challenging, as demonstrated below. This technical note aims todevelop one theory/model to unify and, if necessary, improve these competing laws/models to interpret core soil‐
water processes, recognizing their occurrence in the same, although very complex, soil‐water system.
The core assumption underlying this unified framework is that soil‐water properties follow the tempered stabledensity (TSD) distribution. Unlike the standard stable density, which is a heavy‐tailed (power‐law) distribution
with desirable mathematical properties (such as being applicable to all scales (called “infinitely divisible” inmathematics) (Samorodnitsky & Taqqu, 1994)), TSD truncates the standard stable density to allow a broader
range of probability density functions (PDFs), including power‐law, exponential, and various intermediate states.
1 of 14

Resources: Yong Zhang
Software: Yong Zhang, Dongbao Zhou,
Golden J. Zhang
Supervision: Yong Zhang,
HongGuang Sun
Validation: Yong Zhang, Dongbao Zhou,
Golden J. Zhang
Writing – original draft: Yong Zhang,Golden J. Zhang
Writing – review & editing: Martinus
Th. van Genuchten, HongGuang Sun

10.1029/2023WR035782

This property makes TSD a common distribution for many hydrologic properties (Cvetkovic, 2011; Zhang
et al., 2015).
In the following sections of this technical note, we (a) derive the unified model using the TSD assumption, (b)apply/validate the model by quantifying core equilibrium and kinetic processes for water in soil, and (c) identify
TSD‐distributed soil‐water properties driving these diverse hydrologic processes. Instead of relying on tensof classical laws/models (many of which have limitations), this technical note introduces a unified model preliminarily validated for core hydrologic processes in the soil‐water system. Furthermore, it has the potential for extension to other equilibrium/kinetic processes in the Earth system, as identified in Zhang, Sun,et al. (2017).

2. Methodology Development
We show that assuming a TSD distribution of soil‐material properties, combined with the time subordination
approach (explained below), leads to a unified fractional‐derivative equation (FDE). Standard FDEs provide a
promising tool to model non‐Fickian transport of pollutants in heterogeneous aquifers (Zhang, Sun, et al., 2017)and other natural and social systems (Sun, Zhang, Baleanu, et al., 2018). The unified FDE extends the capabilitiesof standard FDEs by capturing stochastic processes in soil‐water systems, as illustrated in this technical note.
Time subordination is a probability tool for converting physical time (t) into operational time (τ, which represents
the period that the system is in its operating condition) for moving targets (Baeumer et al., 2001; Zhang, Baeumer,
et al., 2017). Here, the operational time denotes the actual time that a water package would have spent in a
dynamic state had it moved or transitioned at the mean rate, resulting in a random operational time for water
packages driven by system heterogeneity. By adopting the TSD assumption, the governing equation for the
resulting operational time density p in time‐subordination takes the form (Zhang, Baeumer, et al., 2017):
β∗∂p(t,τ) ∂ 1 γ,λ ∂p(t,τ)= 1 γ ,λ
,
∂t
∂t
∂τ

(1)

where the subordinator corresponds to the TSD, β* [Tγ 1] is a scale factor, and the symbol ∂1 γ,λ/∂t1 γ,λ denotes
the (exponentially tempered) Caputo fractional derivative defined by Meerschaert et al. (2008) with index 1 γ[ ] and truncation parameter λ [T 1]. Combing Equation 1 with mass conservation, the time‐subordinated water
dynamics has a density function f that solves the following FDE (without source/sink terms):β∗∂f (x,t)∂ 1 γ ,λ= [A] 1 γ,λ f (x,t),
∂t
∂t

(2)

where β* = 1 is used for unit conversion (no longer shown hereafter for description simplicity), and [A] is anoperator (see examples shown in Table 1). Choosing β* = 1 is suitable for all applications in this study; however,its validity for other conditions/systems requires further confirmation. By rearranging the order of the fractional
derivative in Equation 2, we obtain the equation governing the spatiotemporal evolution of f:∂ γ,λ f (x,t)= [A] f (x,t),
∂tγ,λ

(3)

where Γ(·) is the gamma function.
The first term ∂γ,λf/∂tγ,λ in Equation 3 describes the memory impact on f: the status of f at the present time andlocation, denoted as f(x,t), is affected by the historical condition f(x,τ) (where 0 < τ ≤ t) at the same location whose
weight declines as a TSD function in time (whose power‐law exponent is γ). This time‐nonlocal memory impact
exponentially declines for long historic times (when t τ ≫ 1/λ); therefore, the truncation parameter λ controls
the maximum range of this memory impact. The second term [A]f(x,t) in Equation 3 describes the spatial dynamics of function f when [A]f(x,t) changes with f kinetically, or the rate of change of f when [A]f(x,t) defines anequilibrium constant.
In relatively homogeneous systems with local water dynamics uninfluenced by historical events or long‐range
neighbors, the index γ approaches 1 and Equation 3 simplifies to classical models featuring integer‐order
ZHANG ET AL.

2 of 14

19447973, 2024, 3, Downloaded from https://agupubs.onlinelibrary.wiley.com/doi/10.1029/2023WR035782 by Martinus Genuchten - Cochrane Netherlands , Wiley Online Library on [18/03/2024]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License

Water Resources Research

ZHANG ET AL.
Effective
saturation

Temperature

Heat flow in soil
(Zhang et al., 2021)

Wetting/drying:
Equation S12 in
Supporting Information S1

Chemicalconcentration

Pollutant transport(Zhang, Sun, et al., 2017)

Momentum

Turbulent flux (Sun,Zhang, Wei, et al., 2018)

Infiltration
rate

Sediment
concentration

Sediment suspension(Chen et al., 2013)

Water infiltration:
Equation 15

Solid‐phase
concentration

Adsorption isotherm(Zhou et al., 2021)

Water content

Effective
saturation

Water retention:
Equation 12

Moisture movement:
Equation 13

Darcy
velocity

Function: f

Non‐Darcian flow:
Equation 5

Water‐related
processes

Time

Time

Time

Time

Time

Verticaldistance

Water depth

Aqueous‐phase concentration

Pressure head

Pressure
gradient

Independent variable: g

Wetting anddrying rates

Convection & conduction

Advection & dispersion

Maximumchange of rate

(Net)
Water flux

Scaled shear
velocity

Suspension
index

Rate

Rate

Hydraulicconductivity

Operator: [A]Soil surface area &pore connectivity/size
distribution

Multi‐rate thermal
exchange

Multi‐rate mass
exchange

Multi‐rate decay due tosoil texture and soil
collides/air

Statistics of water
residence time and soil
hydraulic conductivity

Vertical velocity

Turbulent (bursting)Multi‐rate
sorption/desorption

Soil pore size
distribution

Water retention & flow paths

TSD property

Not available

Fourier's law: missing non‐Fourier heat transport

Fick's law: cannot capture non‐Fickian diffusion of pollutants

Horton, Philip, and Lewis‐
Kostiakov models: missing thefull evolution of infiltration
rate

Richards equation: cannot
capture non‐Boltzmann scaling

Newton's law: fail to capture
non‐Newtonian dynamics

Rouse equation: missing thespace nonlocal transport

Langmuir model & Freundlichmodel: limited to monolayer or
multi‐layer, power‐law rate
adsorption (neglecting all theintermediate states)

Brooks‐Corey model:
overestimating the pressure
head near saturation (i.e.,
missing the SWRC transition)

Izbash's law, Forchheimer
model, & Swartzendruber
model: missing the full
evolution of non‐Darcian flow

Standard laws/models: Limitation

10.1029/2023WR035782

3 of 14

19447973, 2024, 3, Downloaded from https://agupubs.onlinelibrary.wiley.com/doi/10.1029/2023WR035782 by Martinus Genuchten - Cochrane Netherlands , Wiley Online Library on [18/03/2024]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License

Note. The parameters f, g, and [A] denote the function, independent variable, and operator in Equation 4, respectively. Additionally, sediment, pollutants, turbulent, and thermal equilibrium or transport
processes are included to illustrate the extension of FDE 4 for quantifying mass, momentum, and energy transport in geologic media, including the soil‐water system. The last column lists the classicallaw/model(s) for each water‐related process along with their intrinsic limitations, which can be addressed by FDE 4.

Non‐stable
Kinetics

Stable kinetics

Equilibrium

Dynamics

Table 1
Applications of the Unified FDE 4 to Ten Water‐Related Processes in the Soil‐Water System

Water Resources Research10.1029/2023WR035782

derivatives like Darcy's law or Richards' equation, as discussed below. For fractional derivatives with γ ∈ (0 1), a
smaller γ can capture a stronger nonlocal impact, particularly in more heterogeneous soil‐water systems (deviating further from the local case with γ = 1). This can be observed in saturated high‐permeability media where
random flow paths induce a pronounced spatially nonlocal impact on non‐Darcian flow (refer to Text S1.2 in
Supporting Information S1), and in unsaturated soils with a wider pore size distribution where the soil water
retention curve (SWRC) slope undergoes a more gradual change in the capillary regime (refer to Section 3.2).
When γ ∈ (1 2), a larger γ (deviating further from the case of γ = 1) can capture a stronger nonlocal impact(especially in time), exemplified in Application #1 (Section 3.1) for low‐velocity non‐Darcian flow, where a TSD
memory function is linked to water retention in low‐permeability media.
By adopting the variable index defined by Sun et al. (2019) and generalizing the independent variable, we obtainthe unified FDE:
∂ γ(g),λ(g) f (g)= [A] f (g),∂gγ(g),λ(g)(4)

where g denotes the independent variable (with f becoming a function of g). After defining f, g, and [A],Equation 4 can model core water equilibrium/kinetic processes in soils described below.

3. Method Applications
3.1. Application #1: Equilibrium Non‐Darcian Flow in Saturated Soil
Non‐Darcian flows, characterized by the non‐linear growth of fluid velocity with respect to the increasingpressure gradient in saturated media, are well‐documented in fluid flow with small velocities in unconsolidated
sand or sedimentary rocks (Nimmo, 2021; Soni et al., 1978) and high velocities in fractured media (Şen, 1987).
Tremendous efforts have been dedicated to quantifying non‐Darcian flow for over a century, resulting in well‐known models, including Izbash's law (Izbash, 1935), the Forchheimer model (Forchheimer, 1901), and the
Swartzendruber model (Swartzendruber, 1962). To unify and enhance these models, we simplify the unified FDE
4 as the Equilibrium‐FDE 5:∂ γ,λq = K,∂J γ,λ(5)where q [L/T] is the fluid velocity, J [ ] is the absolute value of the hydraulic gradient, and K [L/T] is the effective
hydraulic conductivity. Equation 5 can also be derived using structure derivatives, as shown in Text S1.1 in
Supporting Information S1. The analytical solution of Equation 5 is:
q = K Eγ1,γ+1 ( λJ) J γ ,(6)∞Γ(γ) zwhere Eγ1,γ+1 (·) denotes the generalized Mittag‐Leffler function: Eγη,μ (z) = ∑ Γ(γ+k)/Γ(ηk+μ)k! . We note thatkk=0

Mainardi (2020) defined the Mittag‐Leffler function as the “queen function” of fractional calculus since it
continuously interpolates between purely exponential law and power‐law‐like functions; those properties arehighly relevant to our study. FDE 5 can also be applied to other equilibrium pairs in hydrology and geoengineering, such as the asymptotic vertical distribution of suspended sediment in rivers/streams (Chenet al., 2013) and the creep of frozen soil (Gu et al., 2023). Its applications in other equilibrium states, including
total or component mass balance, summation, and energy balance (Cantele et al., 2021; Rodriguez‐Robles
et al., 2019), remain to be explored.
Equation 5 uses the truncation parameter λ to regulate the transition from non‐Darcian flow to Darcian or other
non‐Darcian flows. Such a broad transition was overlooked by existing non‐Darcian models. For example, when
λ = 0 (neglecting state transition), Equation 5 reduces to ∂γq/∂Jγ = K, and its solution, Equation 6, reduces to theempirical Izbash's law (Benedikt et al., 2018):

ZHANG ET AL.

4 of 14

19447973, 2024, 3, Downloaded from https://agupubs.onlinelibrary.wiley.com/doi/10.1029/2023WR035782 by Martinus Genuchten - Cochrane Netherlands , Wiley Online Library on [18/03/2024]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License

Water Resources Research10.1029/2023WR035782

Figure 1. Application #1: Non‐Darcian flow of an organic fluid (Soltrol‐130) through sandstone in Core 1 (a), Core 2 (b),Core 3 (c), and Core 4 (d), respectively. The symbols denote measurements by Siddiqui et al. (2016), and the lines represent
the best‐fit solutions using different models, with the best‐fit parameters listed in Table S1 in Supporting Information S1.

q = KJ γ .

(7)Izbash's law 7 captures low‐velocity non‐Darcian flow but is unable to model the transition between Darcian andnon‐Darcian flows, or that between various non‐Darcian flows. To incorporate both Darcian and non‐Darcian
flows, the Forchheimer model (derived from the Navier‐Stokes equation or the capillary tube approach) takes
the form (Takhanov, 2011):J=

μq + ηρq2 ,k(8)

where μ [Pa·T] is the fluid dynamic viscosity, k [L2] is the medium permeability, η [L 1] is the inertial resistance,and ρ [M·L 3] is the fluid density. On the right‐hand side (RHS) of the Forchheimer model 8, the 1st term captures
Darcian flow affected by viscous force at low velocities, and the 2nd term describes quadratic flow due to inertial
resistance at large pressure drops and high velocities (such as flow in the vicinity of a pumping well). Hence,
Equation 8 reduces to Darcy's law at low velocities, missing (a) low‐velocity non‐Darcian flow where the q(J)curve has an initial power‐law slope steeper than 1, and (b) flow in low‐permeability media where the q(J) curveis concave up in a linear‐linear plot, requiring a negative β in Equation 8. To efficiently capture widely observed
low‐velocity non‐Darcian flow, the Swartzendruber model takes the form (Swartzendruber, 1962):
q = K [JZHANG ET AL.

J∗ 1

∗

e J/J )],

(9)

5 of 14

19447973, 2024, 3, Downloaded from https://agupubs.onlinelibrary.wiley.com/doi/10.1029/2023WR035782 by Martinus Genuchten - Cochrane Netherlands , Wiley Online Library on [18/03/2024]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License

Water Resources Research10.1029/2023WR035782

where J* represents the threshold hydraulic gradient. In Equation 9, the change of fluid velocity with respect tothe pressure gradient increases exponentially (whose rate is controlled by the threshold J*) before approaching its
asymptote K at J ≫ J* (i.e., transition to Darcian flow). Therefore, the Swartzendruber model 9 misses high‐velocity non‐Darcian flow.
These models are checked against real‐world data, with one representative result plotted in Figure 1. Additional
details and other cases are presented in Text S1.2 in Supporting Information S1. These applications support theanalysis mentioned above: Izbash's law 7 fails to capture the velocity‐dependent slope of the observed q(J)curves, the Forchheimer model 8 misses low‐velocity non‐Darcian flow, the Swartzendruber model 9 overlooks
high‐velocity non‐Darcian flow, and the FDE 5 emerges as the only viable option capable of capturing the full‐range of non‐Darcian flow.
Non‐Darcian flow has two main triggers. Trigger #1 involves microscopic inertial effects (Forchheimer, 1901;
Ruth & Ma, 1992) and increasing microscopic viscous force (Hassanizadeh & Gray, 1987), generating non‐Darcian flow at high velocities. Trigger #2 includes heterogeneities in sample permeability or flow velocity,
flow streamlining, and pore participation number (Barree & Conway, 2005), non‐Newtonian flow (Kutiĺek, 1972), or boundary layer effects (Wang et al., 2018), causing low‐velocity non‐Darcian flow. Here, weinterpret non‐Darcian flow using the FDE 5. For low‐permeability materials with usually low flow rates, the best‐
fit index γ in FDE 5 is between 1 and 2 (Text S1.2 in Supporting Information S1), where the fractional derivative
can be expanded as:∞ 2
d γ ,λd2 q
1d q λτ
q = 2 ∗ g(J) =
e ν(J
∫
γ,λ
Γ(2 γ) J dτ2dJdJτ)1 γ dτ,(10)

where the symbol “∗” denotes convolution, and g(J) = Γ(21 γ) ∫J ∞ e λτ γ τ1 γ dτ is the TSD memory function.
Equation 10 can be approximated using the one‐shift Grünwald formula for λ → 0 (Meerschaert &Tadjeran, 2004):
dγ,λ q1 N≈∑ fi q[Jγ ,λ(∆J)γ i=0dJ

(i1)∆J] =

1
{q(J + ∆J)
(∆J)γ

N

γq(J) + ∑ fi q[J(i1)∆J] }

(11)i=2

where ∆J is the step size discretizing the pressure gradient, N is a sufficiently large number of grid points of theγ)pressure gradient, and fi = Γ( Γ(iγ)Γ(i+1) is the Grünwald weight (the summation of these weights equals 1 for i = 0,

1,2, ⋯, and fi > 0 for i ≥ 2 when 1 < γ < 2) whose absolute value decreases with an increasing sequence i.
Equation 10 implies that the change of q with J is a nonlocal process. Water moving in low‐permeability soil canbe retarded by fluid‐rock interactions (Oelkers & Schott, 2018), sorption‐desorption at the mineral/water interface(Scheidegger & Sparks, 1996), and mass exchange between the main flow zone and the surrounding less‐mobile
regions or between soil/aquifer zones (Eberts & George, 2000). These less mobile water packages can bereactivated with a remobilization probability proportional to the Grünwald weight |fi| as shown in Equation 11.
This nonlocal impact is described by the summation of all the “less‐mobile water” terms marked by the sequencei ≥ 2 (all the terms except γq(J) on the RHS of Equation 11 are positive, because fi > 0 for i > 2), explaining thepositive concavity of the q(J) curve in a linear‐linear plot, as shown, for example, in Figure S1 in Supporting
Information S1.
In the FDE 5, the index γ defines the overall degree of the nonlocal impact, where a larger γ captures a stronger
nonlocal impact (represented by a relatively larger positive concavity) and a lager power‐law slope (deviating
more from linear flow) of the q(J) curve (e.g., Table S2 and Figure S1 in Supporting Information S1). In addition,
the truncation parameter λ controls the transition from pre‐Darcian to Darcian flow: low‐velocity non‐Darcian
flow occurs for the hydraulic gradient J ≪ 1/λ, gradually shifting to Darcian flow with increasing J ≫ 1/λ.
For the special case of γ = 1, all the “historical” components in Equation 11 disappear (i.e., without the “nonlocal”
impact), and Equation 11 reduces to the forward first‐order difference which defines Darcy flow, as expected.
Therefore, an FDE 5 with an index γ between 1 and 2 characterizes low‐velocity non‐Darcian flow, likely due tothe nonlocal impact with a TSD memory function related to water retention in low‐permeability media. Similarly,
high‐velocity non‐Darcian flow may be due to the spatially nonlocal impact related to random flow paths with a

ZHANG ET AL.

6 of 14

19447973, 2024, 3, Downloaded from https://agupubs.onlinelibrary.wiley.com/doi/10.1029/2023WR035782 by Martinus Genuchten - Cochrane Netherlands , Wiley Online Library on [18/03/2024]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License

Water Resources Research10.1029/2023WR035782

Figure 2. Application #2: Measured (symbols, from Novák and Hlaváčiková (2019)) versus best‐fit model results for the
SWRC of Soils 1 through 7 in (a) and (b). The optimized parameters are listed in (c). The Se(h*) curve for all 7 soils (c), with
Soil 7 plotted again in (d) to show the power‐law slope. Plot (e) shows the specific soil water capacity curve for these 7 soils,
with Soil 7 plotted again in (f) to show the power‐law slope.

TSD distribution in high‐permeability media, and it can be captured by Equation 5 with a γ between 0 and 1 (seeText S1.3 in Supporting Information S1). Therefore, the TSD memory function due to soil‐water heterogeneity
may drive both low‐ and high‐velocity non‐Darcian flows and the subsequent transition, expanding on trigger #2
mentioned above for non‐Darcian flow.
3.2. Application #2: Equilibrium Water Retention in Unsaturated Soil
The unified FDE 4 also links water content θ and pressure head h in unsaturated soil:∂ γ,λSe = k,∂(h∗ )γ,λ(12)where Se = (θ θr)/(θs θr) [ ] is the effective saturation; θ, θr, and θs [ ] are the volumetric, residual, andsaturated water contents, respectively; h* = |hd/h| [ ], with hd [L] representing the air entry value; and the rate k⃒⃒ 1[ ] can be approximated by k ≈ [⃒hd |γ Eγ1,γ+1
λ |hd ⃒)] , or simply calculated by k = Γ(γ + 1) if the index
0 < γ < 1 or λ is small (e.g., λ < 1). The Equilibrium‐FDE 12 (the same form as Equation 5) models a temporally
nonlocal change of water content with the pressure head in multiple rates whose distribution follows the TSD.
This results in an upper‐truncated power‐law decline of Se at large h*, with the corresponding soil water retention
curve (SWRC) exhibiting the shape of a sigmoid curve at small |h|.
One application is depicted in Figures 2a and 2b. Novák and Hlaváčiková (2019) collected the SWRC drainage
branches of typical soils, including nearly monodispersed glass sand (Soil 1), a sandy soil (Soil 3), and a fine‐textured (clay) soil (Soil 7). The index γ is approximated by the log‐log slope of the initial rising limb of theobserved Se(h*) curve (Figure 2c), consistent with Equation 12. This power‐law portion of the Se(h*) curve covers
ZHANG ET AL.

7 of 14

19447973, 2024, 3, Downloaded from https://agupubs.onlinelibrary.wiley.com/doi/10.1029/2023WR035782 by Martinus Genuchten - Cochrane Netherlands , Wiley Online Library on [18/03/2024]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License

Water Resources Research10.1029/2023WR035782

the capillary and adsorption regions, which involve water flow in small pores. The truncation parameter λ, which
captures the upper bound of the heavy‐tailed power‐law function Se(h*), can be approximated by the cross point of
the two power‐law sections of Se(h*) (Figure 2d). Model comparisons in Text S2.1 in Supporting Information S1
show that FDE 12 is closely related to the van Genuchten model, the most reliable model for fitting various
SWRCs (the two models provide almost identical fits), and it improves the Brooks‐Corey model by capturing thetransition of the SWRC slope (Figures 2a and 2b). The close fit observed between FDE 12 and the van Genuchten
model can be attributed to the algebraic relationship of parameters shared between these two models, as illustrated
in Text S2.1 in Supporting Information S1. This connection suggests that the FDE may provide a potential method
for associating these model parameters with soil‐water properties.
The physical meanings of model parameters γ and λ can be interpreted using the concept of capillary tubes.
Equation 12 assumes the TSD distribution of the soil pore size r (usually within the range of 50–0.5 mm), denoted
as f(r) here. According to cylindrical capillary bundle models (e.g., Millington & Quirk, 1961; Novák & Hlaváčiková, 2019), all capillary tubes of radius r > Rj are drained at a given hj (following the capillarity equation
Rj ∝ 1/hj). Since large pores are known to fill or empty at small pressure heads (|h|→0), a PDF f(r) with manylarge pores results in the SWRC dropping quickly at small values of |h|. The section of f(r) for fine pores, on theother hand, controls the SWRC slope, with a wider pore size distribution leading to a more graduate change in the
SWRC slope mainly in the capillary regime (Nimmo, 2004), which can be captured by a smaller index γ in
Equation 12 (see Figure 2c). While the above‐mentioned description ignores the effects of film, corner, and vapor
flow in the very dry range (e.g., Iden et al., 2021), the macroscopic form is realistic for many soils. Therefore, thePDF f(r) for fine pores defines the range of γ, and the transition from fine to large pores in f(r) defines the valueof λ.
We further correlate the soil pore size distribution and the FDE index γ. The specific soil water capacity, c = dθ/dhof a drying θ(h) curve, approximates the PDF of the combined volume of all pores whose effective radius iscentered at r (Figure 2e). For example, Soil 1 (glass sand) had a narrow f(r) (Novák & Hlaváčiková, 2019) (see theline “1” in Figure 2e), resulting in a relatively flat SWRC (like a step function shown by curve “1” in Figure 2a)which can be captured by Equation 12 with the largest index γ and the smallest truncation parameter λ of the sevensoils (listed in Figure 2c). Soil 7 (clay), on the other hand, had a relatively broad f(r) (Figure 2e), whose muchsteeper SWRC (see Figure 2a) can be captured by Equation 12 using small γ and large λ values. Additional
analyses in Text S2.2 in Supporting Information S1 also show that the power‐law slope of the pore volume
distribution is linearly related to the FDE index γ.
3.3. Application #3: Soil Moisture Movement Kinetics
Unsaturated flow can exhibit strong variations in both space and time, involving more soil‐water properties thanequilibrium water retention. When f = θ, g = t, and [A] denotes the differential Buckingham‐Darcy flux, theunified FDE 4 reduces to the Flow‐FDE 13, quantifying unsaturated flow:e λt∂ γ (eλ t θ)
∂
∂θ
∂
∂z
=
[Dw (θ)
]+
[K(θ)
],
∂tγ
∂xi
∂xi
∂xi
∂xi

(13)

where xi [L] denotes the direction of water movement, Dw(θ) = K(θ)dh/dθ [L2/T] is the soil water diffusivity, andK(θ) [L/T] is the unsaturated soil hydraulic conductivity. When λ→0, the Flow‐FDE 13 reduces to the standard
time‐fractional Richards equation as shown by Gerolymatou et al. (2006). When γ = 1 and λ → 0, the Flow‐FDE
13 further reduces to the classical Richards equation (Richards, 1931):
∂θ
∂
∂θ
∂
∂z
=
[Dw (θ)]+[K(θ)].
∂t ∂xi
∂xi
∂xi
∂xi

(14)

Hydrological interpretation of the Flow‐FDE 13 may be related to the multi‐rate, or so‐called “anomalous,” watertransport in unsaturated media proposed by several authors (Gerolymatou et al., 2006; Sun et al., 2013). Water
moves in soil with both fast flow (via local preferential pathways such as structural voids or biological channels)
and slow components (e.g., moving along sinuous paths, trapped by relatively low‐permeability or immobile
zones, or delayed by mass exchange between mobile and stagnant regions), generating a random operational time

ZHANG ET AL.

8 of 14

19447973, 2024, 3, Downloaded from https://agupubs.onlinelibrary.wiley.com/doi/10.1029/2023WR035782 by Martinus Genuchten - Cochrane Netherlands , Wiley Online Library on [18/03/2024]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License

Water Resources Research10.1029/2023WR035782

Figure 3. Application #3: Non‐Boltzmann scaling of wetting fronts (horizontal infiltration). (a) Experimental wetting front
data (El‐Abd & Milczarek, 2004) in a fired‐clay brick, with fitted modeling results using a power‐law diffusivity.
Corresponding times for the five sets of data are t = 419, 2219, 14879, 26099, and 34559 s. The optimized diffusivity
function for the Richards Equation 14 was Dw(θ) = 7.5 × 10 2θ1.75. For the Flow‐FDE 13, the best‐fit data wereDw(θ) = 1.9 × 10 2θ2.25 mm2s γ, γ/2 = 0.60 (γ/2 > 0.5 in the FDE for super‐diffusion), and λτ = 0. (b) Experimental data andmodel fits for a siliceous brick. The times for the five distributions are t = 450, 5370, 24210, 79770, and 170430 s. The
solution curves for Flow‐FDE were obtained using an exponential diffusivity Dw(θ) = 1.78 × 10 4exp(8.4θ) mm2s γ withbest‐fit γ/2 = 0.43 (γ/2 < 0.5 for sub‐diffusion) and λτ = 1.0 × 10 5 s 1. The water content curves for the Richards equation
were reproduced from Figure 9 in Gerolymatou et al. (2006).

(distributed as TSD) for water packages and therefore “anomalous” dynamics. For example, when λ = 0 and Dwis
constant in the FDE 13, the spatial variance of horizontal moisture movement increases nonlinearly in time
(Zhang et al., 2008): σ2 (t) = 2Dw tγ / Γ(γ + 1), which captures super‐diffusion with γ/2 > 0.5 and sub‐diffusion
with γ/2 < 0.5.
Two applications for horizontal water flow are shown in Figure 3. In such instances, solving the nonlinear Flow‐
FDE 13 requires the use of a numerical solver (developed in Text S5 in Supporting Information S1). Results
showed that the classical Richards Equation 14 either underestimated the super‐diffusive water content front(Figure 3a) or overestimated the sub‐diffusive water content front (Figure 3b). This is because Equation 14 assumes that the horizontal travel distance of moisture grows as the square root of time (i.e., Boltzmann scaling),
while real‐world horizontal flow often exhibits non‐Boltzmann or “anomalous” scaling (e.g., Sun et al., 2013)
described by the Flow‐FDE 13. Notably, the non‐Boltzmann scaling can also be effectively quantified through thefractal Richards equation, as demonstrated by Sun et al. (2013). A concise comparison between the fractional andfractal Richards equations is provided in Text S5 in Supporting Information S1. Additionally, the fractional
Richards Equation 13 can be linearized using the TSL and full subordination, as detailed in Text S6 in Supporting
Information S1.
Saturated flow dynamics, the core process in hydrogeology studies, may have similar mechanisms. For example,
transient flow in confined aquifers is “anomalous” and can be captured by a simplified Flow‐FDE 13 by assuming
λ = 0 and replacing the water content with the hydraulic head in Equation 13 (Xia et al., 2021).
3.4. Application #4: Water Infiltration Kinetics
The water infiltration rate i(t), especially during ponding, into a soil usually declines with a variable deceleration
rate in time (Jury & Horton, 2004). It is well‐known that the infiltration rate is high at the beginning of a storm or
ponded irrigation (denoted as i0 [L/T]) when the soil is dry, and then decreases as the soil becomes saturated,
ultimately approaching an asymptote value if at or close to the profile‐averaged saturated hydraulic conductivity.
Many empirical models exist for describing water infiltration into soils, including the widely applied Horton,
Green‐Ampt, and Philip equations (Philip, 1957; Rahmati et al., 2022). To unify and improve these models, we
simplify the unified FDE 4 to the Kinetic‐FDE 15:∂ γ ,λf = i0∂tγ,λif ) f ,with f = i0‐i(t), g = t, and [A] denoting the fluctuation of the infiltration rate i0

ZHANG ET AL.

(15)
if. The solution is:

9 of 14

19447973, 2024, 3, Downloaded from https://agupubs.onlinelibrary.wiley.com/doi/10.1029/2023WR035782 by Martinus Genuchten - Cochrane Netherlands , Wiley Online Library on [18/03/2024]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License

Water Resources Research10.1029/2023WR035782

Figure 4. Application #4: Change in infiltration rate with time. Measured data (symbols) versus the model fit (lines) for sandy
loam (a), loam (b), and clay loam (c). The right plot is a log‐log version of the left plot, emphasizing early‐ and late‐timetailing behavior.

i(t) = i0

i0i f ) [t1 α E1,α2 α ( λt)λα t].

(16)

The Kinetic‐FDE solution 16 captures various evolutions of the infiltration rate, including a power‐law function,
an exponential function, and a constant asymptote. Equation 16 for this purpose contains the Horton and Philipequations as end members. For example, when λ = 0, Equation 16 reduces to the power‐law form i(t) ∝ t1 α,which is analogous to the Lewis‐Kostiakov model (Lewis, 1937):
i(t) = b tc 1 ,(17)where b and c are empirical constants depending on soil type, and the Philip model for vertical infiltration:Si(t) = √̅ + A,
2 t

(18)

over a relatively short time (where S and A are constants). When λ is large, the Kinetic‐FDE solution 16 reduces toa single infiltration rate model analogous to the exponential Horton model 19:
i(t) = i f + i0

i f ) e εt ,(19)

where ε [hr 1] is an empirical parameter describing the rate of decrease of infiltration.
We tested these models by fitting real‐world data (Figure 4). Results revealed that the Kinetic‐FDE solution 16
captures the observed transient infiltration rate slightly better than the other models, with the Lewis model 17 andHorton model 19 producing correct solutions only for early times (t ≪ 1/λ) and late times (t ≫ 1/λ), respectively.
This is because the solution 16 describes an exponentially truncated power‐law function using the power‐law
portion for times t ≪ 1/λ and the exponential portion for times t ≫ 1/λ (we recall here that the Mittag‐Leffler
function can be used to interpolate continuously between Gaussian and Lorentzian functions). By comparison,

ZHANG ET AL.

10 of 14

19447973, 2024, 3, Downloaded from https://agupubs.onlinelibrary.wiley.com/doi/10.1029/2023WR035782 by Martinus Genuchten - Cochrane Netherlands , Wiley Online Library on [18/03/2024]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License

Water Resources Research10.1029/2023WR035782

the Lewis and Horton models are based on a power‐law function and an exponential function, respectively, for thetransient infiltration rate process. In addition, the Philip model 18 assumes a constant power‐law exponent (0.5)
for the full evolution of the infiltration rate for all soils, which is clearly not consistent with all the data.

4. Discussion and Conclusion
4.1. TSD of Soil‐Hydraulic Properties
The TSD assumption underlying the FDE 4 is analogous to the fractal morphology assumed in fractal models forsoil structure. Various models have been proposed to relate the empirical SWRC fitting parameters to the fractal
dimension of soil pore or particle size distributions (e.g., Pachepsky et al., 1995; Zhao, 2020). Soil pores display
self‐similar scaling with an infinite nesting structure of micropores, macropores, and mesopores. The internal
structure, pore size distribution, and pore geometry of soil exhibit fractal features, which can be characterized bydifferent fractal dimensions (Giménez et al., 1997). Independent studies have also demonstrated a close relationship between the fractal dimension and the FDE index (Verma & Viswanathan, 2020). The FDE 4 offers
greater flexibility than standard fractal models in applications, since it allows (a) a single fractional index toencompass multiple fractal dimensions or intrinsic dimensions of soil (Giménez et al., 1997), and (b) captures theimpact of hydraulic properties on water processes. In addition, the truncation parameter λ defines the upper or, insome cases, lower limit of the power‐law distribution of soil and hydraulic properties listed in Table 1. This
truncation is necessary as natural processes are generally bounded, distinguishing TSD from the standard stabledensity.
Significantly, it is essential to recognize that the TSD‐based FDE 4 is not designed as a predictive model,
constituting the principal limitation of this note. Further discussion on this limitation is provided in Text S4 in
Supporting Information S1. Detailed information on the soil‐water properties driving water equilibrium/kinetics(referencing column 6 in Table 1) was not provided by most applications discussed above (explained in Text S2.2
in Supporting Information S1). Further experimental efforts, hence, are needed to detect the distribution of thesoil‐water properties identified theoretically in this note to improve the FDE's predictability. This technical note
focuses on water flow rather than chemical reactions. Chemical reactions in heterogeneous media may involve
properties beyond those related to soil and water. Extending the analysis to include chemical reactions can bechallenging (Bolster et al., 2017) and may necessitate reformulation, particularly in the presence of structural
heterogeneity causing non‐uniform reactivity across exchange rates (Painter, 2021; Roche & Dentz, 2022) or incases where reactions are mixing‐limited (Dentz et al., 2011). Notably, the TSD theory may potentially be
expanded to integrate chemical reactions as an additional term in FDE 4. For instance, the time‐dependentchemical reaction rate, influenced by incomplete mixing of reactants in heterogeneous porous media, may
exhibit a power‐law decline (before reaching its asymptote) (Sanchez‐Vila et al., 2010). This could lead to a
fractional‐order advection‐dispersion‐reaction equation, as proposed by Zhang (2023). We will explore this topic
in future research.
4.2. Model Generalization and Application in Earth Sciences
If a source‐sink term is present (and denoted as r), the FDE 4 can be expressed as:∞∂ γ(g),λ(g) f (g)γ(s) s γ(s) 1
λ(g)g
=
f
+
e
eλ(s) s r(x,s)
ds,
∫
[A]
(g)
γ(g)
,
λ(g)
∂g
Γ(1 γ(s))g(20)

where the last term captures the time‐nonlocal source/sink term (in contrast to the classical, time‐local source/sinkterm). It is important to note that this source‐sink term has not been considered in the applications shown above,and its relevance and formulation necessitate further exploration.
Water flow in soils can be conceptualized as a multi‐rate mass transfer (MRMT) process occurring between
connected pores (containing mobile water) and relatively immobile regions (such as soil aggregates, small or
isolated pores, and interparticle fracture). This water retention/movement in soils is analogous to chemical
sorption/transport in aquifers, which is usually a MRMT process (Zhang, Sun, et al., 2017). The FDE 4 may serve
as a unified model for both processes. Notably, the introduction of chemical heterogeneity into the soil‐water
system leads to more variants in quantification formulas, as evidenced by the 15 chemical adsorption isotherm
ZHANG ET AL.

11 of 14

19447973, 2024, 3, Downloaded from https://agupubs.onlinelibrary.wiley.com/doi/10.1029/2023WR035782 by Martinus Genuchten - Cochrane Netherlands , Wiley Online Library on [18/03/2024]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License

Water Resources Research10.1029/2023WR035782

models reviewed by Foo and Hameed (2010). Therefore, a unified model is needed. As demonstrated by ourrecent work (Zhou et al., 2021), a simplified version of the FDE 4 (see Table 1) outperforms the Langmuir andFreundlich models in capturing various adsorption isotherms of PFAS onto soils. “Multi‐rate equilibrium/kinetics” may for this reason be a common feature in many soil‐water systems, which cannot be reliably captured byclassical models assuming a single‐rate mass transfer.
The FDE 4 with a constant index captures stable equilibrium (exhibiting consistent trends and/or returning to theoriginal status), while the FDE 4 with variable indexes accommodates transient equilibrium, reflecting changing
behaviors in stages or cycles, such as hysteresis. An illustrative example is provided in Text S3 in Supporting
Information S1, where the variable‐index FDE 4 is the only viable model capable of capturing a non‐stable
wetting/drying process. Additionally, the FDE 4 generalizes Newton's law of diffusivity and Fourier's law ofheat conduction, enabling the quantification of momentum and energy transport in geologic media with TSD
properties (Table 1). The FDE 4 reveals transitions in momentum/energy dynamics, controlled by the truncation
parameter, which were overlooked by existing models in Sun, Zhang, Wei, et al. (2018) and Zhang et al. (2021).
Furthermore, it may offer insights into other kinetic processes on land surfaces involving MRMT, such as crystal
growth (which is a nonlocal diffusive‐controlled process whose rates are affected by various factors), rock/mineral weathering (with a spatiotemporal‐scale dependent rate), and pedogenesis (where soil production andtransport rates are affected by interconnected physical and bio‐chemical factors) (Zhang, Sun, et al., 2017).

Conflict of Interest
The authors declare no conflicts of interest relevant to this study.

Data Availability Statement
The data on which this article is based are available in Siddiqui et al. (2016), Novák and Hlaváčiková (2019), andGerolymatou et al. (2006). Our code used for this article is available from Zhang et al. (2023).
Acknowledgments
H.G.S. was partially funded by the
National Natural Science Foundation ofChina (Grants U2267218, 11972148, and41931292). The results of this study do notreflect the view of this finding agency. Y.
Z. was supported partially with federal
funding from the Department of theTreasury under the Resources andEcosystems Sustainability, Tourist
Opportunities, and Revived Economies of
the Gulf Coast States Act of 2012
(RESTORE Act). Y.Z. was also partially
funded by the National Science Foundation(NSF) (Grant 2305141), United States. The
statements, findings, conclusions, andrecommendations are those of the authorand do not necessarily reflect the views of
the Department of the Treasury or ADCNRor NSF.

ZHANG ET AL.

References
Baeumer, B., Benson, D. A., Meerschaert, M. M., & Wheatcraft, S. W. (2001). Subordinated advection‐dispersion equation for contaminanttransport. Water Resources Research, 37(6), 1543–1550. https://doi.org/10.1029/2000wr900409
Barree, R. D., & Conway, M. W. (2005). Reply to discussion of “beyond beta factors: A complete model for Darcy, Forchheimer, and Trans‐
Forchheimer flow in porous media”. Journal of Petroleum Technology, 57(8), 73–74. https://doi.org/10.2118/0805‐0073‐jpt
Benedikt, J., Girg, P., Kotrla, L., & Takac, P. (2018). Origin of the p‐Laplacian and A. Missbach. The Electronic Journal of Differential Equations,16, 1–17.
Bolster, D., Benson, D. A., & Singha, K. (2017). Upscaling chemical reactions in multicontinuum systems: When might time fractional equations
work? Chaos, Solitons and Fractals, 102, 414–425. https://doi.org/10.1016/j.chaos.2017.04.028
Brooks, R., & Corey, A. (1964). Hydraulic properties of porous media. In Hydrology papers. Colorado State University.
Cantele, M., Bal, P., Kompas, T., Hadjikakou, M., & Wintle, B. (2021). Equilibrium modeling for environmental science: Exploring the nexus of
economic systems and environmental change. Earth's Future, 9(9), e2020EF001923. https://doi.org/10.1029/2020ef001923
Chen, D., Sun, H. G., & Zhang, Y. (2013). Fractional dispersion equation for sediment suspension. Journal of Hydrology, 491, 13–22. https://doi.
org/10.1016/j.jhydrol.2013.03.031
Cvetkovic, V. (2011). The tempered one‐sided stable density: A universal model for hydrological transport? Environmental Research Letters,
6(3), 034008. https://doi.org/10.1088/1748‐9326/6/3/034008
Dentz, M., Le Borgne, T., Englert, A., & Bijeljic, B. (2011). Mixing, spreading and reaction in heterogeneous media: A brief review. Journal ofContaminant Hydrology, 120–121, 1–17. https://doi.org/10.1016/j.jconhyd.2010.05.002
Eberts, S. M., & George, L. L. (2000). Regional ground‐water flow and geochemistry in the Midwestern basins and arches aquifer system in partsof Indiana (p. 103). US Department of the Interior, U.S. Geological Survey.
El‐Abd, A. E., & Milczarek, J. J. (2004). Neutron radiography study of water absorption in porous building materials: Anomalous diffusion
analysis. Journal of Physics D: Applied Physics, 37(16), 2305–2313. https://doi.org/10.1088/0022‐3727/37/16/013
Foo, K. Y., & Hameed, B. H. (2010). Insights into the modeling of adsorption isotherm systems. Chemical Engineering Journal, 156(1), 2–10.
https://doi.org/10.1016/j.cej.2009.09.013
Forchheimer, P. (1901). Wasserbewegung durch boden. Zeitschrift des Vereins Deutscher Ingenieure, 45, 1781–1788.
Gerolymatou, E., Vardoulakis, I., & Hilfer, R. (2006). Modelling infiltration by means of a nonlinear fractional diffusion model. Journal of
Physics D: Applied Physics, 39(18), 4104–4110. https://doi.org/10.1088/0022‐3727/39/18/022
Giménez, D., Perfect, E., Rawls, W. J., & Pachepsky, Y. (1997). Fractal models for predicting soil hydraulic properties: A review. Engineering
Geology, 48(3–4), 161–183. https://doi.org/10.1016/s0013‐7952(97)00038‐0
Gu, X. T., Yu, X. N., & Zhang, G. J. (2023). Application of constant‐ and variable‐order fractional and fractal Kelvin models to quantify uniaxial
creep of frozen soil. In prep.
Hassanizadeh, S. M., & Gray, W. G. (1987). High velocity flow in porous media. Transport in Porous Medium, 2(6), 521–531. https://doi.org/10.
1007/bf00192152
Horton, R. E. (1933). The role of infiltration in the hydrologic cycle. EOS Transactions American Geophysical Union, 14(1), 446–460.

12 of 14

19447973, 2024, 3, Downloaded from https://agupubs.onlinelibrary.wiley.com/doi/10.1029/2023WR035782 by Martinus Genuchten - Cochrane Netherlands , Wiley Online Library on [18/03/2024]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License

Water Resources Research10.1029/2023WR035782

Iden, S. C., Diamantopoulos, E., & Durner, W. (2021). Capillary, film, and vapor flow in transient bare soil evaporation (2): Experimental
identification of hydraulic conductivity in the medium to dry moisture range. Water Resources Research, 57(5), e2020WR028514. https://doi.
org/10.1029/2020wr028514
Izbash, S. V. (1935). Construction of Dams by Dumping Stones into Flowing Water. United States Engineer Office, Engineering Division.
Translated by A. Dovjikov. Eastport, Maine. 138 Retrieved from http://resolver.tudelft.nl/uuid:c46a945f‐c0a7‐4d34‐bb3d‐7997e52fc324
Jury, W. A., & Horton, R. (2004). Soil physics (p. 370). John Wiley and Sons, Inc.
Knapp, A. K., Beier, C., Briske, D. D., Classen, A. T., Luo, Y., Reichstein, M., et al. (2008). Consequences of more extreme precipitation regimes
for terrestrial ecosystems. BioScience, 58(9), 811–821. https://doi.org/10.1641/b580908
Kutiĺek, M. (1972). Non‐Darcian flow of water in soils—Laminar region: A review. Developments in Soil Science, 2, 327–340.
Lewis, M. R. (1937). The rate of infiltration of water in irrigation practice. EOS Transactions American Geophysical Union, 18, 361–368.
Mainardi, F. (2020). Why the Mittag‐Leffler function can be considered the queen function of the fractional calculus? Entropy, 22(12), 1359.
https://doi.org/10.3390/e22121359
Meerschaert, M. M., & Tadjeran, C. (2004). Finite difference approximations for fractional advection‐dispersion flow equations. Journal ofComputational and Applied Mathematics, 172(1), 65–77. https://doi.org/10.1016/j.cam.2004.01.033
Meerschaert, M. M., Zhang, Y., & Baeumer, B. (2008). Tempered anomalous diffusion in heterogeneous systems. Geophysical Research Letters,
35(17), L17403. https://doi.org/10.1029/2008gl034899
Millington, R. J., & Quirk, J. P. (1961). Permeability of porous solids. Transactions of the Faraday Society, 57, 1200–1207. https://doi.org/10.
1039/tf9615701200
Nimmo, J. R. (2004). Porosity and pore size distribution. Encyclopedia of Soils in the Environment, 3(1), 295–303. https://doi.org/10.1016/b0‐12‐348530‐4/00404‐5
Nimmo, J. R. (2021). The processes of preferential flow in the unsaturated zone. Soil Science Society of America Journal, 85(1), 1–27. https://doi.
org/10.1002/saj2.20143
Novák, V., & Hlaváčiková, H. (2019). Chapter 7—Soil‐water retention curve. Applied Soil Hydrology, Theory and Applications of Transport in
Porous Media, 32, 77–96. https://doi.org/10.1007/978‐3‐030‐01806‐1_7
Oelkers, E. H., & Schott, J. (2018). In Thermodynamics and kinetics of water‐rock interaction (Vol. 70). Walter de Gruyter GmbH & Co KG.
Pachepsky, Y. A., Polubesova, T. A., Hajnos, M., Sokolowska, Z., & Jozefaciuk, G. (1995). Fractal parameters of pore surface area as influenced
by simulated soil degradation. Soil Science Society of America Journal, 59(1), 68–74. https://doi.org/10.2136/sssaj1995.
03615995005900010010x
Painter, S. L. (2021). On the representation of hyporheic exchange in models for reactive transport in stream and river corridors. Frontiers in
Water, 2, 595538. https://doi.org/10.3389/frwa.2020.595538
Philip, J. (1957). The theory of infiltration: 1. The infiltration equation and its solution. Soil Science, 83(5), 345–358. https://doi.org/10.1097/
00010694‐195705000‐00002
Rahmati, M., Latorre, B., Moret‐Fernández, D., Lassabatere, L., Talebian, N., Miller, D., et al. (2022). On infiltration and infiltration characteristic
times. Water Resources Research, 58(5), e2021WR031600. https://doi.org/10.1029/2021wr031600
Richards, L. A. (1931). Capillary conduction of liquids through porous mediums. Physics, 1(5), 318–333. https://doi.org/10.1063/1.1745010
Roche, K. R., & Dentz, M. (2022). Benthic biolayer structure controls whole‐stream reactive transport. Geophysical Research Letters, 49(5),e2021GL096803. https://doi.org/10.1029/2021GL096803
Rodriguez‐Robles, C. F., Hernandez‐Castro, S., Hernandez‐Escoto, H., Cabrera‐Ruiz, J., Barroso‐Muñoz, F. O., & Terrazas‐Rodriguez, J. E.
(2019). Control analysis of batch reactive distillation column with intermittent fed. Computer Aided Chemical Engineering, 46, 1357–1362.
https://doi.org/10.1016/b978‐0‐12‐818634‐3.50227‐7
Ruth, D., & Ma, H. (1992). On the derivation of the Forchheimer equation by means of the averaging theorem. Transport in Porous Media, 7(3),
255–264. https://doi.org/10.1007/bf01063962
Samorodnitsky, G., & Taqqu, M. S. (1994). Stable non‐Gaussian random processes: Stochastic models with infinite variance. CRC Press.
Sanchez‐Vila, X., Fernàndez‐Garcia, D., & Guadagnini, A. (2010). Interpretation of column experiments of transport of solutes undergoing an
irreversible bimolecular reaction using a continuum approximation. Water Resources Research, 46(12), W12510. https://doi.org/10.1029/
2010WR009539
Scheidegger, A. M., & Sparks, D. L. (1996). A critical assessment of sorption‐desorption mechanisms at the soil mineral/water interface. Soil
Science, 161(12), 813–831. https://doi.org/10.1097/00010694‐199612000‐00002
Şen, Z. (1987). Non‐Darcian flow in fractured rocks with a linear flow pattern. Journal of Hydrology, 92(1–2), 43–57. https://doi.org/10.1016/
0022‐1694(87)90088‐6
Siddiqui, F., Soliman, M. Y., House, W., & Ibragimov, A. (2016). Pre‐Darcy flow revisited under experimental investigation. Journal of
Analytical Science and Technology, 7(1), 1–9. https://doi.org/10.1186/s40543‐015‐0081‐2
Soni, J. P., Islam, N., & Basak, P. (1978). An experimental evaluation of non‐Darcian flow in porous media. Journal of Hydrology, 38(3),231–241. https://doi.org/10.1016/0022‐1694(78)90070‐7
Sun, H. G., Chang, A. L., Zhang, Y., & Chen, W. (2019). A review on variable‐order fractional differential equations: Mathematical foundations,
physical models, numerical methods and applications. Fractional Calculus and Applied Analysis, 22(1), 27–59. https://doi.org/10.1515/fca‐
2019‐0003
Sun, H. G., Meerschaert, M. M., Zhang, Y., Zhu, J. T., & Chen, W. (2013). A fractal Richards’ equation to capture the non‐Boltzmann scaling ofwater transport in unsaturated media. Advances in Water Resources, 52, 292–295. https://doi.org/10.1016/j.advwatres.2012.11.005
Sun, H. G., Zhang, Y., Baleanu, D., Chen, W., & Chen, Y. Q. (2018). A new collection of real world applications of fractional calculus in scienceand engineering. Communications in Nonlinear Science and Numerical Simulation, 64, 213–231. https://doi.org/10.1016/j.cnsns.2018.04.019
Sun, H. G., Zhang, Y., Wei, S., Zhu, J., & Chen, W. (2018). A space fractional constitutive equation model for non‐Newtonian fluid flow.
Communications in Nonlinear Science and Numerical Simulation, 62, 409–417. https://doi.org/10.1016/j.cnsns.2018.02.007
Swartzendruber, D. (1962). Non‐Darcy flow behavior in liquid‐saturated porous media. Journal of Geophysical Research, 67(13), 5205–5213.
https://doi.org/10.1029/jz067i013p05205
Takhanov, D. (2011). Forchheimer model for non‐Darcy flow in porous media and fractures (p. 31). Imperial College.
van Genuchten, M. T. (1980). A closed‐form equation for predicting the hydraulic conductivity of unsaturated soils. Soil Science Society ofAmerica Journal, 44(5), 892–898. https://doi.org/10.2136/sssaj1980.03615995004400050002x
Verma, S., & Viswanathan, P. (2020). Bivariate functions of bounded variation: Fractal dimension and fractional integral. Indagationes Mathematicae, 31(2), 294–309. https://doi.org/10.1016/j.indag.2020.01.006
Wang, F., Liu, Z., Cai, J., & Gao, J. (2018). A fractal model for low‐velocity non‐Darcy flow in tight oil reservoirs considering boundary‐layer
effect. Fractals, 26(05), 1850077. https://doi.org/10.1142/s0218348x18500779

ZHANG ET AL.

13 of 14

19447973, 2024, 3, Downloaded from https://agupubs.onlinelibrary.wiley.com/doi/10.1029/2023WR035782 by Martinus Genuchten - Cochrane Netherlands , Wiley Online Library on [18/03/2024]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License

Water Resources Research10.1029/2023WR035782

Xia, Y., Zhang, Y., Green, C. T., & Fogg, G. E. (2021). Time‐fractional flow equations (t‐FFEs) to upscale transient groundwater flow characterized by temporally non‐Darcian flow due to medium heterogeneity. Water Resources Research, 57(11), e2020WR029554. https://doi.org/
10.1029/2020wr029554
Zhang, Y. (2023). Adjoint models with non‐Fickian reactive transport to identify pollutant sources in water. Journal of Hazardous Materials
Advances, 12, 100331. https://doi.org/10.1016/j.hazadv.2023.100331
Zhang, Y., Baeumer, B., Chen, L., Reeves, D. M., & Sun, H. G. (2017). A fully subordinated linear flow model for hillslope subsurface stormflow.
Water Resources Research, 53(4), 3491–3504. https://doi.org/10.1002/2016wr020192
Zhang, Y., Benson, D. A., & Baeumer, B. (2008). Moment analysis for spatiotemporal fractional dispersion. Water Resources Research, 44(4),W04424. https://doi.org/10.1029/2007wr006291
Zhang, Y., Meerschaert, M. M., Baeumer, B., & LaBolle, L. M. (2015). Modeling mixed retention and early arrivals in multidimensional heterogeneous media using an explicit Lagrangian scheme. Water Resources Research, 51(8), 6311–6337. https://doi.org/10.1002/2015wr016902
Zhang, Y., Sun, H. G., Stowell, H. H., Zayernouri, M., & Hansen, S. E. (2017). A review of applications of fractional calculus in Earth system
dynamics. Chaos, Solitons & Fractals, 102, 29–46. https://doi.org/10.1016/j.chaos.2017.03.051
Zhang, Y., van Genuchten, M. T., Zhou, D. B., Zhang, G. J., & Sun, H. G. (2023). A unified phenomenological model captures water equilibriumand kinetic processes in soil: Codes. December 30, 2023 Release (version 1.0) [Software]. Zenodo. https://doi.org/10.5281/zenodo.10445660
Zhang, Y., Yu, X., Fleckenstein, J. H., Sun, H. G., Lu, C., Yin, M., et al. (2021). Upscaling heat flow in porous media with periodic surface
temperature fluctuation using a one‐dimensional subordinated heat transfer equation. Water Resources Research, 57(7), e2020WR027266.
https://doi.org/10.1029/2020wr027266
Zhao, Y. (2020). Calculating the pore fractal dimension of porous media by using pore size distribution. Data Processing Techniques and Applications for Cyber‐Physical Systems (DPTA 2019), 785–792. https://doi.org/10.1007/978‐981‐15‐1468‐5_92
Zhou, D., Brusseau, M. L., Zhang, Y., Li, S., Wei, W., Sun, H. G., & Zheng, C. M. (2021). Simulating PFAS adsorption kinetics, adsorption
isotherms, and nonideal transport in saturated soil with tempered one‐sided stable density (TOSD) based models. Journal of Hazardous
Materials, 411, 125169. https://doi.org/10.1016/j.jhazmat.2021.125169

ZHANG ET AL.

14 of 14

19447973, 2024, 3, Downloaded from https://agupubs.onlinelibrary.wiley.com/doi/10.1029/2023WR035782 by Martinus Genuchten - Cochrane Netherlands , Wiley Online Library on [18/03/2024]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License

Water Resources Research