Modeling Soil Hydraulic Properties Using Dynamic Variability of Soil Pore Size Distribution

Abstract

The knowledge of temporal variability of soil hydraulic properties (SHPs) in agricultural fields can help in reliable assessment of crop water requirement, thus improving irrigation water usage efficiency. The Fokker–Planck equation (FPE) and its modified forms are popularly used to describe temporal variation in SHPs. These models consider statistical description of soil pore size distribution (PSD) as a probability density function to estimate SHP evolution with time. In this study, we compare four different models to describe the temporal evolution of PSD and SHPs for multiple datasets across the world with different soil types, tillage conditions and crop cover. Further, field experiments were carried out at an experimental agricultural field at IIT Kanpur for rice crops, and the performance of these models was also evaluated for Indian conditions. It is observed that existing models have low accuracy for small pore radii values, and the prediction ability of these models is more affected by soil type rather than tillage conditions. More observations can improve the performance of FPE-based numerical and analytical models. The POWER Model is the least accurate because of its inherent power law assumption of PSD, which results in incorrect values for low pore radii. The FPE analytical model can be reliably used for predicting PSD and SHP evolution at most of the field sites.

1. Introduction

Sustainable management of water resources in an agricultural field requires an accurate assessment of water demand as well as losses. Soil plays an important role in this by not only allowing water and nutrient extraction by roots but also by providing physical support for crop growth [1]. The overall goal of various agricultural management practices, such as tillage, crop rotation, and cropping patterns, is to create optimum characteristics for crop growth and minimize environmental degradation, as well as to optimize water use [2]. Soil hydraulic properties (SHPs) have a significant role in this. SHPs refer to the soil water retention function $\theta \left(h\right)$ and the hydraulic conductivity function $K\left(h\right)$ where $\theta , K$ and $h$ are moisture content, hydraulic conductivity, and suction head, respectively [2,3,4]. The $\theta \left(h\right)$ indicates the ability of the soil to hold and retain water while $K\left(h\right)$ governs the movement of water and solute through the unsaturated zone. SHPs not only govern water movement but also affect plant water and nutrient uptake, thus influencing the dynamics of various hydrological and biological processes within the unsaturated zone [5,6].

SHPs depend upon soil texture, structure, presence of organic matter and chemical and biological properties of the soil. Soil structure-dependent properties are highly dynamic and change with time and depth due to various agronomic practices and environmental factors [7,8]. Various studies have shown that soil physical properties that depend on soil structure, such as moisture content, bulk density, macro-porosity, saturated and near-saturated hydraulic conductivity etc., are heavily influenced by seasonal changes in climate and climate-induced pedological factors [7,9,10]. Furthermore, various researchers have observed that significant changes in soil structure are caused by different soil and crop management practices [11,12,13,14,15,16,17,18]. A common observation from these studies is that the temporal variability of SHPs exceeds the variability induced by crops, tillage, or land use.

These above-mentioned studies have quantified the temporal variation in SHPs, but the findings of these studies are site-specific and vary depending on the climate, soil type, measurement technique and sampling strategy. Furthermore, the SHP measurements in most of the studies have been done at only a few time intervals and have been restricted to the surface soil layer during the study duration.

The SHPs significantly affect various hydrological processes, e.g., evapotranspiration, surface runoff, infiltration, and groundwater recharge. A few studies have looked at the effect of temporal variation in SHP on water movement in agricultural fields. The authors of [19] studied the effect of time-dependent hydraulic conductivity on water flow simulations for an experimental plot located in China under different tillage practices and Maize crop cover. They observed that the incorporation of time-varying hydraulic conductivity improved estimates of percolation losses by 6% under wetting conditions and of water storage by 13–14%. In addition, improvement in estimates of evaporation and transpiration losses was also observed. The authors of [13] investigated the effect of incorporating time-variable SHPs in a soil water simulation and compared the results against measured data for different tillage methods in the near-surface soil profile (0–30 cm). They observed that the inclusion of temporal variability of SHPs led to significant improvement in the performance of the model for all treatments, resulting in average relative errors of less than 13%. Another study [6] attempted to characterize the temporal variation in SHP for an experimental maize plot and analyzed the effect of considering temporally varying SHPs against constant SHPs on simulations of water movement. The incorporation of time-varying SHPs significantly improved soil water predictions up to a 30 cm depth. Another team [20] investigated the temporal and spatial variability of SHPs during a cropping cycle between April and September 2015 for a surface-irrigated maize field in northern Italy and quantified the effect of SHPs on simulations of soil water movement using the FEST- WB model. They observed that a well-parameterized model that accounts for spatiotemporal variation in SHPs improved the accuracy of simulations. Similar observations were made by the authors of [21] in an uncultivated Fluvisol field by using both in situ and laboratory measurements. Recently, the authors of [18] investigated the effect of incorporating seasonal variability of SHPs in hydrological models. They observed the significant impact of the inclusion of short-term temporal variability of SHPs on soil-water dynamics for both ploughed and no-till (NT) condition. Variations up to 44% were observed when the water content was modeled based on SHPs calculated using multiple sampling events vs. that from single sampling events. They postulated that neglecting temporal variability in a hydrological model might lead to an underestimation of the extreme values in a simulation.

The above-mentioned studies highlight the need to incorporate temporal variability in SHPs to improve the simulation of water movement in agricultural fields; however, it is worth noting that SHPs are still considered constant in most of the modeling studies.

