**1. Introduction**

The recharge of groundwater [1,2] is critical in many aspects, for example, natural environments, industry, and agriculture. Therefore, recharging aquifers is urgen<sup>t</sup> in regions with growing demands [3,4] on water supplies that are the key to the local ecosystem and to economic development [5,6]. By using either natural or artificial methods to conduct the recharge, estimating infiltration is usually subject to both uncertainties and multiple types of errors. Moreover, the problem of rainfall-induced shallow landslides represents the most common natural hazard [7] in some areas of the world. These landslides are activated by intense rainfall events where water infiltration causes an increase of both volumetric content and pore pressure, thus worsening the slope stability in landslide prone areas. Thus, effective and economic numerical models are first needed to simulate the movement of water in the vadose zone, especially for large-scale distributed hydrologic applications over a relatively long period [8,9].

For groundwater infiltration there are advanced reservoir modeling methods based on the Richards equation [10] that can provide exact solutions to estimate infiltration [11,12]. However, these methods are computationally expensive especially for a large geographic area such as a city or a large ranch or farm. To simulate the rain-infiltration process over a number of years, they are complicated and computationally costly since they determine where the water moves in space. The Green–Ampt model [13], due to its computational convenience, is widely used in estimating infiltration parameters and states, such as flux, accumulative water content, and infiltration time [14–17]. Nevertheless, it makes some ideal, although unrealistic assumptions, for instance, it assumes the existence of an abrupt wetting front, and uniform water content behind the wetting front. To avoid these limitations, several contributions have modified this model. They can be classified as experimentally based corrections [18,19], and mathematically or physically based optimizations [20–22]. To avoid the drawbacks of the Green–Ampt model and rapidly determine how well ground water and aquifers are recharged only, the Talbot–Ogden model is proposed [23,24] for estimating large-scale surface water infiltration into various unsaturated soil textures [25] over long periods. Valuable features of this model are the relatively low computational cost and the large-scale applicability. These features are essential for integrating the hydrologic–hydrogeologic model into an integrated model for the identification of water related hazards as well as supporting an early warning system for the reduction of hydrogeological risk.

The Talbot–Ogden model is derived from the unsaturated Darcy's law and conservation of mass for water moving through a variably saturated porous media. It quickly simulates the infiltration in the water content-depth (θ-*<sup>z</sup>*) domain as its new perspective. In Figure 1, the fundamental idea in the Talbot–Ogden model is presented. In this model, bins are constructed by discretizing the moisture content domain as shown in Figure 1a. According to their water content values, they are independently arranged in parallel, not in series. Bins represent a collection of pore sizes corresponding to a specific range of moisture content θ*i* ≤ θ ≤ θ*e* in a soil. Within a particular bin, this range of moisture content can be found throughout the soil over the vertical domain (Figure 1b). Since the model has only one spatial dimension, this assumption is valid regardless of how su fficiently small the horizontal discretization is. The Green–Ampt equation is transformed and applied to compute the depth infiltration independently in each bin. A process called redistribution, which is invoked at every time step immediately after infiltration, governs the horizontal inter-bin flow along the θ-axis according to the capillary pressure associated with each bin. This process will take into account all the saturated bins but not only restricted to the local neighborhoods of di fferent bins. The infiltration and redistribution are respectively driven by gravitational and relative capillary forces in each bin. During the infiltration, the capillary pressure and hydraulic conductivity become dynamic, and the wetting front as in the Green–Ampt model (Figure 1a) may not exist. The discretized water content domain has also been extended to be a ffine multi-dimensional [26] for depicting more complicated pore size distributions. It is an intrinsically mass conservative model that can be applied to various soil textures [8,23,27]. However, its suitability is directly related to the uncertain number of bins because the predicted flux is highly nonlinear with respect to the discretization in the moisture content domain. Therefore, this uncertainty plays an important role under di fferent soil conditions in this model. The convergence test for choosing a proper number of bins by Talbot and Ogden [23] is more rigorously analyzed in this work, and its physical meaning, a greater variety of pore sizes leading to a larger infiltration rate, can be naturally explained from this study. It will also be quantitatively estimated how an infinite water content discretization a ffects the flux variation through an asymptotic analysis by linearly fitting the wetting front. It directly indicates that a particular depth ratio of the deepest to the shallowest bin fronts can maximize the infiltration flux.

**Figure 1.** The mechanism in the Talbot–Ogden model. (**a**) Bins corresponding to different porosity θ-values compared to the wetting front in the Green–Ampt model; (**b**) Different moisture content in a soil over the vertical domain.

