On the three‐parameter infiltration equation of Parlange et al. (1982): Numerical solution, experimental design, and parameter estimation

Soil water infiltration is central to hydrologic studies and lends itself for detailed experimentation and mathematical–physical modeling. The most rigorous of these approaches numerically solve Richards’ equation, possibly coupled through a heterogeneous soil to a surface water routine and groundwater model. This approach is computationally expensive and prone to mass balance errors and overparameterization. Analytic solutions of the infiltration process obviate the need for specification of hydraulic functions and simplify computation and inverse determination of soil properties. This paper investigates the usefulness of Parlange's three‐parameter infiltration equation for forward and inverse modeling of vertical infiltration experiments. This quasi‐exact implicit solution of Richards’ equation, also credited to Haverkamp and coworkers in recent literature, is valid for the entire infiltration event, matches cumulative infiltration data from different soils and its (super‐)parameters, S, Ks, and β, exhibit a solid mathematical–physical underpinning. Nonetheless, Parlange's equation has not entered mainstream use for infiltration simulation and data analysis in the absence of a robust, exact and efficient numerical solution. This paper builds upon the recent work of Jaiswal et al. and presents theory, algorithms, and source codes of two numerical procedures for forward and inverse modeling of Parlange's infiltration equation. We illustrate the procedures using measured infiltration data from the Soil Water Infiltration Global (SWIG) database. Our findings highlight the potential of Parlange's equation for infiltration modeling and hydraulic characterization of the soil sorptivity, S [L T−1/2], saturated hydraulic conductivity, Ks [L T−1], and unitless coefficient, β. Parlange's infiltration equation provides a powerful alternative to mathematically more convenient explicit infiltration equations that suffer physical underpinning and/or a limited time validity.


INTRODUCTION AND SCOPE
Infiltration is the process by which water enters the soil and moves downward under the influence of gravity and/or capillary action, thereby soaking and/or filling up the pore space, replenishing the root zone, and possibly seeping into rocks through cracks. Infiltration is not only an important component of the hydrologic budget but also exerts a large control on surface runoff, erosion, and root zone soil moisture content, which, in turn, governs slope stability, plant water uptake, evaporation, groundwater recharge, and heat exchange between the Earth's surface and the atmosphere (Liu et al., 2011). As a result, infiltration is studied in different fields by water resource and hydraulic engineers, hydrogeophysicists, hydrogeologists, plant scientists, agronomists, and ecologists (Philip, 1969). Soil water infiltration lends itself for detailed experimentation, mathematical-physical modeling, and spatiotemporal analysis. Numerous mathematical-physical approaches have been developed to describe water infiltration into variably saturated soils (Assouline, 2013). The most rigorous of these approaches would use Richards' equation (Richards, 1931), coupled, if necessary, through a heterogeneous soil with a surface water routine and groundwater model. As this modeling approach is computationally expensive and prone to mass balance errors and, possibly, overparameterization, preference is usually given to an analytic solution of the infiltration process. Such analytic functions simplify computation, obviate the need for specification of the water retention and unsaturated soil hydraulic conductivity functions, and allow for a rapid estimation of their unknown coefficients using least squares regression methods. These coefficients may be viewed as super-parameters of the hydraulic functions.
A quintessential example of a super-parameter is the soil sorptivity, [L T −1/2 ], which measures the capacity of a soil to absorb water by capillarity. This parameter is innate to analytic solutions of Richards' equation for horizontal and/or vertical infiltration and may be defined as follows (Parlange, 1975;: where θ, θ i , θ s , and θ r [L 3 L −3 ] denote the current, initial, saturated, and residual volumetric soil moisture contents, respectively, w (θ) [L 2 T −1 ] signifies the soil water diffusivity function of Mualem (1976) and van Genuchten (1980), α [L −1 ] is a proxy of the soil's air-entry value, m is a unitless coeffi-