In the last two decades, limited efforts have been directed towards investigating the temporal dynamics of SHP [5,6,13,16,17,22,23]. Predictive models for SHPs can be developed using the statistical description of soil pore space as a probability density function (PDF) [24,25,26,27,28]. The authors of [24] proposed a stochastic modeling framework to describe the temporal variability of soil PSD. The Fokker–Planck equation (FPE) was used to combine the probabilistic nature of PSD along with physically based soil deformation models to account for the dynamics of the mean soil pore radius and its variance, as well as total porosity. Using the numerical solution of FPE, they modeled the changes in soil PSD due to the wetting—drying cycles and the resultant Soil Water Retention Curve (SWRC) and $K\left(h\right)$. The initial soil PSD was taken as lognormal.

Another study [25] presented two analytical solutions for modeling PSD evolution as per the modeling framework presented by the authors of [24]. For the first solution, it was assumed that the drift and degradation coefficients depend on time; however, for the second solution, both were dependent on pore size as well as time. Both solutions were applied to model the soil PSD evolution for Millville silt loam soil. The predicted soil PSD was used to obtain $\theta \left(h\right)$, which was found to be in good agreement with experimental data. The initial soil PSD was taken as a lognormal distribution. A further attempt was made by [26] to present an analytical solution for the PSD for the case where drift and degradation coefficients depend on time, and the dispersion coefficients are proportional to the drift coefficients. These coefficients could be estimated from PSD or SWRC data as well as those from mechanistic models.

Following previous work by the authors of [24,25,26], an alternate modeling framework was proposed to describe the temporal evolution of the PSD while considering processes such as tillage, consolidation, and changes in organic matter [28]. Instead of FPE, a modified transport equation was used that neglected the diffusion term and a time-varying power-law PSD was obtained as the solution. The obtained solution is different from earlier studies and has a clear connection to well-known power law for SWRC [29,30]. They derived expressions for soil properties from the model and linked its parameters to soil processes and management activities such as tillage, soil consolidation, and changes in SOM. It was found that the model in this study accurately captured the physical consistency of key soil parameters as they changed with soil biogeochemistry over time and exhibited realistic patterns. It was noted that the proposed model did not account for very short-term changes in soil properties, such as those due to wetting and drying cycles.

The proposed numerical and analytical models attempt to quantify management-induced changes in soil structure and incorporate these results into mathematical functions describing temporal changes in the PSD. They can improve our capacity to assess the overall impacts of different soil and crop management practices on soil physical properties and water balance. However, as of now, the full potential of these models in predicting the temporal variation in SHPs is yet to be explored.

The authors of [31] attempted to characterize the alteration in soil structure due to a crop–pasture rotation at an experimental site at Lincoln University. The measurement results were incorporated into the analytical models proposed in [25,26]. It was observed that the adapted pore size evolution model was suitable for predicting management-induced changes in soil structure. However, a reasonable agreement between the measured and predicted PSD was obtained only when an assumed time-dependent degradation of pores was incorporated into the model. The authors of [32] investigated the effect of crop species with different root systems on field soil pore properties. They conducted a field experiment with twelve species from different families. The parameters of Kosugi’s PSD model were determined inversely from the tension infiltrometer data. The measured root traits were related to pore variables by regression analysis. Then, they used the FPE solution proposed by [25,26] to analyze if observed pore dynamics followed a diffusion-like process. However, they observed that the diffusion-type pore evolution model could only partially capture the observed PSD dynamics.

The authors of [2] reviewed various studies related to the temporal dynamics of SHPs. They applied the analytical solution of FPE found in [25,26] to two SWRC datasets available in the literature to evaluate its suitability to predict the evolution of soil PSD following tillage as well as when there is a change in tillage regime. They observed that the model performed better for change in tillage regimes as compared to the seasonal evolution of soil PSD. Another study [3] further investigated the applicability of the model in [24,25] in terms of the retention data obtained from various studies for different tillage treatments across the globe. The datasets were categorized for two cases: (a) for change in tillage regime and (b) for the temporal dynamics of soil pore space in the months following tillage. They observed that the model provides a reasonable fit. Further, they observed that the applicability of the model was limited due to the lack of adequate datasets for different management practices, soil types, and climate regimes. The availability of such data can be used in calibrating and including new coefficients and sink/source terms for the model.

While the above studies provide useful insight into the performance of the FPE-based models [24,25,26], none of the existing studies have evaluated the performance of all the existing models for a common dataset. Therefore, the prime objective of this study is to perform comparison studies for different models to describe the temporal evolution of soil PSD and, subsequently, the soil hydraulic properties.

2. Governing Equations and Models

In this section, we briefly describe FPE, and the modified transport equation used to describe the time-evolution of the soil PSD.

2.1. Fokker–Planck Equation

The FPE solution describes the evolution of the probability density function (PDF) for the attribute [24]. FPE has been used to model the PDF associated with stochastic processes. Representation of soil PSD in terms of statistical PDF and considering aggregate deformation and porosity loss in wet soils as a solid diffusion process enables modeling soil pore space changes as the evolution of soil PSD is governed by external forcing and soil deformation processes [24].

The FPE is given as follows:

$$\frac{\partial f\left(r,t\right)}{\partial t}=\frac{1}{2}\frac{{\partial }^2}{\partial r^2}\left[D\left(r,t\right)f\left(r,t\right)\right]+\frac{\partial }{\partial r}\left[V\left(r,t\right)f\left(r,t\right)\right]-M\left(r,t\right)f\left(r,t\right)$$

(1)

where $f$ is the PSD, $V\left(r,t\right)$ is the drift coefficient $\left(L{T}^{-1}\right)$ or the infinitesimal mean, $D\left(r,t\right)$ is the diffusion coefficient $\left(L^2{T}^{-1}\right)$ or the infinitesimal variance, and $M\left(r,t\right)$ is the source/sink term or first-order pore degradation factor representing instantaneous pore loss $T^{-1}$. The terms $r$ and $t$ refer to pore radius $\left(L\right)$ and time $\left(T^{-1}\right)$, respectively. The application of FPE to modeling PSD evolution was also motivated by its similarity with the advection-dispersion equation and a similar decay term that may allow the use of existing analytical solutions. The boundary conditions of the above equation used are.