In Section 2, the Talbot–Ogden model is first introduced, which allows us to create a fast solver to simulate infiltration using finite water content discretization. A detailed sensitivity analysis is then made on the quantitative change of the flux as a function of the number of bins, and this work is generalized to an asymptotic analysis that gives an upper bound for infiltration flux variation. In Section 3, the infiltration flux variation using this model is firstly estimated from the physical parameters of a variety of real soil textures. Infiltration simulations for examples in both coarse and fine soil textures are presented using different numbers of bins to validate the theoretical analysis. Section 4 is the conclusive part.

#### **2. Theory and Methodology**

#### *2.1. The Talbot–Ogden Model*

Before presenting the Talbot–Ogden model, let us briefly review the Green–Ampt model relevant to it. The Green–Ampt equation [13], which is based on Darcy's law, provides a very simple model to describe the infiltration of water into the subsurface soil. By neglecting the depth of ponded water, the Green–Ampt equation for vertical infiltration is given by:

$$f = K\_s \left( \frac{(\theta\_d - \theta\_i)H\_c}{F(t)} + 1 \right) \tag{1}$$

where *f* is the infiltration rate, *Ks* is saturated hydraulic conductivity, θ*i* is initial moisture content, θ*d* is the maximum moisture content during infiltration, *Hc* is the effective capillary drive at the wetting front, and *F*(*t*) is the cumulative infiltration depth at time *t*.

The Talbot–Ogden model discretizes the entire water content domain into bins that flow in soils based on the porosity [23] and there is only one vertical spatial dimension (Figure 2) in this model. From the initial moisture content θ*i* to the effective porosity θ*e* there are *n* bins indexed by *j* with equal bin width Δθ. The midpoint value (*j* − 1)Δθ + Δθ/2 represented by θ*j* is the moisture content of the *j*-th bin and the depth of its saturated wetting front is *zj*. The *j*-th bin is assumed to be either fully saturated or dry at any depth. The residual moisture content is θ*<sup>r</sup>*. The rightmost saturated bin is θ*d*. The bins between θ*d* and θ*e* are unsaturated but can become saturated later. In the Green–Ampt model, the

cumulative infiltration function *F*(*t*) is defined as *F*(*t*) = *<sup>z</sup>*(θ*d* − θ*i*). If we let *f* = *dF*/*dt* = (θ*d* − θ*i*)*dz*/*dt*, the vertical infiltration formula by substituting this expression into Equation (1) is obtained.

$$\frac{dz}{dt} = \frac{1}{\left(\theta\_d - \theta\_i\right)} \left(\frac{K\_s H\_\ell}{z} + K\_s\right) \tag{2}$$

**Figure 2.** The infiltration in the Talbot–Ogden model.

Based on Equation (2), the front of the *j*-bin grows as:

$$\frac{dz\_j}{dt} = \frac{1}{\theta\_d - \theta\_i} \left( \frac{K(\theta\_d)\psi(\theta\_d)}{z\_j} + K(\theta\_d) \right) \tag{3}$$

Here, ψ represents the capillary pressure.

