SPorDyn: A Python code for modeling the evolution of soil pore size distribution after tillage

Graphical abstract


Method details
Surface soil structure is responsive to both natural and anthropogenic changes and consequently impacts soil hydraulic properties (SHP). Generally, the structural pore space with pore radius (r) ! 5 mm is expected to be the most affected by management practices such as tillage and natural stresses (e.g. rainfall). The need for inclusion of soil pore dynamics in hydrological models to improve their reliability and accuracy has been stressed in recent times [1][2][3]. In this context, Or et al. [4] proposed to use a partial differential equation (PDE) known as the Fokker-Planck Equation to capture the dynamics of soil pore size distribution (PSD) following tillage with respect to time and pore radius. Knowledge of the PSD will pave way to predict unsaturated SHP which can then be incorporated in hydrological models. The PDE comprises physically-based coefficients (drift (V), degradation (M) and dispersion (D) coefficients) and they encompass our perception of the mathematical behavior of soil PSD in response to tillage practices. The coefficients are subject to an initial condition as well as upper and lower boundary conditions. Based on these conditions, an analytical solution to the proposed PDE was provided by Leij et al. [5,6]. In Chandrasekhar et al. [1], the model was applied to different water retention parameter data sets around the world to evaluate its capability in predicting the temporal dynamics of soil pore space for two cases (1) when there is a change in the tillage regime/land-use and (2) in the months following tillage. In the present contribution, we share the Python code used in Chandrasekhar et al. [1] to capture the evolution of soil pore space following tillage. This paper is organized as follows: As a first step, a skeletal framework of the mathematical model is briefly described. For a detailed overview of the model and its application, the reader is directed to Chandrasekhar et al. [1]. The required input data is listed. Brief descriptions of the remaining aspects of the model such as the coefficients as well as the accompanying code are then provided. Finally, the optimization process and steps to obtain the analytical solution are given.

Mathematical model
In this section, a skeletal overview of the mathematical model and its coefficients is provided. The following partial differential equation (PDE), also known as the Fokker-Planck equation, is used to describe the evolution of soil PSD following tillage [4]: where f is the PSD or frequency [L À1 ] of pores as a function of time t [T] and pore radius r and V quantify the changes with time of the variance of the PSD and mean pore radius, respectively, while M is a first-order degradation factor representing instantaneous pore loss, i.e., the fraction of pores that are lost due to instantaneous collapse. The initial and boundary conditions for solving Eq. (1) are: Vf À D @f @r f 0 in Eq. (2) is the initial PSD for the parameterization of which the lognormal distribution function of [7] is used. The lower boundary condition (Eq. 3) stipulates a zero-probability flux meaning that only positive pore sizes are allowed while the upper boundary condition (Eq. 4) necessitates a zero gradient for infinitely large pore radii. The PDE subject to the initial and boundary conditions yields the following analytical solution [5,6]: where t and j are dummy integration variables. Due to the fact that tillage treatment cannot be readily converted to time as an independent variable, the cumulative drift term T is used as an independent variable instead of time, meaning that the evolution of PSD is predicted based on the gradual changes in the pore radii [5].

Input data
The input data for the modeling approach are water retention parameters obtained at different temporal stages. These parameters may be parameterized according to the Kosugi [7] water retention model. If Van Genuchten parameters [8] are available, they can be converted to the Kosugi parameters which is explained in Chandrasekhar et al. [1]. The required input water retention parameters are listed in Table 1.

Initial pore size distribution
Eq. (6) shows the lognormal distribution function for the pore radii [7]: where is the total initial porosity, r m [L] is the initial median pore radius or geometric mean and s [-] is the standard deviation of the log-transformed pore radius. Leij et al. [5] assumed f 0 to be the difference between saturated and residual water contents: f 0 = u s -u r . Finally, r m is calculated from the Young-Laplace equation where r = A/h, where A is a proportionality constant obtained from the variables of the equation, A = À0.149 cm 2 and h [L] is the pressure head. The code for the initial PSD is shown in Fig. 1.

Coefficients of the PDE
The coefficients of the PDE are calculated following the approach of [4] who used moment analysis of the PSD to yield the definitions for mean and variance. Moments are defined by integrating the PSD  with respect to the pore size: Normalized moments (M n ) are calculated through division by the zero-order moment (m 0 ). The code to obtain the zero order moment is shown in Fig. 2.

Degradation term
The degradation term is expressed by means of an exponentially decaying function [9]: Here, c and d are empirical coefficients which are obtained from the zero-order moment. In our study using moment analysis, we evaluated m 0 values to check if the probability was preserved. If the  probability was not preserved, we included the degradation term in the analytical solution. The code for the degradation term is shown in Fig. 3.

Drift term
The first-order normalized moment M 1 characterizes the mean pore size hri [L]: The first-order moment can also be obtained from other independent models [1]. For instance, the rather popular expression from Thornley [10] is used for the drift term.
The code for both these scenarios are shown in Fig. 4.

Dispersion term
Finally, the second-order centralized moment m 2 characterizes the variance [11]: However, our lack of knowledge on how the dispersion behaves with respect to the pore size paved the way to use the dispersivity l.

Analytical solution
As a final step, the pore size distribution (analytical solution in Eq. 5) is obtained using the code in Fig. 5a and b.
The optimized lambda value should now be used in the code ( Fig. 5a and b) instead of lam to get the predicted pore size distribution curves. The reader is directed to Chandrasekhar et al. [1] for the output curves.

Outlook
The possibility of using other independent models for the coefficients exists and is being looked into. Further, the optimization process for lambda necessitates the final water retention parameter values making the model redundant. However, lack of sufficient data sets that capture the temporal dynamics of soil structure hampers our quest to establish a range of values for the coefficients as well as for the validation process. In conclusion, the model is able to capture very well the evolution of soil pore size distribution in its current state with scope for improvement in how we estimate the coefficients.