$$f\left(r,0\right)=f_0\left(r\right), 0<r<\infty$$

(2)

$$Vf-D\frac{\partial f}{\partial r}=0, r=0,t>0$$

(3)

$$\frac{\partial f}{\partial r}=0, r\to \infty ,t>0$$

(4)

with $f_0$ as the initial PSD. The lower boundary condition stipulates a zero-probability flux, i.e., pores cannot assume a negative radius. To ensure the conservation, the upper boundary condition requires a zero gradient, i.e., a zero-probability flux for infinitely large pores. The homogeneous zero flux conditions imply that any loss of probability (PSD) is due to degradation. The coefficients $V\left(r,t\right), D\left(r,t\right)$ and $M\left(r,t\right)$ can be obtained as defined in [32].

The mean change in pore radius and the variance of changes in pore radii with time are represented by $V\left(r,t\right)$ and $D\left(r,t\right),$ respectively. The moments are defined by integrating the PSD with respect to the pore radius.

$$m_n={{\displaystyle \int }}_0^{\infty }{r}^{n}f\left(r,t\right)dr$$

(5)

For example, a zero-order moment $\left(m_0\right)$ is given as:

$$m_0={{\displaystyle \int }}_0^{\infty }f\left(r,t\right)dr$$

(6)

Normalized moments $\left(M_n\right)$ are calculated through the division of $m_n$ by $m_0$. The first-order normalized moment $M_1$ characterizes the mean pore size $\langle r\rangle$ and the second-order centralized moment ${\mu }_2$ characterizes the variance $\left({\sigma }^2\right),$ and they can be obtained as:

$$M_1=\langle r\rangle =\frac{m_1}{m_0}=r_m\exp \left(\frac{{\sigma }^2}{2}\right)$$

(7)

$${\mu }_2=r_m^2\exp \left({\sigma }^2\right)\left[{\sigma }^2-1\right]$$

(8)

$V$ can be obtained by:

$$V\left(t\right)=\frac{d}{dt}\langle r\rangle =a\left(1-\frac{\langle r\rangle }{b}\right)\langle r\rangle ;\mathrm{where}\langle r\rangle =\frac{b\langle r_0\rangle }{\langle r_0\rangle +\left(b-\langle r_0\rangle \right)\exp \left(-at\right)}$$

(9)

The cumulative drift term (T) is given by:

$$T=\langle r\rangle -\langle r_0\rangle ={{\displaystyle \int }}_0^{t}V\left(\tau \right)d\tau$$

(10)

$D$ and $V$ are related as:

$$\lambda =\frac{D\left(t\right)}{\left|V\left(t\right)\right|}$$

(11)

2.1.1. Analytical Solution

The analytical solution of FPE (Equation (1)) subjected to the initial and boundary conditions as per Equations (2) and (3) yields the following solution as found in [25,26]:

$$\begin{array}{ll}f\left(r,t\right)=& \exp \left({{\displaystyle \int }}_0^T\frac{M\left(\tau \right)}{V\left(\tau \right)}\right){{\displaystyle \int }}_0^{\infty }{f}_0\left(\xi \right)\\ & \times \{\frac{1}{\sqrt{4\pi \lambda T}}\left[\exp \left(-\frac{{\left(r-\xi +T\right)}^2}{4\lambda T}]+\_\right)+\exp \left(-\frac{r}{\lambda }-\frac{{\left(r+\xi -T\right)}^2}{4\lambda T}\right)\right]\\ & +\frac{1}{2\lambda }exp\left(-\frac{r}{\lambda }\right)erfc\left(\frac{\left(r+\xi -T\right)}{\sqrt{4\lambda T}}\right)\}d\xi \end{array}$$

(12)

where $\tau$ and $\xi$ are the dummy integration variables. As the tillage treatment cannot be easily transferred to $t$ as an independent variable, $T$ is used in its place, which means that the evolution of PSD is projected based on the gradual changes in $r$.

$M$ was taken zero for modeling the soil PSD evolution [25,26]. The same was considered in [33] while investigating the effect of crop species with different root systems on field soil pore properties. However, authors of [31] attempted to characterize $M$ by describing it as an exponential term which was then used by another study [2]. The same is used in this study as well.

2.1.2. Numerical Solution

The numerical solution of FPE has been used by [24] to describe the temporal variability of SHPs. In this study, the FPE equation was discretized by using the Crank–Nicolson scheme in time and centered finite differences in space to write the equation in the form of AX = B. It was assumed that the coefficients $V$ and $D$ vary with time only and degradation,$M=0$. Further, the LM algorithm was used to determine the value of $\mathsf{\lambda }$.

To verify our solution, both the analytical and numerical solutions of FPE were used to model the PSD evolution for a time-dependent $V$ and $M=0$ as per test problem (3) (see Figure A1). It was observed that except at extremely low $r$, the numerical model was able to simulate and match the analytical solution of FPE. This might be due to the distribution not truly being lognormal and the analytical solution becomes less accurate when $r=0$ [25].

2.2. Transport Equation

The authors of [28] proposed a time varying power law PSD as the solution of the generic form of the transport equation as given below to model PSD evolution.

$$\frac{\partial f\left(r,t\right)}{\partial t}=\frac{\partial }{\partial r}\left[V\left(r,t\right)f\left(r,t\right)\right]-M\left(r,t\right)f\left(r,t\right)$$