The vertical infiltration of water in each bin is governed by Equation (3). The horizontal movement of water through bins is shown in Figure 3. By Equation (3), the infiltration rate *dzj*/*dt* is inversely proportional to *zj* if the other parameters are constant for the *j*-th bin. Hence, bins to the right tend to have greater front depths than the left ones especially in the beginning of infiltration (Figure 3a).

**Figure 3.** The redistribution in the Talbot–Ogden model: (**a**) pre-redistribution and (**b**) post-redistribution.

Due to the capillary effect in soil, water in bins with large θ-values tends to flow to those bins with smaller θ. This horizontal flow is referred to as the redistribution. It reorganizes the redundant water collected from those protruding wetting fronts in all saturated bins proportional to the values of the capillary pressure ψ(θ*j*) of every bin participating in the redistribution [23]. In the redistribution process (Figure 3b), the last deeper bin is defined as the first saturated bin found over the moisture content domain along the negative direction of the θ-axis whose front depth is higher than the shallowest saturated bin. Note that the wetting front depths gradually decrease from left to right after the redistribution since the capillary pressure in the bins on the left acts immediately on water found at depth in the bins to the right in this model [23].

The reasons for choosing a fixed moisture content θ*d* for both *K*(θ) and ψ(θ) in Equation (3) are found in [23]. Water tends move downward through the saturated bins with large θ-values, which is the reason that *K*(θ*d*) is chosen for every bin. However, the capillary pressure in a soil with only lower moisture content than the current saturated bin is always satisfied prior to the rightmost saturated bin, which is the reason that ψ(θ*d*) is chosen for every bin. The functions *K*(θ) and ψ(θ) are from Brooks and Corey [28], but other soil hydraulic models [29–31] can be used without affecting the analysis and conclusions in this paper.

Due to the uncertain discretization of the moisture content domain into *n* bins, a dynamical system [32] is generated

$$z\_{j}^{t+\Delta t} = z\_{j}^{t} + \frac{dz\_{j}^{t}}{dt} \times \Delta t + r \text{dist}\_{j}^{t}, j = 1 \ldots n, t \in [0, T], \tag{4}$$

and

$$\text{redist}\_{j}^{t} = \text{redist}\_{j}^{t} \Big( \psi(\theta\_{1})\_{\prime} z\_{1\prime}^{t} \psi(\theta\_{2})\_{\prime} z\_{2\prime}^{t} \dots \psi(\theta\_{\text{last\\_deper\\_bin}})\_{\prime} z\_{\text{last\\_deper\\_bin}}^{t} \Big) \tag{5}$$

In Equations (4) and (5), the term redist represents the redistribution as the inter-bin movement of water.

The second term of the right hand side in Equation (4) comes from Equation (3) as infiltration. The computation is inexpensive using Equation (4), as for every time step only *n* ordinary different equations are solved for the infiltration. The simulation of redistribution takes *O*(*n*) operations for all *n* bins. In the numerical convergence test by Talbot and Ogden, the largest *n* is 400 with a time step size 2.5 *s* for the coarsest soil sand, which provides a fast, accurate solution.

#### *2.2. Instantaneous Infiltration Rates Analysis*

The most important quantity in the Talbot–Ogden model is the bin width Δθ. Therefore, the range of Δθ must be selected to match the unsaturated flow in specific soil systems because the flux is nonlinear with respect to Δθ. There are two simple assumptions made before analyzing the flux: (1). θ*d* is assumed to be fixed and independent of the number of bins, which is reasonable since θ*d* corresponds to the rightmost saturated bin that is changing during the infiltration; (2). All bins with θ*i* ≤ θ ≤ θ*d* are already saturated. If there are empty bins, it means that the surface water can be absorbed in the next time step so that the instantaneous infiltration rate equals the precipitation rate.

The number of bins and its nonlinear influence on the predicted flux can be considered under these assumptions. Now the only uncertain quantity in the Talbot–Ogden model is the number of bins, which means either a finer or a coarser discretization. All other parameters remain constant. The unsaturated flow divides into two steps within every time step: infiltration and redistribution (Figure 4).

**Figure 4.** A numerical example: (**a**) infiltration and (**b**) redistribution.

The infiltration governs the vertical downward movement of water in every bin. The redistribution governs the horizontal inter-bin flow. The amount of water redistributed to each bin is proportional to the local capillary pressure. In Figure 4, redundant water in the protruding bin caused by infiltration is redistributed to the bins to its left. In the next time step, the process begins with another infiltration followed by redistribution (Equation (4)). These two processes alternate until the simulation completes. The strength of the peak in that particular protruding bin with more infiltration will depend on the bin width Δθ due to its nonlinearity effect in Equation (4). The redistribution phase moves no surface water downward into the soil by assumption, but only rearranges the wetting fronts in different bins. Due to the properties of the decreasing front depths by the redistribution, an analysis of the instantaneous infiltration rates in the Talbot–Ogden model can be made.

#### 2.2.1. One Bin versus Two Bins

Figure 5 presents the model of one bin and its division into two bins. The rectangle *ACHF* represents *bin*1. Its wetting front is *FH*, its width is Δθ1, and its depth is *z*1. Similarly, the two bins in Figure 5 are *bin*2 (*ABJI*) and *bin*3 (*BCED*) with wetting fronts *IJ* and *DE*, respectively. Both their bin widths are equal to Δθ2, which means Δθ1 = 2Δθ2. Their front depths are *z*2 and *z*3. It is assumed that the entire water content for the two cases is the same. So, the water in *bin*1 equals that contained in the union of *bin*2 and *bin*3, that is *z*1 × Δθ1 = *z*2 × Δθ2 + *z*3 × Δθ2. These two cases can be compared. Let *VOnebin* and *VTwobins* denote the instantaneous infiltration rates for each case. *VOnebin* is calculated by

$$\begin{array}{rcl} V\_{\text{Conbin}} &= V\_{\text{bin}\_1} = \frac{dz\_1}{dt} \times \Delta \theta\_1\\ &= \frac{2 \times \Lambda \theta\_2 \times \mathcal{K}(\theta\_d)}{\theta\_d - \theta\_1} \Big(\frac{\psi(\theta\_d)}{z\_1} + 1\Big) \end{array} \tag{6}$$

**Figure 5.** One bin versus two bins.

*VTwobins* is calculated by

$$\begin{array}{rcl} V\_{T \text{uvbits}} &= V\_{\text{lin}\_2} + V\_{\text{bin}\_3} = \left( \frac{dz\_2}{dt} + \frac{dz\_3}{dt} \right) \times \Delta\theta\_2\\ &= \frac{\Delta\theta\_2 \mathcal{K}(\theta\_d)}{\theta\_d - \theta\_i} \left( \frac{\psi(\theta\_d)}{z\_2} + \frac{\psi(\theta\_d)}{z\_3} + 2 \right) \end{array} \tag{7}$$

Based on the capillary pressure of the soil, a relationship among depths exists: *z*2 > *z*1 > *z*3, since left bins always have deeper wetting fronts than the bins to their right. Moreover, if *z*2 = *z*3, then splitting one bin into two bins reverts back to the one bin case. Now, the difference between the instantaneous infiltration rates can be computed by subtracting (6) from (7). Thus,

$$\begin{array}{rcl} V\_{\text{Tuvbits}} - V\_{\text{Conbin}} &=& V\_{\text{bin}\_2} + V\_{\text{bin}\_3} - V\_{\text{bin}\_1} \\ &=& \frac{\Lambda \vartheta\_2 \times K(\vartheta\_4) \times \psi(\vartheta\_4)}{\vartheta\_4 - \vartheta\_i} \bigg( \frac{\left(z\_2 - z\_3\right)^2}{z\_2 z\_3 \left(z\_2 + z\_3\right)} \bigg) \\ & \ge 0 \text{ (with equality if } z\_2 = z\_3\text{)}. \end{array} \tag{8}$$

Equation (8) uses *z*1 × Δθ1 = *z*2 × Δθ2 + *z*3 × Δθ2. Therefore, when one bin is split into two bins and all possible pores are saturated, additional water infiltrates into the soil faster. Its increment can be bounded by (8). Doubling the number of bins leads to a more continuous flow among bins. Hence, more bins mean a higher infiltration rate. One bin is increased to *n* bins to prove this assumption.

#### 2.2.2. One Bin versus *n* Bins

Similarly, *n* bins can be obtained by equally discretizing the moisture content domain of the one bin case into *n* pieces. All other parameters remain constant. In Figure 6, the rectangle *ABCD* represents one bin *binX*. Its depth and bin width are *zX* and Δθ*X*, respectively. Then *binX* is split into *n* bins, marked by *bin*1, *bin*2, ... , *binn*−1, and *binn*. All these new bins have the same width Δθ. Their depths from left to right are *z*1, *z*2, ... , *zn*−1 and *zn*. It is similar to the previous example to have *zX* × Δθ*X* = (*<sup>z</sup>*1 + *z*2 + ... + *zn*) × Δθ, *n*Δθ = Δθ*X* = θ*d* − θ*i*, and *z*1 ≥ *z*2 ≥ ... ≥ *zn*. By geometry there exists an index *l* such that *z*1 ≥ *z*2 ≥ ... ≥ *zl* ≥ *zX* ≥ *zl*+1 ≥ ... ≥ *zn*. The infiltration rate of *binX* is

$$\begin{array}{rcl} V\_{\text{Conbin}} &= V\_{\text{bin}\chi} = \frac{dz\chi}{dt} \times \Delta\theta\_X\\ &= \frac{\Lambda\theta\_X \mathcal{K}(\theta\_d)}{\theta\_d - \theta\_i} \Big(\frac{\psi(\theta\_d)}{z\_\chi} + 1\Big) \end{array} \tag{9}$$

**Figure 6.** One bin versus *n* bins.

For the *n-bin* case, the summed infiltration rate is

$$\begin{array}{rcl} V\_{\text{nbins}} &=& \sum\_{j=1}^{n} V\_{\text{bin}\_j} = \sum\_{j=1}^{n} \frac{dz\_j}{d\theta} \times \Delta\theta \\ &=& \sum\_{j=1}^{n} \frac{\Delta\theta K(\theta\_d)}{\theta\_d - \theta\_i} \left(\frac{\psi(\theta\_d)}{z\_j} + 1\right) \end{array} \tag{10}$$

By Equation (9) and Equation (10), the infiltration difference between these two cases is

$$\begin{split} V\_{\text{nibins}} - V\_{\text{Conbin}} &= \left( \sum\_{j=1}^{n} V\_{\text{bin}\_j} \right) - V\_{\text{bin}\_X} \\ &= \frac{\Delta \theta \times K(\theta\_d) \times \psi(\theta\_d)}{\theta\_d - \theta\_i} \left( \sum\_{j=1}^{n} \frac{1}{z\_j} - \frac{n}{z\_X} \right) \\ &= \frac{\Delta \theta \times K(\theta\_d) \times \psi(\theta\_d)}{\theta\_d - \theta\_i} \left( \sum\_{j=1}^{n} \frac{1}{z\_j} - \frac{n^2}{\sum\_{j=1}^{n} z\_j} \right) \ge 0 \end{split} \tag{11}$$

Note that the last step of Equation (11) uses

$$\frac{\sum\_{j=1}^{n} z\_j}{n} \ge \left(\prod\_{j=1}^{n} z\_j\right)^{\frac{1}{n}} \ge \frac{n}{\sum\_{j=1}^{n} \frac{1}{z\_j}}.\tag{12}$$

Therefore from Equation (11), which means *Vnbins* − *VOnebin* ≥ 0, more bins in the model result in a greater infiltration rate.

Now, *n*1 bins is increased to *n*2 bins, where *n*2 and *n*1 are integers and *n*2 > *n*1. This case resembles the extension from one bin to *n*2/*<sup>n</sup>*1 bins. The difference in infiltration rates between these *n*1 and *n*2 bins is bounded by *Vn*2*bins* − *Vn*1*bins* ≈ *<sup>V</sup>n*2/*<sup>n</sup>*1*bins* − *VOnebin*. By Equation (11), we can conclude that if a soil texture can be fitted by a finer discretization, its overall conductivity becomes higher in the Talbot–Ogden model. We found the infiltration rate using asymptotic analysis.

#### *2.3. Asymptotic Analyses and Its Physical Meaning*

Asymptotic analysis is used when the number of bins increases toward infinity. Equation (11) and the assumption *n*Δθ = θ*d* − θ*i* are used first:

$$V\_{n\text{bins}} - V\_{\text{Conclev}} = K(\theta\_d) \times \psi(\theta\_d) \left( \frac{1}{n} \sum\_{j=1}^{n} \frac{1}{z\_j} - \frac{n}{\sum\_{j=1}^{n} z\_j} \right) \tag{13}$$

Before asymptotic analysis, the wetting front curve is loosely assumed to approximate a straight slanted line (Figure 7) in the θ-*<sup>z</sup>* domain. This assumption can also be demonstrated by the water content profiles in the solutions of the Richards Equation using Hydrus-1D [27] under many circumstances. Although it is only an approximation, all the wetting fronts in this model decrease from *z*1 to *zn* (Figure 6) after redistribution in every time step. When *n* → <sup>∞</sup>, the effect of the number of bins on the instantaneous infiltration rate in the Talbot–Ogden model is:

$$\begin{split} &= \lim\_{n \to \infty} (V\_{n\text{bins}} - V\_{\text{Conbin}}) \\ &= \lim\_{n \to \infty} \left[ K(\theta\_d) \times \psi(\theta\_d) \times \left( \frac{(z\_1 - z\_n) \sum\_{j=1}^n \frac{1}{z\_j}}{n(z\_1 - z\_n)} - \frac{z\_1 - z\_n}{\frac{z\_1 - z\_n}{n} \sum\_{j=1}^n z\_j} \right) \right] \\ &= \lim\_{n \to \infty} K(\theta\_d) \times \psi(\theta\_d) \times \begin{pmatrix} \frac{1}{z\_1 - z\_n} \int\_{z\_n}^{z\_1} \frac{1}{z} dz - \frac{z\_1 - z\_n}{\int\_{z\_1}^{z\_1} z dz} \\ \frac{1}{z\_1 - z\_n} \int\_{z\_n}^{z\_1} -\frac{1}{z\_1 + z\_n} \end{pmatrix} \\ &= \lim\_{n \to \infty} K(\theta\_d) \times \psi(\theta\_d) \times \left( \frac{1}{z\_1 - z\_n} \ln \frac{z\_1}{z\_n} - \frac{2}{z\_1 + z\_n} \right) . \end{split} \tag{14}$$

**Figure 7.** Approximately linearly decaying wetting fronts.

Equation (14) is obtained by the fundamental law of integration, which is the Newton–Leibniz rule. Summing the reciprocals and then taking the average results in a positive number. Therefore, in the integration the lower bound is 1/*zn*, and the upper bound is 1/*<sup>z</sup>*1.

Moreover, *z*1 corresponds to the deepest saturated bin, whereas *zn* corresponds to the shallowest one. For convenience let *zi* and *zd* denote these two depths, respectively, because of the moisture content they represent. Thus, when *n* → <sup>∞</sup>,

$$\begin{split} \|V\_{n\text{bins}} - V\_{\text{Oucbin}}\|\_{L\_1} &= \lim\_{n \to \infty} (V\_{n\text{bins}} - V\_{\text{Oucbin}}) \\ &= K(\theta\_d) \times \psi(\theta\_d) \times \left(\frac{1}{z\_i - z\_d} \ln \frac{z\_i}{z\_d} - \frac{2}{z\_i + z\_d}\right) \end{split} \tag{15}$$

What Equation (15) tells us is that the variation in infiltration rate depends only on the depths corresponding to the initial porosity θ*i* and the moisture content θ*d* of the rightmost saturated bin. These two special bins are actually the deepest and the shallowest ones, respectively. An interpretation is that the whole wetting front is pushed downward by water in the largest saturated porosity θ*d* and is prevented from progressing by water in the initial moisture content θ*i*. Therefore, the advancement of the wetting front is a compromise between these two bins. Interestingly, the unit of length does not count in the expression, and the quantity 1 *zi*−*zd ln zi zd* − 2 *zi*+*zd* is dimensional and is meaningful in physics. However, its value needs some further research, which is crucial for this model. The following practical case is presented: suppose the two wetting fronts satisfy

$$\exists \delta > 0 \; and \; \hat{\mathcal{N}} \in \mathbb{Z}^+ \; such \; that \; \; z\_d(t) + \hat{\mathcal{N}}\delta \ge z\_i(t) \ge z\_d(t) + \delta, \; for \; all \; t.$$

If *lim <sup>t</sup>*→∞*zd*(*t*) → <sup>∞</sup>, then the following limits hold:

$$\lim\_{z\_d \to \infty} \left( \frac{1}{z\_i - z\_d} \ln \frac{z\_i}{z\_d} - \frac{2}{z\_i + z\_d} \right) \le \lim\_{z\_d \to \infty} \left( \frac{1}{\delta} \ln \left( 1 + \frac{\hat{\mathcal{N}} \delta}{z\_d} \right) - \frac{2}{2z\_d + \hat{\mathcal{N}} \delta} \right) = 0 \tag{16}$$

and

$$\lim\_{z\_d \to \infty} \left( \frac{1}{z\_i - z\_d} \ln \frac{z\_i}{z\_d} - \frac{2}{z\_i + z\_d} \right) \ge \lim\_{z\_d \to \infty} \left( \frac{1}{\text{N}\delta} \ln \left( 1 + \frac{\delta}{z\_d} \right) - \frac{2}{2z\_d + \delta} \right) = 0. \tag{17}$$

Therefore, it follows that

$$\lim\_{z\_d \to \infty} \left( \frac{1}{z\_i - z\_d} \ln \frac{z\_i}{z\_d} - \frac{2}{z\_i + z\_d} \right) = 0. \tag{18}$$

Equation (18) means that if the wetting fronts corresponding to the initial moisture content θ*i* and the e ffective porosity θ*e* are within some distance from each other in depth, the infiltration in the Talbot–Ogden model will eventually become a steady state flow. However, a transformation will make this discussion easier.

Without loss of generality, let *zi* = *rzd* with *r* > 1, then a function *D*(*r*) characterizing the difference between *Vnbins* and *VOnebin*, is defined by *D*(*r*) = *ln*(*r*)/(*r* − 1) − 2/(*r* + 1). It follows that

$$\|V\_{\rm nbits} - V\_{\rm Oncbin}\|\_{L\_1} = \frac{K(\theta\_d) \times \psi(\theta\_d)}{z\_d} \times D(r). \tag{19}$$

The potential maximum or minimum of *D*(*r*) can be attained by taking its derivative

$$\frac{dD(r)}{dr} = -\frac{\ln(r)}{(r+1)^2} + \frac{1}{r(r-1)} + \frac{2}{(r+1)^2} \tag{20}$$

which has a zero point at *r* ≈ 8.16. There are also two limits for two extreme cases:

$$\lim\_{r \to 1^{+}} D(r) \quad = \lim\_{r \to 1^{+}} \frac{(r+1)\ln(r) - 2(r-1)}{r^2 - 1}$$

$$= 0\_{\epsilon}$$

and

$$\lim\_{r \to +\infty} D(r) \quad = \lim\_{r \to +\infty} \left[ \frac{1}{r-1} \ln(r) - \frac{2}{r+1} \right] \tag{22}$$

Equation (21) corresponds to the Green–Ampt model. Equation (22) means that if the deepest wetting front is too far away from the shallowest one, then the infiltration rate also tends to the one bin case (Figure 8).

**Figure 8.** Ratio r of the deepest wetting front to the shallowest one and its influence on *D*(*r*).

The curve in Figure 8 very clearly shows the change in *D*(*r*). With the maximum value of *D*(*r*), it is able to bound *Vnbins* − *VOnebin<sup>L</sup>*1 . Note that

$$\begin{aligned} \|V\_{n\text{bias}} - V\_{\text{Conbin}}\|\_{L\_1} &= \frac{K(\theta\_d) \times \psi(\theta\_d)}{z\_d} \times D(r) \\ &\leq \frac{0.0748 \times k\_\* \times \psi\_b \times \left(\frac{\theta\_d - \theta\_i}{\theta\_b - \theta\_i}\right)^{3 + \frac{1}{\lambda}}}{z\_d} (by \text{ Brooks and Corry}) \end{aligned} \tag{23}$$

where λ denotes the pore size distribution index [14]. If θ*d* = θ*<sup>e</sup>*, then the hydraulic models used in inequality (23) are unimportant because *K*(θ*d*) = *Ks and* ψ(θ*d*) = ψ*b* always hold. In this situation, hence,

$$\|V\_{n\text{bins}} - V\_{\text{Cuechin}}\|\_{L\_1} \le 0.0748 \times \frac{K\_s \times \psi\_b}{z\_d},\tag{24}$$

the upper bound of the infiltration rate error can be calculated. Note that if a soil texture is chosen, then *zd* can grow as a function of time and the rainfall rate, therefore making this upper bound decrease with respect to time.

## **3. Numerical Experiments**

In this section, the flux variation upper bound is calculated from the physical parameters of a variety of soil textures. Numerical tests are then carried out to simulate infiltration and verify our analysis in the previous section.

If *zd* is fixed in Equation (24), what is important is the product of saturated hydraulic conductivity *Ks* and the bubbling pressure ψ*<sup>b</sup>*. Therefore, 0.0748·*Ks*ψ*<sup>b</sup>* gives a bound for discrepancy of infiltration rates. In Table 1 [23], *Ks*ψ*b* decreases from coarser to finer soils. The one exception is between clay loam and silty clay loam. In general, the upper bound of *Vnbins* − *VOnebin<sup>L</sup>*1 becomes larger with coarser soils. Significantly, this finding means that if the soil is finer, then the infiltration will not be distinguishable based on a change in the number of bins. However, in coarse textured soil systems, the number of bins is more important than in finer ones. This is the reason why more bins may be required to test the coarser soil: the outcome varies in a wider range.


**Table 1.** Parameters of different soil systems.

From the last column of Table 1, if *zd* is fixed and not too deep nor too shallow, then the value 0.0748·*Ks*ψ*<sup>b</sup>* is in fact a certain bound for the variation in infiltration rate. This condition means that if the rainfall rate is around that value or of the same scale for a specific soil system, then the choice of a proper number of bins to make the simulation realistic is needed. Hence, the choice of the number of bins should be determined at least by both the soil and rainfall rates. Note that when *zd* is very large in inequality (24), the bound becomes small, which means the Talbot–Ogden model acts as the Green–Ampt model or the Richards model for the steady state flow.

Using the Talbot–Ogden method it is possible to simulate this infiltration process by choosing some parameters, including an appropriate number of bins. If the proper number of bins is relatively large, then the modeled soil lets water quickly pass through it, which indicates that an underground reservoir recharges quicker. On the other hand, if fewer bins are needed, then the model indicates that this soil has a good capacity to retain water, which is beneficial for growing plants and for agriculture. In general, the number of bins may indicate the pore size categories that can be modeled in the soil for infiltration.

The comparison of instantaneous infiltration rates computed by Talbot–Ogden method and Hydrus-1D is tested. In this work, the main focus is the influence of the number of bins on the infiltration rates in the Talbot–Ogden method. The soil parameters used are listed in Table 2 [14]. The time step Δ*t* is set as 10 *s* for all the simulations to guarantee the comparison consistency. Figure 9 shows the infiltration rate curves corresponding to different numbers of bins for three soil textures. These three types are sand, sandy clay, and silt loam. Different rainfall rates are used because each soil

type can generate comparable data only under the proper rainfall rate. One continuous rainfall is one pulse. There are two pulses, and each lasts for 1.5 h. The numbers of bins used in testing are 25, 125, and 250.

**Figure 9.** Comparison of infiltration rates in (**a**) sand, (**b**) silt loam, and (**c**) sandy clay using different numbers of bins.


**Table 2.** Soil textures and parameters used in model evaluation.

With an increase in the numbers of bins, the infiltration rate in any soil system rises. During the second precipitation, this phenomenon is even more obvious than in the first one. The reason is explained by the analysis given in Section 2. At the beginning of the second rainfall there are already some dry pores with larger pores that correspond to θ near θ*<sup>e</sup>*, which is to the right of the moisture content domain. The depth *zd* is very small at this moment so that the change in infiltration rate is more sensitive to the perturbation of the numbers of bins than that in the first pulse. Mathematically, if *zd* < 1 cm (the unit is consistent with that in Table 1), the upper bounds given by inequality (24) and Table 1 are relatively easier to approach during the second pulse, thus verifying the analysis in Section 2.

The quadratic or root mean squared (*RMS*) difference in the instantaneous infiltration rates calculated at every time step with 25 bins and 125 or 250 bins, are defined as

$$RMS\_1 = \sqrt{\frac{\sum\_{time \ step=1}^{N\_1} \left(f\_{25 \ binis}^{time \ step} - f\_{125 \ binis}^{time \ step}\right)^2}{N\_1}} \tag{25}$$

and

$$RMS\_2 = \sqrt{\frac{\sum\_{\substack{t=2 \\ t \text{ times}}}^{N\_2} \left(f\_{25 \text{ bits}}^{\text{time step}} - f\_{250 \text{ bits}}^{\text{time step}}\right)^2}{N\_2}},\tag{26}$$

where *f* represents the instantaneous infiltration rates as computed for different numbers of bins, respectively, and *N* is the number of *f* values in the interval for which *RMS* is calculated. In Table 3, the columns of *RMS* values and their comparisons show that the effect on the infiltration rate becomes smaller with an increase in the number of bins. In the fourth column of Table 3, the decreasing *RMS*1/*RMS*2 values indicate that the coarser the soil texture is, the more sensitive the infiltration rate is to the change in the number of bins, especially at the beginning of this change. This phenomenon can be attributed to the relatively large range of infiltration change for coarser soils as shown in Table 1. Note that the root mean square values in Table 3 also satisfy *RMS*1,2 < 0.0748·*Ks*ψ*<sup>b</sup>*.

**Table 3.** Influence of the number of bins on the infiltration in the Talbot–Ogden model.


Root mean sqaure1 = The difference of infiltration rates between 25 and 125 bins; Root mean square2 = The difference of infiltration rates between 25 and 250 bins.

Figure 10a shows how infiltration behaves when the number of bins increases rapidly. There are three arrays of data for comparison: one with 100 bins, another with 500 bins, and the rest with bins.The rainfall rate is set to 2 cm/h, and two pulses are selected with each lasting 1.5 h. Figure 10b shows similar results but with a coarser soil system silt loam. Let *N* denote the number of time steps where

infiltration discrepancies appear corresponding to different numbers of bins, the following equation is used to calculate the influence of the rapid increase in the number of bins:

$$\text{Influence of large number of bins} (\mathbf{x}) = \sqrt{\frac{\sum\_{i=1}^{N} \left(\frac{V\_{i,100\text{ bin}} - V\_{i,x\text{ bin}}}{V\_{i,100\text{ bin}}}\right)^2}{N}} \times 100\%. \tag{27}$$

**Figure 10.** Infiltration in (**a**) sandy clay loam, and (**b**) silt loam with a large increase in the number of bins.

For the three cases in Figure 10a,b, the use of (27) yields Table 4. In Table 4, for sand clay loam, the difference in average infiltration rates using 100 bins and bins is

$$2.0\text{ cm}/h \times 8.56\% = 0.1712\text{ cm}/h < 0.0748 \cdot \text{K}\_{\delta} \psi\_{b} = 0.630\text{ cm}/h. \tag{28}$$

**Table 4.** Influence of large number of bins on the infiltration rate compared with 100 bins as the baseline.


Similarly, for silt loam, their di fference is

$$0.3.5\,\mathrm{cm}/\mathrm{h} \times 24.32\% = 0.8512\,\mathrm{cm}/\mathrm{h} < 0.0748 \cdot \mathrm{K}\_{\mathrm{s}}\psi\_{\mathrm{b}} = 1.058\,\mathrm{cm}/\mathrm{h}.\tag{29}$$

In Equations (28) and (29) show that if the number of bins increases by 100 times, then the increment in the infiltration rates is still bounded by the data in Table 1. In fact, the increment of infiltration rate is more significant from 100 bins to 500 bins. These results show that infiltration in the finer soil system is less sensitive to the rapid increase in the number of bins than the coarser one, which is consistent with the analysis in Equation (24) and Table 1.