Core Ideas
• Methods for exact and robust forward and inverse modeling of Parlange's three-parameter infiltration equation are presented. • Methods allow for a rapid inverse estimation of , s , and β from cumulative infiltration measurements. cient, s [L T −1 ] denotes the saturated hydraulic conductivity, a measure of the soil's ability to transmit water under the influence of gravity, λ [−] is a pore-connectivity parameter, and e = (θ − θ r )∕(θ s − θ r ) is the unitless degree of saturation. The sorptivity, , not only depends on the soil's hydraulic properties as characterized by the Mualem-van Genuchten parameters, θ r , θ s , α, = 1∕(1 − ), s and λ, but is also a function of the soil's initial moisture content. Hence, the sorptivity is not an invariant soil property, which obfuscates its use in hydraulic functions. For notational convenience, we omit the explicit functional dependence of the soil sorptivity, , on the initial moisture content.
In this paper, we focus our attention on the three-parameter infiltration equation of Parlange et al. (1982). This quasi-exact implicit solution of Richards' equation describes cumulative vertical infiltration, ( ) ≥ 0 [T], into a homogeneous soil at uniform initial moisture content where ≥ 0 [T] denotes time, i [L T −1 ] signifies the soil hydraulic conductivity at the initial moisture content, θ i , and, 0 < β < 2 is a unitless parameter. This range of β values exceeds the unit interval stipulated by Parlange et al. (1982) and satisfies the theoretical definition of β from Fuentes et al. (1992) stated explicitly in Haverkamp et al. (1994). For i = 0 and β → 0, the above expression simplifies to the wellknown infiltration equation of Green and Ampt (1911) and to soil types for which the diffusivity, w (θ), increases much more rapidly with θ than d ∕dθ (Assouline, 2013). As the above expression originates from the water content form of Richards' equation, its use is limited to nonponded conditions. Parlange et al. (1985) and Haverkamp et al. (2015) generalized Equation 2 to ponded conditions with an explicit form presented by Barry et al. (1995). Haverkamp et al. (1994) renewed interests in Equation 2 in pursuit of an analytic solution for three-dimensional infiltration from a disc infiltrometer. Some publications in the vadose zone literature even credit Equation 2 to Haverkamp et al. (1994), reinforcing Stigler's law of eponymy (Stigler, 1980), a supposed tendency of eponymous expressions to be named after people other than their original discoverer(s). As quasi-exact analytic solution of Richards' equation, Parlange's infiltration equation has several desirable characteristics. Equation 2 is valid for the entire duration of the infiltration event, matches cumulative infiltration data from a wide range of soils, and the parameters, , s , and β, exhibit a solid mathematical-physical underpinning, adjustable to the initial and boundary conditions of the infiltration event (Haverkamp et al., 1994;Rahmati et al., 2020). Despite these desirable qualifications, the two-term infiltration equation of Philip (1957), has remained the default choice of many soil hydrologists for analyzing measured cumulative infiltration curves, {̃,̃} =1 , consisting of different (̃,̃) data pairs. This infiltration equation is easy to use, supported by detailed mathematicalphysical analysis, and the optimum values of and the product of the dimensionless multiple and s can be obtained from linear least squares via a closed-form expression. One has to be careful, however, in fitting Philip's two-term equation to measured infiltration data as the function has a limited time validity, valid [T]. This issue is often overlooked and severely compromises the practical applicability of Equation 3. Indeed, we should not use (̃,̃) data pairs beyond valid , as this will corrupt the least squares estimates of and s (Jaiswal et al., 2021). Thus, the choice for Philip's two-term infiltration equation is not given in by its superior ability in describing the infiltration process, but rather by numerical convenience.
Parlange's infiltration equation does not suffer a limited time validity. This benefit comes with several practical advantages and, if nothing else, allows determination of and s from any collection of measured (̃,̃) data pairs. The reason so as to why Parlange's infiltration equation has not entered into mainstream use for infiltration modeling and data analysis was articulated in an earlier publication by Jaiswal et al. (2021). The designation of time, , as an independent variable in Parlange's infiltration function honors the causal relationship between and but prevents a closed-form solution of Equation 2 in pursuit of the ( ) relationship. Indeed, at each time , the cumulative infiltration, , needs to be determined iteratively so that the left-hand side of Equation 2 matches exactly the two terms on the right-hand side. In essence, this amounts to an inverse problem and necessitates a proper numerical solver. Unfortunately, most published papers about Parlange's infiltration equation in the vadose zone literature did not document a numerical solution to the ( ) form (Kargas & Londra, 2020;Lassabatere et al., 2009Lassabatere et al., , 2014Wang et al., 2017), nor consider a recipe for parameter estimation. Other researchers turned their attention instead to the simplified solution presented by Haverkamp et al. (1994) and implement three separate equations for cumulative infiltration at very short, short, and long times, respectively (Bouarafa et al., 2019;Di Prima et al., 2020;Haverkamp et al., 1994;Liu et al., 2018;Moret-Fernández et al., 2020;Rahmati et al., 2020;Stewart & Abou Najm, 2018;Vandervaere et al., 2000). These three equations are much easier to evaluate yet trade numerical issues of Equation 2 with problems of time validity and continuity.
In an earlier publication, Jaiswal et al. (2021) introduced an exact, robust, and computationally efficient solution of Parlange's infiltration equation. This solution uses root finding with the secant method to find the zero-points of the residual form of Equation 2. This amounts to a finite-difference approximation of Newton's method and uses consecutive secant lines to produce a sequence of iterates intended to converge to solutions of the cumulative infiltration, ( ), relationship. Several authors have used the bisection method to find the interval of , wherein the residual function changes sign (Latorre et al., 2015(Latorre et al., , 2018Moret-Fernández et al., 2019). This interval-halving approach is not particularly efficient nor does it optimally exploit the characteristics of the residual function. In this paper, we refine the numerical solution of Jaiswal et al. (2021) and present theory, algorithms, and source codes of two different procedures for forward and inverse modeling of Parlange's infiltration equation. Of the two forward modeling approaches, one involves only explicit equations, which obviates the need for an iterative solution of the ( ) relationship and allows use of the basic features of a spreadsheet for forward modeling of Equation 2. Both numerical solutions admit application of gradient-based parameter estimation with the Levenberg-Marquardt (LM) algorithm (Levenberg, 1944;Marquardt, 1963) for least squares estimation of the Parlange parameters, , s , and β, from measured cumulative infiltration data. We benchmark the two numerical solutions and present the results of our preliminary studies using a large sample of experiments and soils from the Soil Water Infiltration Global (SWIG) database of Rahmati et al. (2018). To the best of our knowledge, this is the first paper that fully exploits the potential of Parlange's infiltration equation for inverse estimation of hydraulic properties. Latorre et al. (2015) used inverse estimation of Equation 2 to infer the soil sorptivity, , and saturated hydraulic conductivity, s , from measured cumulative infiltration data, but they fixed β to the common value of 0.6. A similar approach with β = 0.6 was used in other studies (Dohnal et al., 2010;Latorre et al., 2018;Moret-Fernández et al., 2020;Rahmati et al., 2019).
The remainder of the paper is organized as follows. Section 2 expands on the iterative solution of Jaiswal et al. (2021) and introduces two different numerical solutions of Parlange's infiltration equation, the so-called infiltration and time forms of Equation 2. In this section, we are especially concerned with the partial derivatives of the two forms of Parlange's equation. This is followed by Section 3, which describes briefly the measured infiltration data of the SWIG database, and Section 4, which verifies the accuracy, robustness, and computational efficiency of the two numerical solutions and evaluates their ability to match cumulative infiltration data of a large number of experiments of the SWIG database. In this section, we benchmark the least squares values of the Parlange parameters and their confidence intervals against their counterparts derived from Bayesian analysis using Markov chain Monte Carlo (MCMC) simulation with the DiffeRential Evolution Adaptive Metropolis [DREAM (ZS) ] algorithm with sampling from past states (Vrugt et al., 2008(Vrugt et al., , 2009(Vrugt et al., , 2016. Finally, Section 5 concludes our paper with a summary of our main findings.

PARLANGE'S INFILTRATION EQUATION
The cumulative infiltration function of Parlange in Equation 2 is a quasi-exact implicit solution of Richards' equation. We conveniently refer to this original formulation as the infiltration form of Parlange, in anticipation of an alternative formulation, the so-called time form, detailed in Section 2.2. In this section, we present numerical solutions of both forms of Parlange's infiltration function and introduce recipes for least squares estimation of their parameters, including the sorptivity, , saturated hydraulic conductivity, s , dimensionless coefficient, β, and/or hydraulic conductivity, i , at the initial moisture content, θ i . In Section 4, we illustrate the two forms of Parlange's infiltration function by application to measured infiltration data from the SWIG database of Rahmati et al. (2018). As in Jaiswal et al. (2021), we follow the theoretical definition of β from Fuentes et al. (1992) and Haverkamp et al. (1994) and assume that 0 < β < 2 and β ≠ 1.

Forward modeling
The classification of Parlange's infiltration function as a semiimplicit equation may seem superfluous but nevertheless, has important practical implications, as it means that we cannot simply compute the cumulative infiltration, , as a function of time, . Rather, at each time, , we must use an iterative recipe to determine the value of so that the left-hand side of Equation 2, ( , , s , i , β), matches exactly its right-hand side, ( , , , s , i , β), as follows: , We can express the above equation in two new variables, Δ = s − i [L T −1 ], and ξ = Δ ∕ 2 [L −1 ]. This change of variables simplifies the definition of the residual function: If we divide all terms by ξ > 0, then we yield which may be further rearranged to The wording residual function is appropriate as ( , , Δ , i , β, ξ) measures the deviation between an observable true value (zero) and the estimated value of this quantity. At any given time, , the zero-point of ( , , Δ , i , β, ξ), equals a solution, ( ), of the cumulative infiltration curve. Thus, we must find the value of the cumulative infiltration, , at which the residual function equals zero. To provide insights into the functional shape of the residual function, please consider Figure 1, which presents a plot of ( , , Δ , i , β, ξ) at = 3 h using = 2.0 cm h −1/2 , i = 0.0 cm h −1 , s = 1.0 cm h −1 , and β = 1.5. The red cross at the intersection of the horizontal and vertical dashed red lines equals the root of the residual function, ( , , Δ , i , β, ξ).
To understand whether the properties of the residual function generalize to other Parlange parameter values and measurement times, , Figure 2 visualizes the residual function, ( , , Δ , i , β, ξ), for 12 USDA soil types using i = 0 and values of , s , and β documented in Table 5 of Jaiswal et al. (2021) and shown in Table 1. The right-most column reports the time, 5 , it takes for each soil type to infiltrate = 5 cm of water. For each soil type, we consider = 10 equally spaced measurement times, , between = 5 ∕ and = 5 . This sequence of 10 measurement times is color coded with a palette of light to dark blue and matches the color bar of the clayey soil in Figure 2a. The different graphs confirm our earlier findings. The residual function, ( , , Δ , i , β, ξ), is smooth, continuous, and strictly decreasing for all soil types but (d) loamy sand and (e) sand. Specifically, ( , , Δ , i , β, ξ) increases with if β > 1 and decreases for β < 1. This strict monotonicity proofs the existence of a unique root for all soil types and measurement times. Support for this claim is provided by the different graphs, which demonstrate that the residual function crosses the y axis only once within the domain, 0 ≤ ≤ 6 cm, of cumulative infiltration values, . These properties of the residual function are highly desirable and will help guide a search method to the zero-point of ( , , Δ , i , β, ξ) from any initial guess, (0) , of this root. This guarantees a robust and exact numerical solution of the infiltration form of Parlange's equation.
In the remainder of this paper, we suppress the dependence of the residual function, ( , , Δ , i , β, ξ), on soil-specific constants, Δ = s − i , i , β, and ξ = Δ ∕ 2 , and write ( , ) instead. The residual function in Equation 9 warrants application of an arsenal of root-finding algorithms to solve for a collection of different ( , ) data pairs that make up the cumulative infiltration curve, ( ). Consider the residual function,  Table 1 below with = 20 measurement times each. Each time, , is coded with a different blue tint to b . Latorre et al. (2015) used the bisection method to find the interval of , wherein the residual function changes sign. This interval-halving approach is an improvement over bracketing and uses the function value, ( c , ), at the midpoint, c = a + ( b − a )∕2, of the root interval, [ a , b ], to verify in which of the two subintervals, [ a , c ] or [ c , b ], the function changes sign. If ( a , ) and ( c , ) have opposite signs, the subinterval [ a , c ] must contain a zero-point, and, consequently, we set b = c ; otherwise, there must be a sign change in the interval [ c , b ], and we set a = c . We then compute the midpoint of the new [ a , b ] interval, and so forth. Thus, at each successive iteration, bisection selects the subinterval (left or right) with opposite signs, thereby reducing the width of the root interval by 50%. This process is continued until the interval (bracket) surrounding the root is sufficiently small. The bisection method may be easy to use but requires a relatively large number of residual function evaluations to determine an accurate value of its zero-point at each time .
In our earlier publication (Jaiswal et al., 2021), we used the secant method to locate the root of ( , ). This approach exploits much better the monotone nature of the residual function and, as a result, demands only a handful of iterations to find its zero-point. In this paper, we further refine our iterative procedure and resort to Newton's method instead: where ′ ( , ) signifies the derivative of the residual function with respect to and denotes iteration counter. The derivative, ′ ( , ), of the residual function, ( , ), will help determine the direction of movement to the zero-point. For example, if [ ( ) , ] < 0 and ′ [ ( ) , ] > 0 at the kth iterate, ( ) , then the root of the residual function must lie to the right of ( ) , and the next iterate, ( +1) , should exceed, ( ) . Equation 10 satisfies this analogy. The derivative of the residual function, ′ ( , ), can be derived by analytic means: We now have all ingredients to implement Newton's method and solve for the root of Parlange's residual function at any time, > 0. Figure 3 illustrates the application of Equation 10 to the residual function, ( , ), at = 1 h using = 2.0 cm h −1/2 , i = 0.0 cm h −1 , s = 1.0 cm h −1 , and β = 1.5. Suppose our initial guess of the root, (0) , is (0) = 0.45. Based on the value of the residual function, ( , ), and its derivative function, ′ ( , ), at = (0) , the recurrence relationship of Equation 10 yields a new guess, (1) ≈ 4.956, of the zero-point of ( , ). If we repeat this recipe for (1) , then we yield our next iterate, (2) = 2.485. After one more iteration, we attain (3) = 2.279, in close proximity of the actual root, r = 2.274, of ( , ). Thus, after only a handful of iterations at = 1 h, we yield the approximate zero point of the residual function, ( , ), a testament to the efficiency of Newton's method. In practice, however, we may need more iterations to improve solution accuracy and satisfy convergence thresholds specified by the user. This may include termination tolerances on the root, tol , and its function value, tol , which promote solution accuracy and precision, and conditions such as the maximum allowable number of iterations, max , which safeguard computational efficiency. Note that the distance to the root and/or magnitude of the residual function, ( , ), do not always decrease in the sequence of iterates, (0) , (1) , …, ( ) . Indeed, the solution may deteriorate temporarily between successive iterations en route of Newton's method to the zero point of the residual function.
In the examples considered herein, the first iterate, (1) , is worse than (0) , or | [ (1) , ]| > | [ (0) , ]|, but subsequent solutions improve upon the initial guess. Smoothness and monotonicity come in handy as these properties of the residual function virtually guarantee that the sequence of iterates from Newton's method will converge to the zero-point of ( , ), independent of the choice of the initial guess, (0) . The constant sign of the derivative function, ′ ( , ), will help guide the iterates to the root, r , of ( , ) from any reasonable initial guess, (0) > 0. Care should always be exercised that the residual function, ( , ), is differentiable in the interval around the root.
If so desired, we can substitute the expressions of ( , ) and ′ ( , ) in Equation 10. This produces the following recurrence relation: We do not use this equation in our numerical solution but rather evaluate Equation 10  > −1 for all = (2, 3, … , ). The user is free to specify the termination tolerances on the function value of the root, tol > 0, and the maximum allowable number of iterations, max , of Newton's method. Our numerical experiments have shown that the default values of tol = 10 −12 and max = 20 guarantee a robust, accurate and efficient numerical solution of the ( ) relationship. To simplify application and use, Supplemental Section S1 presents the implementation of Algorithm 1 in MATLAB. The subroutine Parlange_I.m may be executed directly from the MATLAB command prompt by typing I = Parlange_I(eta,t) and is equal to a vector-valued function, = (η, ), of the × 1 parameter vector, = [ s β i ] ⊤ , and the × 1 time vector, = [ 1 2 ⋯ ] ⊤ . The values of tol and max may be submitted as third and fourth input arguments in the function call to Parlange_I, otherwise their default values of 10 −12 and 20, respectively, will be assigned instead.
We would be remiss not to comment on two important aspects of our numerical solution of Parlange's infiltration form. In the first place, the choice of the initial guess, (0) [named (0) in Algorithm 1], and its impact on algorithmic robustness and efficiency has been studied to considerable extent. Details of these studies are presented in Supplemental Section S2. Secondly, the exponential and logarithmic . This limit of the exponential function may already be surpassed at early infiltration times for large values of ξ = Δ ∕ 2 . Very small values of ξ, on the other hand, may provoke numerical underflow. These numerical artefacts are unavoidable but can be evaded to some extent by imposing a lower limit of 10 −4 on the values of , s , and β (see Table 2).
As i does not elicit numerical problems, its lower bound may be set to zero. According to the mathematical definition of β given by Haverkamp et al. (1994), this dimensionless quantity cannot exceed two. Furthermore, for β = 1, the residual function, ( , ), will equal zero for all > 0. As a result, we yield, β ∈ [10 −4 , 1) and β ∈ (1, 2). The upper bounds of the soil sorptivity, , and initial and saturated hydraulic conductivities, i and s , respectively, are diffuse and not as clearly defined. Their values should just be set large enough to accommodate all soil types of the textural triangle. We will revisit this topic of floating-point arithmetic in a later section and present a simple workaround (patch) that rectifies simulation of the ( ) relationship at late times and/or for unusual parameter values.
where ( ,̃) is the scalar-valued form of Parlange's infiltration Equation 2, ( ,̃), signifies the × 1 vector of cumulative infiltration residuals, , and the symbol ⊤ denotes transpose. The formulation, ( , ), implies that time, ≥ 0, is treated as exogenous (independent) variable and, consequently, the function will return as output argument the cumulative infiltration, ≥ 0. The vector-valued form of Equation 2, ( , ), returns a vector of cumulative infiltration values, = [ 1 2 ⋯ ] ⊤ , corresponding to the entries of the time vector, = [ 1 2 ⋯ ] ⊤ . Figure 5 illustrates the application of the infiltration form of Parlange's equation to a data record of = 10 measured (̃,̃) data pairs. The vertical lines portray the residuals, (η,̃) = [ 1 (η,̃1) 2 (η,̃2) ⋯ 10 (η,̃1 0 )] ⊤ . In Jaiswal et al. (2021), we used the DREAM (ZS) algorithm to find the optimum values of the Parlange parameters and their associated posterior uncertainty. This MCMC method uses repeated evaluation of the infiltration form of Parlange's equation to delineate the d-dimensional parameter space of statistically acceptable solutions. In this paper, we focus our attention on the class of local search methods, which seek iterative improvement starting from a single initial point, (0) , in the d-dimensional parameter space, ∈ ⊆ ℝ . Well known methods such as Newton, Gauss-Newton, and LM rely on a so-called Hessian (sensitivity) matrix, ( ), and gradient vector, (η), to find the most productive search direction in pursuit of the objective function minimum. These methods share in common their expression for the new iterate, ( +1) , of Parlange parameter values but differ in their computation of the × 1 shift vector, Δ ( ) . We use the LM algorithm (Levenberg, 1944;Marquardt, 1963), also known as damped least squares: where f [ ( ) ,̃] denotes the × Jacobian matrix of Parlange's infiltration form, [ ( ) ,̃], is the × 1 vector of cumulative infiltration residuals for the current iterate, ( ) , of Parlange parameter values and, Λ ( ) , is the × damping matrix where λ > 0 is a strictly positive damping factor and the binary operator ⊙ works as follows: and produces a × diagonal matrix with entries λ ( ) ∀ = (1, 2, … , ) on the main diagonal and zeros elsewhere. The damping factor, λ, is strictly positive and adjusted after each iteration, .
Thus, by adaptively changing the value of λ between successive iterations, the LM method can switch between Gauss-Newton (uses curvature of response surface) and steepest descent (uses gradient only).
The × Jacobian matrix, f ( ,̃), is a generalization of the gradient vector, ( ,̃), to vector-valued functions and stores the first-order derivatives of the infiltration form of Parlange's Equation 2 with respect to its d parameters, , s , i , and β: is the vector-valued form of Equation 2 with the measured infiltration times,̃, as exogenous variable, and = (1, 2, … , ). Thus, the jth column, f , of f ( ,̃), measures the derivative of ( ,̃1), ( ,̃2), …, ( ,̃), with respect to the jth parameter, η , with the other − 1 parameters held fixed at their respective values. For the time being, we treat the initial hydraulic conductivity, i , as an unknown parameter. In practice, however, this parameter is often set to zero. This assumption may not seem appropriate for soils with a rather high initial moisture content, nevertheless, is justified by our earlier findings presented in Jaiswal et al. (2021).
The implicit infiltration form of Parlange's function in Equation 2 admits analytic expressions for its partial derivatives. Symbolic differentiation (see Appendix A) leads to following expressions for the parameters: and an equation for the infiltration rate, ∂ ∕∂ : where α is a dimensionless variable The × Jacobian matrix, f ( ,̃), of the infiltration form of Parlange's equation stores the partial derivatives of Equations 19a,19b,19c,and 19d, as follows: with column-wise entries written in numerator layout where ( , ) signifies the vector-valued form of the infiltration expression in Equation 2. The analytic form of the Jacobian matrix, f ( ,̃), can thus be computed directly and without difficulty, thereby avoiding the need for less accurate and computationally more expensive numerical differentiation in our search for the most productive descent direction, Δ ( ) , at = ( ) . Figure 6 visualizes the four columns of the Jacobian matrix, f (η,̃), for 0 ≤ ≤ 20 cm, using = 2.0 cm h −1/2 , s = 1.0 cm h −1 , i = 0.1 cm h −1 , and β = 1.5.
Notice the excellent agreement between the analytic values of (a) ∂ monotonically with the cumulative infiltration. The partial derivatives of the soil sorptivity, , and initial hydraulic conductivity, i , have a striking similarity. The two functions increase synchronously from zero, reach a maximum magnitude after infiltrating about = 10 cm of water, and stay constant thereafter. This sigmoidal relationship of the partial derivatives with respect to is not innate to ∂ ∕∂ and ∂ ∕∂ i but is also observed for the parameter β, albeit opposite in magnitude and sign. Indeed, the partial derivative, ∂ ∕∂β, decreases monotonically from zero at = 0 to reach a constant magnitude of about −0.65 halfway through the infiltration experiment. The partial derivative of the saturated hydraulic conductivity, ∂ ∕∂ s , is negligible at early times of infiltration but increases linearly with and thereafter to attain relatively large values at late infiltration times.
The Jacobian matrix is not only a principal input of the LM algorithm but also conveys important information about Parlange parameter sensitivity. Indeed, the larger the magnitude of ∂ ( , )∕∂ , ∂ ( , )∕∂ s , ∂ ( , )∕∂β, and ∂ ( , )∕∂ i , with respect to measured cumulative infiltration data,̃, the easier it is to pinpoint the values of , s , β, and i . From the partial derivatives of , β, and i , we conclude that cumulative infiltration measurements at intermediate to late infiltration times contain most information about the soil sorptivity, , coefficient β, and initial hydraulic conductivity, i . As the partial derivative of s continues to grow with time, late infiltration measurements will be most informative for the saturated soil hydraulic conductivity, s . We conclude that the sensitivities to the Parlange parameters are not segregated in time but rather congregate at intermediate to late infiltration times. This finding, albeit useful for experimental design, explains the strong correlation between the posterior samples of s , s , and β reported in our earlier study (Jaiswal et al., 2021) using synthetic infiltration data simulated with HYDRUS-1 D. This correlation not only increases parameter uncertainty but also complicates the search for the least squares values of , s , and β. What is more, as ∂ ∕∂ , ∂ ∕∂ s , ∂ ∕∂ i , and ∂ ∕∂β are maximized at intermediate to late infiltration times, this places a premium on the length of the infiltration experiment. As most infiltration experiments typically last a few hours only (see Figure 7), the measured dataset, {̃,̃} =1 , of (̃,̃) pairs may not contain enough information to warrant an accurate determination of the sorptivity, initial and saturated soil hydraulic conductivities, and coefficient β. Certainly, short infiltration experiments will provoke problems with an accurate estimation of the Parlange parameters.
A detailed algorithmic recipe of the LM optimization method is presented in Supplemental Section S3. A comparison of the objective function at the proposed point, p , and the current iterate, ( ) , is used to adjust the value of λ. Furthermore, we describe an amendment to the LM algorithm that guarantees that the least squares estimates of , s , β, and i stay within the parameter space delineated by the half-closed and closed intervals of Table 2.
In the next section, we will discuss an explicit solution of Parlange's infiltration equation. This so-called time-form has a direct analytic solution and serves as a prelude to Section 4 in which we will compare the numerical solutions of the infiltration and time forms of Parlange's infiltration equation and benchmark their performance against measured data from the SWIG database.

Forward modeling
The infiltration form of Parlange's infiltration equation, ( , ), in Equation 2 requires a root-finding method to solve for the cumulative infiltration, , at each time, . Our iterative procedure with Newton's method, albeit exact, efficient, and fast, may not convince all readers to use this numerical solution of Parlange's infiltration equation for forward modeling and parameter estimation. We, therefore follow Lassabatere et al. (2009) and Kargas and Londra (2020) and consider an alternative solution of Equation 2. This so-called time form of Parlange's infiltration equation does not demand an iterative procedure for forward modeling of its ( ) relationship for given values of , s , and β. This obviates the need for a programming language such as MATLAB and allows use of a spreadsheet for infiltration modeling.
If we set the hydraulic conductivity, i , at the initial moisture content, θ i , equal to zero, then Parlange's infiltration function in Equation 2 simplifies to The designation of time, , as an independent variable in the above expression is common to science and engineering yet necessitates use of an iterative recipe such as articulated in Section 2.1.1 for simulation of the ( ) relationship. This includes a matching notation of ( , ) and ( , ) for the scalar-and vector-valued forms of Parlange's infiltration function and extraneous variables, = [ s β] ⊤ . With a swap of the dependent and independent variables, the above formulation of Parlange's infiltration equation admits a direct solution for the infiltration time, , as follows: , if the input argument, 2βξ , of the exponential function exceeds the value of 709.783. We utilize the remedy of Table 2 and specify a lower limit of 10 −4 on the values of , s , and β. An antidote to numerical under-and overflow of the time and infiltration form of Parlange's equation will be presented in the results section.

Inverse modeling
The parameters, , s , and β, of the time form, ( , ), of Parlange's infiltration Equation 2 may be estimated from a measured record, {̃,̃} =1 , of cumulative infiltration values. Supplemental Section S5 details our implementation of the LM algorithm for constrained least squares estimation of the parameters , s , and β of the time form of Parlange's infiltration function. This algorithmic recipe is similar to that used for least squares estimation of , s , β, and i using the infiltration form of Parlange's equation in Supplemental Section S3 but with a different formulation of the objective function and Jacobian matrix. The time form of Parlange's equation treats time, (and not as in the infiltration form), as the dependent variable. This swap of the dependent and independent variables implies the use of time residuals in evaluating the goodness of fit of the simulated infiltration curve (see Figure 9). If we make the convenient assumption that the time residuals are uncorrelated and with a constant variance, then we yield the least squares objective function: The × Jacobian matrix of the time form of Parlange's infiltration equation, f ( ,̃), organizes column-wise the partial derivatives of ( ,̃) with respect to , s , and β as follows: We can simplify the above equations by expressing the partial derivatives as a function of Δ = s and ξ = Δ ∕ 2 . This change of variables would imply different entries of the parameter vector, , which is not desired in our formulation of the LM method. To provide insights into the different columns of f ( , ), please consider Figure 10  First and foremost, notice the excellent agreement between the analytic values of ∂ ∕∂ , ∂ ∕∂ s , and ∂ ∕∂β and their estimates derived from numerical differentiation. The analytic derivatives are preferred in practice as they are not only exact, but also easy to compute. Secondly, the traces of the partial derivatives of , s , and β are in close agreement with their respective counterparts of the infiltration form of Parlange's equation depicted previously in Figure 6. A side-by-side comparison of ∂ ( , ∕∂η ) and ∂ ( , ∕∂η ) for = (1,2,3) demonstrates that the partial derivatives of , s , and β hardly differ in magnitude but only in sign. The traces of ∂ ∕∂ and ∂ ∕∂β exhibit the characteristic "S"-shaped behavior of ∂ ∕∂ and ∂ ∕∂β but with opposite signs. Their magnitude increases monotonically from zero at early times of infiltration and reach a constant magnitude after infiltrating approximately 10 cm of water. The partial derivative of s is negligible at early times of infiltration but decreases linearly with and thereafter to attain relatively small values at late infiltration times. Lastly, the partial derivatives suggest that time measurements taken at intermediate to late times contain the most information about the Parlange parameters. This overlap confirms reported problems with the joint estimation of , s , and β and reiterates the importance of a sufficiently long infiltration experiment. Now that the forward and inverse solutions of the infiltration and time form of Parlange's equation have been presented, we move on to the experimental data that are used to benchmark both solutions.

EXPERIMENTAL DATA
We evaluate the infiltration and time forms of Parlange's equation using measured data from the SWIG database of Rahmati et al., 2019). This database documents more than 5,000 infiltration datasets involving experiments on all 12 soil types (listed in Table 1)

BENCHMARK ANALYSIS
Before we can proceed with application of the infiltration form of Parlange's infiltration equation to the experiments of the SWIG database, we must first benchmark its numerical solution. This will inspire confidence in the accuracy of the simulated cumulative infiltration values. The time form of Parlange's infiltration equation comes in handy as its numerical solution of the ( ) curve is explicit and, thus, exact. We set i = 0, and draw at random = 1,000 parameter vectors from a = 3-dimensional cube using Latin hypercube sampling with lower limits of 10 −4 (see Table 2) and upper bounds of 50, 100, and 2 for , s , and β, respectively. Next, we create a vector of = 100 cumulative infiltration values spaced equally between = 0.01 cm and = 10 cm. We coin this vector, true = [0.01 0.02 … 9.99 10.00] ⊤ , and repeat the same steps for each parameter vector. Specifically, we use true as input to the time form of Parlange's Equation 26 and store the corresponding infiltration times in the vector true ; thus, true = ( , true ). Next, we admit the vector true to the infiltration form of Equation 4 to yield a × 1 vector of simulated cumulative infiltration values, sim = ( , true ). The numerical error of each parameter vector, err, is now set equal to the Manhattan distance or 1 norm of the true and simulated cumulative infiltration curves, err = ∑ =1 | true, − sim, |, in units of cm. Figure 11 presents a frequency distribution of the logarithmic err values (base 10) of the thousand parameter vectors. The frequencies are divided by the sample size to yield normalized values between 0 and 1.
The Manhattan distance of the true and simulated cumulative infiltration curves is almost always smaller than 10 −6 cm. This provides support for the claim that our numerical solu-F I G U R E 1 1 Histogram of the numerical error of the infiltration form of Parlange's equation. We use a termination tolerance of tol = 10 −12 on the absolute function value at the root and max = 20 iterations in Newton's method tion of Equation 4 is exact and robust. The solution is also efficient requiring only a few milliseconds to complete a single ( ) relationship on an Intel Core i7-10700 T @ 2.00 GHz desktop computer. Note that in about 2% of the cases (about 20 parameter vectors), the MATLAB codes of the infiltration and time form of Parlange's equation terminate prematurely. This comes as no surprise. The argument, 2β s ∕ 2 , of the exponential function in the infiltration (and time) form of Parlange's equation can exceed the upper limit of 709.783 on a 64-bit computer before the simulated infiltration reaches true = 10 cm. This numerical overflow of the exponent results in not-a-number for the respective ( ) and ( ) data pairs.
To provide insights into these so-called dissident (or aberrant) parameter vectors, Figure 12 presents scatter plots of the F I G U R E 1 2 Scatter plots of the parameter pairs, (a) ( , s ), (b) ( , β), and (c) ( s , β) used in our Monte Carlo experiment, where is sorptivity, s is saturated hydraulic conductivity, and β is a unitless coefficient. Dissident parameter vectors are highlighted with a blue color and return not-a-number for one or more entries of the vectors of infiltration time, , and cumulative infiltration, . The bottom-right graph (d) displays the simulated infiltration curves of the dissident parameter vectors bivariate parameter samples of (a) ( , s ), (b) ( , β), and (c) ( s , β) used in our Monte Carlo experiment. The blue-colored squares correspond to dissident parameter vectors and, thus, those that are subject to numerical overflow.
The aberrant parameter vectors have a common denominator. They all congregate at the lower bound of the soil sorptivity, , irrespective of the values of s and β. The closer the sorptivity is to zero, the larger the value of the exponent, 2β s ∕ 2 , and the more susceptible the forward, ( ) or ( ), solution will be to numerical overflow. As the exponent grows linearly with , numerical overflow is inevitable. The infiltration experiment just needs to be long enough. A simple workaround may help salvage the simulated ( ) or ( ) relationships of Parlange's infiltration equation. The cumulative infiltration curves of the dissident parameter vectors depicted in the bottom-right graph (Figure 12d) exhibit a constant slope well before the time of overflow, o . In other words, the soils of the dissident parameter vectors are fully saturated prior to exponent overflow. We can exploit this finding and protect Parlange's infiltration equation against numerical over-flow by switching to saturated flow at = o . The MATLAB functions, Parlange_I_patch and Parlange_t_patch, in Supplemental Section S6 implement this workaround for the infiltration and time form of Parlange's equation, respectively. These patched functions complete the ( ) or ( ) relationships of all thousand soils (parameter vectors). As overflow is not prevalent in the short experiments of the SWIG database, we do not make use of the patched functions. To prevent numerical overflow, one can also discard the second term, β − 1, in the exponent of Equations 9 and 26. This provides a nearly exact solution of ( )or ( ) for ≥ , equivalent to our numerical solution.
We now move on to the collection of 646 measured infiltration datasets and use the LM algorithm to determine the least squares values of , s , and β for the infiltration and time forms of Parlange's equation. We execute Algorithms S3.1 and S5.1 10 times using different initial guesses of , s , and β. In the absence of detailed information about the upper values of the soil sorptivity and saturated hydraulic conductivity, we did not impose upper limits for and s in our successive trials of the LM algorithm. As a result, the optimum Parlange parameters will satisfy, ≥ 10 −4 cm h −1/2 , s ≥ 10 −4 cm h −1 , 10 −4 ≤ β < 2, and β ≠ 1. Figure 13 compares the least squares values of the (a) soil sorptivity, , (b) saturated hydraulic conductivity, s , and (c) coefficient, β, derived from the infiltration (x axis) and time (y axis) forms of Parlange's infiltration equation. For convenience, we assign the subscripts "I" and "T" to the parameters of the infiltration and time forms, respectively. Each square corresponds to a different infiltration experiment. The solid black line going from left to right across the graphs equals the 1:1 line of the plotted quantities.
It is important to keep in mind that we did not inspect nor appraise the infiltration data. Thus, the three scatter plots present only our preliminary findings for the collection of vertical infiltration experiments of the SWIG database. The least squares values of and s center on the 1:1 line and exhibit a relatively small dispersion. This is particularly true for the soil sorptivity and inspires confidence in the ability of the LM algorithm to converge on robust values of the Parlange parameters. The ( s,I , s,T ) pairs in the bottom-left corner of the scatter plot demonstrate considerable variation around the 1:1 line, as is emphasized in the small inset in the top-left corner of the middle graph. It is not immediately clear what causes these differences in the least squares values of the saturated hydraulic conductivity of the infiltration and time forms of Parlange's equation. Without doubt, data quality and quantity exert a large control on the estimates of s,I and s,T , and so will the duration of an infiltration experiment. The scatter of the ( s,I , s,T ) data pairs around the 1:1 line simply articulates different optima of the saturated hydraulic conductivity for the cumulative infiltration and time residuals, respectively. The infiltration and time form of Parlange's equation may differ in their sensitivity to outliers and/or other (̃,̃) data points that breach assumptions of soil homogeneity and/or a constant initial water content. The data of a small number of experiments raise concern as their (̃,̃) relationship violates monotonicity (e.g. Sample IDs 3646 and 3647). The data of some other experiments is at odds with the default dimensions of cm and hour used in the SWIG database. For example, Sample IDs 3710-3745 present high-quality infiltration data, yet their clayey soils supposedly infiltrate between 30 and 50 cm of water in an experiment lasting about 0.05 h (3 min). This raises concerns about data units and calls for a more exhaustive appraisal of data quality. These dissonant soils are not too difficult to single out, as Parlange's equation is unable to describe reasonably well their measured infiltration curves. This is confirmed by the rather large values for the root mean square error of their least squares parameter values, which exceeds by far the typical values of 0.006-0.8 cm (infiltration form) and 2.5 × 10 −4 to 0.1 h (time form). By and large, however, the experimental data of the SWIG database appear robust and accurate. For these soils, discrepancies between s,I and s,T simply articulate a different optimum value of the saturated hydraulic conductivity for the time and cumulative infiltration residuals.
The right-most scatter plot illustrates large differences in the least squares values of β for the infiltration and time forms of Parlange's equation. With the exception of a small cluster of points in the center of the graph, a large proportion of the (β I , β T ) pairs deviate considerably from the 1:1 line and congregate at the edges of the domain. In quite a number of cases, the optimal β I and β T values appear at opposing sides of the parameter space. This is a disturbing result and may signal an insufficient sensitivity of the simulated ( ) relationships to parameter β. However, we should be careful in drawing conclusions about these discrepancies in the inferred β I and β T F I G U R E 1 4 First-order approximation of the parameter confidence intervals of coefficient β for each of the 646 infiltration records of the Soil Water Infiltration Global (SWIG) database. The horizontal and vertical gray lines portray the 95% confidence intervals of the inferred β values for the infiltration (β I ) and time (β T ) form of Parlange's equation, respectively values of the infiltration and time forms, respectively, as the scatter plot portrays only the optimal, least squares, values of β without recourse to their underlying uncertainty.
To provide insights into parameter uncertainty, please consider Figure 14 which plots error bars on the inferred β values. The horizontal and vertical gray lines portray the 95% confidence intervals of the β I and β T values of the infiltration and time form, respectively. These confidence intervals can be computed from the diagonal entries of the covariance matrix, (̂) = (̂, χ)∕( − )[ f (̂, χ) ⊤ f (̂, χ)] −1 , at the least squares solution,̂, of , s , and β, where input argument, χ, signifies the exogenous variable,̃or̃, for the infiltration and time form of Parlange's equation, respectively. The so-obtained confidence intervals are only an approximation of the actual parameter uncertainty.
The least squares β I and β T values exhibit a large uncertainty with 95% confidence intervals that cover a large part of the feasible search space. Indeed, the gray lines (horizontal and vertical) extend far beyond a small neighborhood surrounding the least squares, (β I , β T ), data pairs and go from left to right and top to bottom across the scatter plot. Thus, parameter β appears poorly defined by calibration against the measured infiltration data. This finding does not come as a surprise. The partial derivatives of the infiltration and time forms of Parlange's equation have warned against the use of short infiltration experiments. Such datasets provoke problems with the joint estimation of , s , and β as the information for the Parlange parameters congregates at intermediate to late infiltration times. As a large proportion of the infiltration experiments last a few hours only (see Figure 7), this will complicate an accurate estimation of β I and β T . Similar findings were reported by  for upward infiltration experiments, wherein β governs the exfiltration at long times. In fact, Latorre et al. (2018) concluded that infiltration experiments must exceed 10,000 s (∼2.8 h) to warrant an accurate determination of β. Our findings suggest that 3 h may be optimistic. One should therefore be careful in interpreting the least squares values of β of the infiltration and time forms of Parlange's equation.
To benchmark the least squares estimates and first-order confidence intervals of the Parlange parameters, we separately analyzed the collection of infiltration experiments using Bayesian analysis with the DREAM (ZS) algorithm. This MCMC method explores exhaustively the search space in pursuit of the optimal parameter values and their multivariate posterior distribution. This distribution provides an exact characterization of parameter uncertainty. We use a uniform prior parameter distribution, ( ), and Gaussian likelihood function, ( |χ) = ( , χ) − ∕2 , where ( , χ) signifies the objective function of Equations 13 and 27 for the infiltration, χ =̃, and time, χ =̃, forms of Parlange's equation, respectively. Figure 15 compares the DREAM (ZS)derived values of I , s,I , and β I with their counterparts of Algorithm S3.1 using gradient-based search with the LM algorithm. We observe similar results for the time form of Parlange's equation and thus do not present the findings here.
The optimum soil sorptivity and saturated hydraulic conductivity derived from the DREAM (ZS) algorithm are in nearly perfect agreement with their counterparts of the LM algorithm. The ( DREAM , LM ) and ( s,DREAM , s,LM ) data pairs lay almost exclusively on the 1:1 line, inspiring confidence in the ability of the LM algorithm to locate successfully the least squares values of the Parlange parameters. The (β DREAM , β LM ) data pairs, on the contrary, exhibit considerable dispersion around the 1:1 line. The β DREAM values are found interior to their search space and do not congregate at the edges as their counterparts, β LM , of the LM algorithm do. This manifests our earlier problems with an accurate estimation of coefficient β I . The 95% confidence intervals of DREAM and LM are remarkably similar for the soil sorptivity but differ substantially for the saturated hydraulic conductivity. The first-order approximation characterizes rather poorly the uncertainty of the least squares s,I values of the LM algorithm with 95% confidence intervals that go from top to bottom in the graph and extend far beyond their nonlinear counterparts. Thus, the saturated soil hydraulic conductivity is much better defined by calibration against cumulative infiltration data than one would expect from the linear confidence intervals. Fortunately, the linear and nonlinear confidence intervals are in a much better agreement for parameter β I , as they almost always envelop the search domain. Altogether, we conclude that the LM algorithm successfully infers the least squares values of the Parlange parameters. The firstorder approximation provides an adequate characterization of the confidence intervals of and β but significantly overestimates the uncertainty of the saturated hydraulic conductivity, s . Bayesian analysis coupled with the DREAM (ZS) algorithm will yield an exact characterization of Parlange parameter uncertainty.
We are now left with a demonstration of the quality of fit of the infiltration and time forms of Parlange's infiltration equation. Figure 16 displays the measured infiltration data (red circles) of eight different soil types from our collection of = 646 experiments of the SWIG database, including (a) sand, (b) sandy loam, (c) sandy clay loam, (d) loam, (e) silty loam, (f) clay loam, (g) silty clay loam, and (h) clay. The solid lines portray the least squares fit of the infiltration (in blue) and time (green) forms of Parlange's equation using the values of , s , and β derived from the LM algorithm. The bottom-right corner of each graph reports the SWIG code of each soil. T A B L E 3 Least squares estimates of the parameters, sorptivity ( ), saturated hydraulic conductivity ( s ), and coefficient β, for the infiltration (subscript "I") and time (subscript "T") forms of Parlange's equation Note. Each row corresponds to a different soil type of Figure 16 and documents the length, , of the data record, Soil Water Infiltration Global (SWIG) code, and values of I , s,I , β I , T , s,T , and β T , respectively.
The simulated cumulative infiltration curves are in excellent agreement with the measured data. For some of the soils, we observe small differences between the observed data and the Parlange function at early times of infiltration. These apparent discrepancies are rather small and may be the result of model misspecification, a nonuniform initial water content, and/or heterogeneous soil. Fortunately, this model-data mismatch should not corrupt too much the Parlange parameter estimates as the partial derivatives of the infiltration and time forms are relatively small at early times of infiltration (see Figures 6 and 10), and thus parameter sensitivity is relatively low. Infiltration measurements at late times, on the contrary, will exert a much larger control on the values of , s , and β.
To understand the close correspondence between the timeseries plots of the infiltration and time forms of Parlange equation for the different soil types, please consider Table 3, which reports the least squares parameter estimates of the eight infiltration data sets of Figure 16. The first three columns list the soil type, number of observations, , of the data record, and corresponding SWIG code followed by the optimized values of , s , and β for the infiltration and time forms of Parlange's equation.
The infiltration and time form of Parlange's equation return fairly similar values of the Parlange parameters. This is true for all soil types, except for sandy loam and sandy clay loam. The differences in the parameter estimates do not seem to affect much the simulated ( , ) relationships (see Figure 16), which are in excellent agreement with each other. We have noticed similar findings for other soils of the SWIG database, yet in some cases using values for β I and β T at opposite sides of the parameter space. This is a somewhat troubling finding and demonstrates the lack of sensitivity of the simulated infiltration curve to coefficient β. Notice that if I < T then s,I > s,T . This is true for all tabulated soils and implies the presence of a negative correlation between and s . If the sorptivity increases, then the saturated hydraulic conductivity must decrease to yield an approximately similar cumulative infiltration curve.

SUMMARY AND CONCLUSIONS
Analytic solutions of the infiltration process obviate the need for specification of the hydraulic functions and allow a rapid characterization of soil physical properties using least squares regression methods. Such functions are preferred over numerical solutions of Richards' equation (Richards, 1931), which are computationally expensive, prone to mass balance errors and subject to overparameterization. In this paper, we focused attention on the three-parameter infiltration function of Parlange et al. (1982). This quasi-exact implicit solution of Richards' equation is valid for the entire duration of the infiltration event, describes accurately measured infiltration data from a wide range of soils, and the parameters, , s , and β, exhibit a solid physical and/or mathematical underpinning, adjustable to the initial and boundary conditions of the infiltration event. Despite these desirable capabilities, Parlange's equation has not entered into mainstream use for infiltration modeling and data analysis, as researchers and practitioners struggle with its numerical solution. This paper builds on our recent work published in Jaiswal et al. (2021) and presents theory, algorithms, and source codes of two separate approaches for forward and inverse modeling of Parlange's infiltration equation. The first approach, coined the infiltration form, treats time, , as an independent (exogenous) variable and uses root finding with Newton's method to simulate the ( ) relationship. This iterative solution may not persuade and/or convince all readers to use Parlange's equation for infiltration modeling and data analysis. We, therefore, considered a second and arguably easier numerical solution of Parlange's infiltration equation. The so-called time form evokes an erroneous causal relationship between and , in favor of an analytic solution of the ( ) relationship. This enables use of a spreadsheet for forward simulation. Benchmark experiments confirmed that the two numerical solutions of Parlange's infiltration equation are exact and efficient yet suffer from numerical overflow of the exponential function. The use of a constant infiltration rate at the time(s) of overflow offers robust protection of Parlange's equation and its numerical solutions against numerical overflow.
The infiltration and time form of Parlange's infiltration equation admit application of gradient-based search with the LM algorithm for least squares estimation of the Parlange parameters, , s , and β, from measured cumulative infiltration data. Preliminary investigations, carried out on more than 600 infiltration experiments of the SWIG database, demonstrated a high level of agreement in the optimal values of and s for the infiltration and time forms of Parlange's equation. Furthermore, the so-obtained least squares values of the sorptivity and saturated hydraulic conductivity were shown to correspond almost perfectly with their counterparts of and s derived separately from an exhaustive exploration of the parameter space using MCMC simulation with the DREAM (ZS) algorithm. This inspires confidence in the ability of the LM algorithm to locate globally optimal values of the Parlange parameters. We witnessed considerable differences in the least squares values of β inferred by the two numerical solutions and parameter estimation algorithms. This somewhat troubling finding demonstrates a lack of sensitivity of the simulated infiltration curve to coefficient β, a result provoked by the rather short lengths of the infiltration records in our collection of experiments. Linear and nonlinear confidence intervals of , s , and β were computed and compared to better understand parameter uncertainty. The first-order approximation provided a reasonable description of the uncertainty of and β but was deficient for the saturated hydraulic conductivity. The relatively short lengths of the experiments of the SWIG database and limited duration contribute to the large uncertainty of the inferred s , and β values as the infiltration process is dominated by the capillary effect.
Theory, algorithms, and source codes presented herein should encourage the use of Parlange's equation for infiltration modeling and data analysis and provide an easy enough alternative to functions suffering from a limited time validity and lack of physical underpinning. For those unfamiliar with MATLAB, we provide a spreadsheet with forward and inverse modeling capabilities of Parlange's infiltration equation.

A C K N O W L E D G M E N T S
We greatly appreciate the constructive comments of the anonymous reviewers. The first author acknowledges interactions with Dr. Mehdi Rahmati on Parlange's infiltration equation and the SWIG database and conversations with Profs. Jan Vanderborght and Harry Vereecken from the Forschungszentrum, Germany, on the definition of length and/or times scales for infiltration. This led to definition of the sorptivity as a super parameter of the hydraulic functions and initial conditions. This paper is in memory of Prof. Scott James from Baylor University, TX, whose sudden passing has left a hole in the scientific community. Many of us talk about interdisciplinary research. He actually did it.

C O N F L I C T O F I N T E R E S T
The authors declare no conflict of interest.

D A T A A N D S O F T WA R E AVA I L A B I L I T Y
The data and software are available upon request from the corresponding author (jasper@uci.edu) and can be downloaded at the following URL https://www.dropbox.com/ sh/xpgpe8mu9pquyva/AAASIRsnJU7Q2FO6RAF_MD4Ta? dl=0. This includes an Excel spreadsheet with forward and inverse modeling capabilities of Parlange's infiltration equation.