(13)

The drift term $V\left(r,t\right)$, and source/sink term $M\left(r,t\right)$ of the evolution equations are defined as:

$$V\left(r,t\right)=\frac{r}{a\left(t\right)b\left(t\right)}\left(a\left(t\right)b^{\prime \left(t\right)}ln\left(r\right)-a^{\prime }\left(t\right)\right)$$

(14)

$$M\left(r,t\right)=\frac{b^{\prime }\left(t\right)}{b\left(t\right)}\left(1+\ln \left(r\right)\right)-\frac{a^{\prime }\left(t\right)}{a\left(t\right)b\left(t\right)}$$

(15)

where the time-derivative of the parameters is indicated with a prime (′) symbol. This was done to modify the given form of evolution equation which then was solved using the method of characteristics to obtain a power law solution. The initial condition for this problem is:

$$f\left(r,0\right)=f_0\left(r\right)=a\left(0\right)r^{-b\left(0\right)}, 0\le r\le R_m\left(0\right); t=0$$

(16)

The time-varying power law PSD was obtained as

$$f\left(r,t\right)=a\left(t\right)r^{-b\left(t\right)}, 0\le r\le R_m\left(t\right)$$

(17)

where, $a\left(t\right)$ is a scaling parameter, $b\left(t\right)$ is the power law exponent, and $R_m\left(t\right)$ is the maximum effective pore radius. In addition, $R_m$ is related to the air-entry or bubbling pressure, which can be defined as the matric pressure at which air enters the soil pores. Porosity $(\varnothing )$ was obtained by the integral of the PSD over $r$.

$$\varnothing \left(t\right)=\frac{a\left(t\right)R_m{\left(t\right)}^{1-b\left(t\right)}}{1-b\left(t\right)}$$

(18)

Obtaining the expression of the scaling parameter $a\left(t\right)$ and putting in $f\left(r,t\right)$, the resultant expression was:

$$\begin{gathered} f\left(r,t\right)=\left(\frac{\varnothing \left(t\right)\left(1-b\left(t\right)\right)}{R_m{\left(t\right)}^{1-b\left(t\right)}}\right)r^{-b\left(t\right)} \\ b\left(t\right)<1 \end{gathered}$$

(19)

The authors of [28] proposed that the power law parameters can be modified and expanded to account for the various processes in the soil. The power law description of the PSD is similar to the power law model and has a basis in the fractal fragmentation of soils [29,30].

2.3. Modeling SWRC Using Soil PSD

The evolution of PSD is reflected in SWRC. The relationship between soil pore radius $\left(r\right)$ and suction head $\left(h\right)$ is given as:

$$r=\frac{A}{h}, \mathrm{where} A=-0.149 {\mathrm{cm}}^2 \mathrm{if} h \mathrm{is} \mathrm{in} \mathrm{cm}$$

(20)

The slope of the SWRC, water capacity function $C\left(h\right)$ is given as (21)

$$C\left(h\right)=\frac{d\mathsf{\theta }}{dh}=f\left(r\right)\frac{dr}{dh}$$

(21)

The soil moisture $\theta \left(h\right)$ can be obtained as:

$$\begin{array}{ll}\theta \left(h\right)=& {{\displaystyle \int }}_{-\infty }^{h}C\left(h\right)+{\theta }_r={{\displaystyle \int }}_0^{r}f\left(r\right)dr+{\theta }_r\\ & \mathrm{where} {\theta }_r=\underset{h\to \infty }{\lim }\theta \left(h\right)=\underset{r\to 0}{\lim }f\left(r\right)\end{array}$$

(22)

In this study, the initial soil PSD was assumed to be lognormal. Kosugi’s parameters [34] can be obtained from the van Genuchten (VG) parameters [35] using the following equations.

$$m=1-\frac{1}{n}$$

(23)

$$h_0=\frac{m^{1-m}}{\alpha }$$

(24)

$${\sigma }^2=\left(1-m\right)\ln \frac{2^{\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$m$}\right.}-1}{m}$$

(25)

$$h_m=h_0\exp \left({\sigma }^2\right)$$

(26)

For the POWER model [27], the resultant SWRC curve is given as

$$h=-\left(\frac{A}{R_m\left(t\right)}\right){\left(\frac{\theta \left(t\right)}{\varnothing \left(t\right)}\right)}^{-1/\left(1-b\left(t\right)\right)}$$

(27)

3. Methodology and Data

3.1. Methodology

To meet the objective of this study, the following methodology will be followed: (a) formulation of different models for PSD evolution; (b) estimation of Kosugi’s parameters from initial SWRC datasets; (c) determination of coefficients for FPE for analytical and numerical models as well as those of power model; (d) use of the LM algorithm to perform curve-fitting to observed final PSD data to estimate $\lambda$ for the FPE_1, FPE_2, and FPE_NM models and to estimate $b$ and $R_m$ for the power model; and (e) a comparison of different models over different datasets.

The steps used in the methodology have been described in Figure 1. The different models used in this study are as follows:

(a)

FPE_1: It represents the analytical solution of FPE found in [24,25]. In this model, $M$ was neglected, due to which the first part of the analytical solution, i.e., the exponential term (refer to Equation (12)), was taken as 1. To obtain $\lambda$, the Levenberg–Marquardt (LM) method was used to fit the observed PSD to the one given by Equation (12).

(b)

FPE_2: It represents the analytical solution of FPE used by [3,30]. In this, $M$ was described as $M=d\times \exp \left(ct\right)$. $\lambda , c$ and $d$ were obtained using the LM method, same as that of the FPE_1 model.

(c)

NM: It is the numerical model obtained from the finite difference-based discretization of the FPE equation used by [24]. In this model, $M$ was neglected.

(d)

POWER: It represents the power law model given by [28]. The model coefficients were obtained by curve fitting using the LM method.

3.2. Data

In this study, to compare the performance of four different models to predict the temporal evolution of PSD and SWRC, datasets from various studies across the world have been collected. In addition, the comparison was performed on the SWRC data collected from an experimental agricultural plot at IIT Kanpur, India, for the rice season in 2022. The WP4C Dewpoint Potentiometer from the Meters’ group was used to estimate SWRC from the soil samples taken from the experimental plot. Figure A2 shows a world map showing the various datasets used in the study as well as the experimental site at IIT Kanpur.

Table 1. Summary of different soil datasets and their classification upon soil type, agriculture practice, cover, crop grown, etc., used in the study.

S. No.	Study Area	Study Period	Soil Type	Field and Tillage Condition	Data Collection and Summary	Reference
1	Germany	1984	Silt loam	No tillage and annual moldboard tillage	The study was conducted in Germany over a period of 18 years. SWRC were determined for the pressure head range (h = −20, −50, −100, −300, −1000, −15,000 cm) on undisturbed soil samples.	[36]
2	West-central Alberta	2001	Inseptic and oxyaquic, haplocryalfs, and typic dystrocryepts	Empty and loaded pass	The study was conducted on boreal forest soils over different traffic cycles for fourteen sites. At every site, skidding was performed for 3, 5 and 7 cyclic passes, and after each pass, soil samples were collected from 4 randomly selected sampling points. SWRC were determined for the pressure head range (h = −20, −50, −100, −300, −1000, −15,000 cm) on undisturbed soil samples.	[37]
3	Washington, USA	2001–2002	Chernozem	CT with winter wheat (Triticum aestivum L.) and spring pea (Pisum sativum) and NT with winter wheat	Primary tillage operations were performed in October 2000 and 2001 and soil samples were taken on 16 May (46 days after tillage) and 30 November in 2001 and on 22 April, 26 June and 18 September, 18 in 2002. Undisturbed soil sample were used to determine the SWRC at different pressure heads (h = −1, −2, −4, −8, −10.5, −12.5, −30, −40, −100, −200, −300, −400, −600, −1000 cm).	[38]
4	University of Chile, Antumapu, Chile	2003	Sandy clay alluvial soil	CT and NT with spring wheat (Triticum turgidum L. var. durum)	Two sites were prepared, one for 4 years and another for 7 years. The soil samples were collected between July 2003 and January 2004 at depths 0–2, 2–5, and 5–15 cm at two sites per plot. Undisturbed soil cores were used to obtain SWRC at different pressure heads (h = −100, −330, −1000, −3000 and −15,000 cm) by fitting the data to the VG model.	[39]
5	Raasdorf, Austria	2008–2010	Silt loam	CT, RT and NT with winter wheat (Triticum aestivum)	Ten infiltration measurements were made using tension disc infiltrometer between August 2008 and June 2010 at h = 0, −1, −4, −10 cm for the topsoil layer (0–30 cm) to inversely estimate the parameters of the VG model.	[14]
6	Canterbury, New Zealand	2008–2009	Haplic Cambisol	CT with wheat (Triticum aestivum L.) and pasture (G) with sheep and beef cattle grazing	Undisturbed soil cores were used to determine SWRC at different pressure heads (h = −10, −31, −100, −316, −1000, −15,800 cm) and fitted to the Kosugi model.	[31]
7	Raasdorf, Austria	2009–2012	Chernozem	CT with fallow and with mustard (Sinapis alba L.) and rye (Secale cereal L.)	Three different soil cover treatments were integrated during the fallow period between cash crops of a three-year crop rotation. Infiltration experiments were conducted twelve times between September 2009 and July 2012. Tension disc infiltrometer method was used (h = 0, −1, −4, −10 cm) for the topsoil layer (0–30 cm) to inversely estimate the parameters of the VG model using the 2D Richards’ equation.	[7]
8	Lamothe, France	2012	Stagnic Luvisol with an illuvial clay horizon	Spring moldboard plowing with maize	Undisturbed soil samples were collected at sowing, maize flowering and one month before harvest. Tension disk infiltrometers were used between pressure heads −15 and −1 cm at the soil surface and at 20 cm depth. SWRC was made using measurements at different pressure heads (h = 0, −16, −33, −100, −330, −1000, −6300, −10,000 and −16,000 cm) using pressure plates.	[6]
9	Zhongghou Catchment, China	2012–2013	Silt loam	CT at the beginning. The end area was covered by black locust (Robinia pseudoacacia L.) plantation.	The field measurements were conducted under black locust plantation, grass, and rape cultivation in June and September 2012 and March 2013. For the cropland, the measurements were carried out 2 weeks after harvesting, 2 weeks after sowing, and during the flowering phase of rape. Undisturbed soil cores were used to determine the SWRC by means of evaporation method (h = −10, −20, −31, −63, −100, −200, −316, and −15,000 cm).	[40]
10	Attica, Greece	2011–2013	Eutric Fluvisol	RoT	SWRC were determined in the laboratory using a sand-kaolin box for the pressure head range h = 0 to −200 cm (h = 0, −40, −80, −120, −160, −200 cm) on undisturbed soil samples.	[15]
11	Zaragoza, Spain	November 2011– November 2012	Hypercalcic Calcisol	CT and NT in the long fallow phase	Five different soil samplings during the fallow period were performed (a) pre-tillage (b) post-tillage before rain (c) post-tillage and rain; (d) late fallow; and (e) end-fallow. The undisturbed soil samples were saturated, and decreasing pressure heads were sequentially applied (h = −5, −15, −30, −100, −500, −1000, −5000 and −15,000 cm) to obtain respective SWRC.	[41]
12	North Dakota State University, western North Dakota, USA	2014–2019	Lihen sandy loam (sandy, mixed, frigid Entic Haplustoll).	CT and NT with crop rotation of corn (Zea mays L.)-soybean (Glycine max [L.])	Undisturbed soil cores were collected from 0–15 cm and 15–30 cm depths in corn rows at one sample per plot on 8 October 2014; 26 October 2015; October 2016; and 16 October 2017. Soil core sampling was replicated five times. Measured SWRC obtained by the HYPROP.	[42]
13	Indian Institute of Technology Kanpur, Uttar Pradesh, India	September–December 2022.	Silt loam	NT with wheat	This study was conducted on an experimental plot at the Indian Institute of Technology Kanpur, Uttar Pradesh, India from September–December 2022. Soil samples were taken at depths of 10 cm, 25 cm and 50 cm at monthly intervals. Saturated hydraulic conductivity was measured using a falling head permeability test, gravimetric soil moisture content (θ) was determined from the oven drying method, and suction head (ψ) was obtained using WP4C Dew Point Potential Meter (manufactured by METER Group, Inc. Pullman, WA, USA).

provides the summary of different soil datasets and their classification based on soil type, agriculture practice, land use/cover, and crop type, which have been used in this study for comparison of different models. From the datasets under consideration, the tillage regimes considered are conventional tillage (CT), usually using a moldboard plow, no tillage (NT) with none or very minimal disturbance to the soil, and rototillage (RoT).

4. Results and Discussion

The ability of the four different models described earlier to predict PSD evolution, and SWRC was evaluated in this section. The SWRC was obtained for both initial and final conditions using the VG equation for the observed data. The suitability of the models was checked for various agricultural and management practices, such as the crop used, tillage conditions, effects of land-use change and soil type.

Figure 2 shows the prediction of PSD and SWRC evolution over time by different models for chernozem soil in Washington, USA [38]. All the models predicted the PSD evolution very well within the range of observed values. The POWER model over-predicts the PSD values for both the initial (CT) as well as final condition (NT) when $r<1$ μm (see Figure 2a,c). FPE_1, FPE_2 and NM perform well when predicting PSD evolution. However, for predicting SWRC, it was observed that the FPE_1, FPE_2 and NM models underpredict compared to the observed values, while the POWER model is only able to partially predict the values (see Figure 2b,d). The increase in performance of the POWER model can be attributed to higher values of observed $r$.

Figure 3 shows the comparison of performances of different models for predicting the evolution of both soil PSD and SWRC against the observed values for Haplic Cambisol soil from the New Zealand dataset [31]. It was observed that all the models performed well except for the NM model, which underpredicted the soil PSD when $r<1 \mathsf{\mu }\mathrm{m}$ (see Figure 3a,c). Similarly, for predicting SWRC, the FPE_1 and FPE_2 models performed well despite FPE_1 slightly underpredicting. Both the NM and POWER models slightly underpredicted and overpredicted the final SWRC values, respectively (see Figure 3b,d).

Figure 4 shows the performances of different models for predicting PSD and SWRC evolution from CT to the plantation for Eutric Fluvisol soil in China [40]. It was observed that NM performed best compared to other models for predicting soil PSD. The FPE_1 and FPE_2 models predicted similar values for both soil PSD as well as SWRC evolution. The PSD values predicted by FPE_1 and FPE_2 was higher as compared to the observed PSD values, while POWER overpredicted soil PSD (see Figure 4c). For predicting SWRC, POWER performed well as compared to others. The FPE_1 and FPE_2 models overpredicted SWRC near saturation and underpredicted it near dryness. NM predicted the least SWRC values for both initial as well as the final condition (see Figure 4a–d).

Figure 5 shows the comparison of different models for predicting PSD and SWRC evolution for Eutric Fluvisol soil from Greece [18]. It was observed that within the observed range of pore radii, all the models accurately predicted PSD as well SWRC (see Figure 5a–d); however, when $r<1 \mathsf{\mu }\mathrm{m}$, the POWER model highly overpredicts soil PSD as compared to other models for both initial as well as final condition. For SWRC, when $r<1 \mathsf{\mu }\mathrm{m}$, POWER overpredicted the values for the initial state, while for the final condition, the model performed well like other models.

Figure 6 shows the performance of different models for predicting PSD and SWRC evolution from CT to NT for Hypercalcic Calcisol soil in Spain [41]. The performance of the NM model was best while predicting both soil PSD as well as SWRC evolution. The lowest performance was observed for the POWER model (see Figure 6a–d). FPE_1 and FPE_2 predicted similar values, but they underestimated both soil PSD and SWRC values for the final condition (see Figure 6a–d).

The performances of different models were compared for predicting PSD and SWRC evolution from CT to NT for Lihen sandy loam soil from the Dakota, USA dataset [42] (Figure 7). It was observed that the FPE_1, FPE_2 and NM models performed well for predicting soil PSD evolution, whereas the POWER overpredicted soil PSD when $r<10 \mathsf{\mu }\mathrm{m}$ and underpredicted soil PSD values when $r>10 \mathsf{\mu }\mathrm{m}$ for both initial as well as final condition (see Figure 7a,c). For predicting SWRC evolution, the POWER model performed well as compared to the others. FPE_1 and FPE_2 predicted similar results, which were lower than observed SWRC values for both initial as well as final conditions, while NM highly underpredicted SWRC values (see Figure 7b,d).

Figure 8 shows the performance of different models for predicting PSD and SWRC evolution from t = 1 day to t = 30 days under NT conditions for silty loam soil in the experimental agricultural plot at IIT Kanpur, India. Both the FPE_1 and FPE_2 models predicted similar results. However, we were unable to fit the NM and POWER into the observed data. While the performance of the models was poor, they were able to show the temporal evolution of both PSD and SWRC. Modeling of seasonal, temporal variability of SHPs as well as the inclusion of SHPs information for different subplots, may lead to improvement in the performance of the models.

The performance of the models was compared for several more datasets, the results of which are not shown for brevity (see Appendix A). Different models for predicting PSD evolution and SWRC were compared for the German dataset of silty loamy soil [34] (see Figure A3). It was observed that the FPE_1 and FPE_2 models well predicted PSD and SWRC evolution. The PSD and SWRC curves by NM and POWER models were found to be underpredicted and overpredicted, respectively. NM predicted PSD and SWRC with higher accuracy when $r>1 \mathsf{\mu }\mathrm{m}$ and the model performance decreases when $r<1 \mathsf{\mu }\mathrm{m}$. POWER model underpredicted the initial soil PSD and SWRC at initial conditions when $r>1 \mathsf{\mu }\mathrm{m}$ and overpredicted the final PSD and SWRC when $r>1 \mathsf{\mu }\mathrm{m}$.

The four models were compared for datasets from Boreal Forest soils due to compaction by harvesting equipment [35] for the 0th and 3rd skidding cycles (see Figure A4). It was observed that both FPE_1, FPE_2 and POWER models underpredict the final PSD while the final SWRC predicted by FPE_1 and FPE_2 was quite close. The performance of the POWER model was extremely poor for this dataset. The prediction by NM was close to the final soil PSD as well as that of SWRC values (Figure A4b,d).

Figure A5 shows the performance of different models in predicting PSD and SWRC evolution for the 4-year dataset in Chile at depths 0–2 cm and 2–5 cm for sandy clay alluvial soil for initial (CT) and final conditions (NT) [39]. It was observed that the soil PSD as well as SWRC values predicted by FPE_1 and FPE_2 were quite similar and close to the observed values of PSD and SWRC for both depths 0–2 cm and 2–5 cm. FPE_1 and FPE_2 had the best performance among all the models except for predicting SWRC at a depth 0–2 cm, where they overpredicted. For soil PSD prediction, it was observed that the POWER model overpredicts the values for low pore radii for depths 0–2 cm and 2–5 cm. NM overpredicts the values at depths of 0–2 cm, and between 2–5 cm depth, the accuracy of the NM model was improved. However, for predicting SWRC values, it was observed that at the initial condition, the POWER model over and underpredicted at depths 0–2 cm and 2–5 cm, respectively; however, at the final condition, the values predicted by POWER were close to observed values (see Figure A5g,h). For predicting SWRC, NM performed best as compared to others (see Figure A5c,d,g,h). It was observed that with an increase in depth from 0–2 cm to 2–5 cm, the performance of the models improved, and they were able to model soil PSD evolution well.

The PSD and SWRC evolution for the Chile 7-year dataset [39] at depths 0–2 cm and 2–5 cm for sandy clay alluvial soil for initial (CT) and final condition (NT) were best described by the NM (see Figure A6a–e). Both NM and POWER models predicted the final SWRC better as compared to the FPE_1 and FPE_2 models at both 0–2 cm and 2–5 cm depths (see Figure A6g,h). Furthermore, it was observed that, like in the previous case, an improvement in the performance of the FPE_1 and FPE_2 models was observed for the prediction of soil PSD upon an increase in depth from 0–2 cm to 2–5 cm. In addition, the POWER model was found to overpredict the SWRC as well as soil PSD values for the initial condition, but its performance improved for the final condition as well as with an increase in depth.

For the silt loam soil dataset from Austria [14], soil PSD and SWRC evolution were modeled for initial (CT) and final conditions (NT). It was observed that the performances of all the models were satisfactory when $r$ lies between the range of observed $r$ values (see Figure A7a,b). However, when $r<100 \mathsf{\mu }\mathrm{m}$, the performance of the models varies greatly. The POWER model predicted very large PSD values, while the FPE_1, FPE_2, and NM models predicted lower PSD values in correspondingly decreasing order (see Figure A7a,b). All the models performed well in terms of predicting SWRC. POWER model captured the final SWRC despite overpredicting the initial SWRC. In addition, for the final condition, the SWRC curves from FPE_1 and NM models are higher than that of observed SWRC values, while the FPE_2 model predicted SWRC values close to the observed values (see Figure A7b,d).

Different models were compared in Figure A8 for predicting soil PSD and SWRC evolution against observed values for chernozem soil in the Austria-II dataset [7]. It was observed that all the models performed well while predicting PSD evolution from initial to final state for the observed range of $r$ values; however, when predicting SWRC, the performance of the models varies greatly. When $r<100 \mathsf{\mu }\mathrm{m}$, both FPE_2 and POWER models predict high values of the soil PSD while FPE_1 predicts the low PSD values, the predicted PSD values by NM lie in between (see Figure A8b). While predicting SWRC, the POWER and NM models were the most accurate and FPE_2 highly overpredicted, and FPE_1 underpredicted the SWRC values (see Figure A8d). Furthermore, the POWER predicted high values of SWRC as well as soil PSD for the initial condition (see Figure A8a,b), which might be due to its inherent power law assumption.

Different models were compared for Stagnic Luvisol soil with Spring moldboard plowing with maize in France [6]. It was observed that within the observed range of pore radii, POWER and NM models exhibit the best performance in predicting PSD as well SWRC, whereas FPE_1 and FPE_2 models overpredicted (see Figure A9a–d); however, when $r<1 \mathsf{\mu }\mathrm{m}$, it was observed that all the models highly overpredict PSD and SWRC.

4.1. Effect of Tillage

In the above section, it was observed that the performance of different models varies across different sites across the world. In this section, we attempt to discuss how the performance of models varies for similar tillage practices across the world.

For the USA-I (Figure 2), Chile (Figure A4 and Figure A5), Austria-I (Figure A6), Spain (Figure 6) and USA-II (Figure 7), the agricultural practice changes from CT to NT [37,38,40,41]. It was observed that NM performed well in most of the cases while the performance of FPE_1 and FPE_2 was good, but they overpredicted the SWRC values, whereas the POWER performed poorly in most cases except for the Chile dataset, where it was able to predict PSD and SWRC values close to that observed at final conditions despite predicting lower values at initial conditions. Overall, it can be postulated that POWER performed poorly.

The effect of the change of tillage practice from CT to pasture; Moldboard ploughing to NT; CT to the forest were observed in Germany (Figure A3), New Zealand (Figure 3), and China (Figure 4), respectively. It was observed that the FPE_1 and FPE_2 models did not perform as expected; however, the NM and POWER models performed well in those cases.

Overall, it is observed that the FPE_1, FPE_2 and NM Models were able to predict PSD and SWRC evolution for change in tillage as well as the continuation of the same tillage. For better performance of models, higher numbers of observations are required. NM model performed poorly for low pore radii but better as compared to the POWER model.

4.2. Effect of Soil Type

The effect of soil types on models’ performance efficiency was observed for silt loam in Austria-I (Figure A8), Germany (Figure A3) and China (Figure 4) [13,36,40], for sandy clay alluvial soil in Chile (Figure A5 and Figure A6); for Chernozem soil in USA (Figure A2) and Austria-II (Figure A9); Hypercalcic Calcisol soil in New Zealand (Figure 3), for Eutric Fluvisol soil in Greece (Figure 5), respectively. It was observed that for silt loam soil, the values of soil PSD and SWRC for FPE_1 and FPE_2 are quite close to each other as well as close to the observed values except for the China dataset. The POWER and NM models performed well for predicting both PSD and SWRC values.

For sandy clay alluvial soil in Chile, 4- and 7-year datasets from two stations were used to evaluate the models’ performance. It was observed that for silt loam soil, FPE_1 and FPE_2 predicted similar values of PSD and SWRC, which were close to the observed values. NM had the next best performance for both sets of datasets, and while the POWER model underpredicts PSD and SWRC values at the initial state, however, its performance improves for the final condition. Further, it was observed that with an increase in depth, the models’ performance improved (see Figure A5 and Figure A6). For Chernozem soil, it was observed that the FPE_1, FPE_2 and NM models predicted close values of soil PSD but underpredicted SWRC; the POWER model performed best in predicting both PSD and SWRC (see Figure 2 and Figure A9). Furthermore, it was observed that FPE_1 and FPE_2 models predicted different values for both soil PSD and SWRC, unlike previous cases.

It is observed that the performance of the models is highly dependent upon the soil type as it determines the soil PSD as well as available observation data. FPE_1 and FPE_2 predicted similar values and were able to predict the soil PSD and SWRC well for most soil types. NM performed similarly to the FPE_1 and FPE_2 models when predicting both SWRC and PSD. However, the performance of POWER varied highly from the initial condition to the final condition depending upon soil type. POWER performed well if a higher number of observations were present. It was also able to improve its prediction of the final condition despite poor estimation of initial PSD and SWRC values.

5. Conclusions

SHPs exhibit significant spatiotemporal variability and have a strong influence on flow in unsaturated soils. SHPs depend upon the soil structure, which can be described using the soil PSD. The soil PSD and its statistical description as a probability density function are used in the existing predictive models to describe the temporal variation in SHPs. In this study, we compared the performance of four models to describe the temporal evolution of soil PSD and SHPs for various datasets across the world. In addition, the model’s performance was compared to SWRC data collected from an experimental agricultural field at IIT Kanpur. These models were compared for different tillage conditions as well as changes in tillage conditions, soil types, and crop types. Further, an attempt was made to compare the performance of the models for both short-term as well as long-term temporal variability.

Among the four models, two were obtained from the analytical solution of the FPE equation (FPE_1 and FPE_2), and the other two from the POWER and NM models were obtained from the numerical solution of the FPE and power law solution of the generic transport equation. It was observed that when pore radii values were low, the accuracy of the models decreased. The models perform better when the soil pore radii are comparatively high, especially for the POWER and NM models. This correlates to the earlier assumption by [25] that pore-size evolution may affect larger pores as compared to smaller pores. Further, the NM model did not perform well when the number of observed matric potential values were low. The predictive ability of different models was not much affected by changes in tillage condition as that by soil type. The performance of the models was not affected by the crop type. The NM model’s performance can be improved by having more observations, whereas the POWER model is not recommended due to its inherent assumption that soil PSD should follow power law which leads to erroneous values for low soil pore radii. The POWER model performed well when a low number of observations were present, while the FPE_1 and FPE_2 models performed well for most of the cases. For the SWRC dataset obtained from an experimental plot IIT Kanpur, the performance of the POWER and NM models was poor, while the FPE_1 and FPE_2 models performed similarly to each other.

It was observed that models based on the analytical solutions of FPE (FPE_1 and FPE_2) were the best among the models considered, and between them, FPE_1 should be preferred to model the temporal variability of SHPs and PSD. Furthermore, to improve the performance of models, more observations are required for modeling both short and long-term temporal variation in SHPs. Future studies should focus not only on modeling the temporal variability of SHPs but also on relating the parameters of FPE with the various agricultural practices in the field, which can be useful for effective water management.