**Carlos Fuentes <sup>1</sup> and Carlos Chávez 2,\***


Received: 7 September 2020; Accepted: 26 September 2020; Published: 27 September 2020

**Abstract:** The aim of this study is the deduction of an analytic representation of the optimal irrigation flow depending on the border length, hydrodynamic properties, and soil moisture constants, with high values of the coefficient of uniformity. In order not to be limited to the simplified models, the linear relationship of the numerical simulation with the hydrodynamic model, formed by the coupled equations of Barré de Saint-Venant and Richards, was established. Sample records for 10 soil types of contrasting texture were used and were applied to three water depths. On the other hand, the analytical representation of the linear relationship using the Parlange theory of infiltration proposed for integrating the differential equation of one-dimensional vertical infiltration was established. The obtained formula for calculating the optimal unitary discharge is a function of the border strip length, the net depth, the characteristic infiltration parameters (capillary forces, sorptivity, and gravitational forces), the saturated hydraulic conductivity, and a shape parameter of the hydrodynamic characteristics. The good accordance between the numerical and analytical results allows us to recommend the formula for the design of gravity irrigation.

**Keywords:** Saint-Venant equations; Richards' equation; Parlange equations; optimal irrigation flow; soil parameters; analytical representation

#### **1. Introduction**

Gravity irrigation is the water supply at the head of a channel or inclined ditch built on a plot, as a border or a furrow, in order to take advantage of the gravitational field to provide the necessary amount of water for optimal development of cultivated plants. In continuous gravity irrigation, three phases are distinguished in the surface water movement [1,2]. The first begins when water flow is provided on the dry border until the water wave reaches the end part of the same; it is known as the advance phase. The second starts from the arrival of the wave at the end of the border until the water supply is cut off, known as the storage phase. Finally, the third phase, known as the recession, is composed of two sub-phases, one is the vertical recession starting from cutting off the water supply until the depth at the head of the border disappears, and the other is the horizontal recession starting from the disappearance of the depth at the head and ends when the depth at the end of the border disappears.

According to the principles used in the modeling, studies reported in literature can be grouped in the context of four approaches [2]: (1) The modeling of surface and underground movements is addressed in a full empirical way [3]; (2) the surface movement is modeled with the Barré de Saint-Venant equations and their simplifications, and the underground movement with empirical equations as those of Kostiakov and Mezencev [4–11]; (3) surface movement is modeled by Saint-Venant simplified equations (kinematic wave, diffusion wave or null inertia and hydrological model) and in the modeling of the underground movement there is the possibility of using rational equations [12–23]; (4) the surface movement is modeled with Barré de Saint-Venant equations and the underground movement with Richards' equation [24] in [2,25–28].

In the study to model the border irrigation developed by Schmitz et al. [29], the complete equations of Barre de Saint-Venant are solved with the method of characteristics, and for the underground movement, the analytical solution for infiltration obtained by Parlange et al. is used [30]. Other studies use the Green and Ampt infiltration equation [31], which is a special case of Richard's two-dimensional equation [32], and other group coupled Saint-Venant and Richards equations in border irrigation. The first ones are resolved with a Lagrangian method and the second one with the finite element method [25–28]. In Saucedo et al. [27], using the full hydrodynamic model, the optimal irrigation flow is obtained using numerical methods, with Saint-Venant equations coupled internally with Richards' equation.

The aim of this study is the deduction of an analytic representation of the optimal irrigation flow depending on the border length, hydrodynamic properties, and soil moisture constants, with high values of the coefficient of uniformity.

#### **2. Materials and Methods**

#### *2.1. Water Surface Flow*

The continuity and amount of movement equations in a border—considering that the effects of the borders were negligible and that water was the shallow or hydraulic hypothesis—were known as equations of Barre de Saint-Venant for border irrigation and written as follows:

$$\frac{\partial \mathbf{J}}{\partial \mathbf{t}} + \frac{\partial \mathbf{J}}{\partial \mathbf{x}} + \frac{\partial \mathbf{J}}{\partial \mathbf{t}} = \mathbf{0} \tag{1}$$

$$\frac{1}{\hbar} \frac{\partial \mathbf{q}}{\partial \mathbf{t}} + \frac{2 \mathbf{q}}{\hbar^2} \frac{\partial \mathbf{q}}{\partial \mathbf{x}} + \left(\mathbf{g} - \frac{\mathbf{q}^2}{\hbar^3}\right) \frac{\partial \mathbf{h}}{\partial \mathbf{x}} + \mathbf{g} (\mathbf{J} - \mathbf{J}\_o) + \lambda \frac{\mathbf{q}}{\hbar^2} \frac{\partial \mathbf{J}}{\partial \mathbf{t}} = \mathbf{0} \tag{2}$$

where q(x, t) = U(x, t) h(x, t) is the flow per width unit of border, x is the spatial coordinate in the main direction of movement of water in the border; t is time; U is the mean velocity; h the water depth; J<sup>o</sup> is the topographic slope of the border; J is the friction slope; V<sup>I</sup> = ∂I/∂t is the infiltration flow, that is, the volume of water infiltrated in the time unit per width unit and per length unit of the border, I is the infiltrated depth; g is the gravitational acceleration; the dimensionless parameter λ = UIX/U, with UIX the projection in the movement direction of the output speed of the water body due to infiltration.

The system of equations of Barré de Saint-Venant was not closed since the evolution in time and space of infiltrated depth and friction slope, were unknown. The first was provided by the Richards equation [24] and the second by a law of resistance to the flow that related the friction slope with the hydraulic variables q and h, which were discussed below [2].

#### *2.2. Water Flow in the Soil*

If the hypothesis that the irrigation was carried out in flat parallels to the development of the border was accepted, then it is possible to use the two-dimensional form of Richards' equation [24], which results from combining the continuity equation and Darcy's law [33]:

$$\mathbf{C}(\boldsymbol{\psi})\frac{\partial\boldsymbol{\psi}}{\partial\mathbf{t}} = \frac{\partial}{\partial\mathbf{x}}\bigg[\mathbf{K}(\boldsymbol{\psi})\frac{\partial\boldsymbol{\psi}}{\partial\mathbf{x}}\bigg] + \frac{\partial}{\partial\mathbf{z}}\bigg[\mathbf{K}(\boldsymbol{\psi})\Big(\frac{\partial\boldsymbol{\psi}}{\partial\mathbf{z}} - 1\Big)\bigg] \tag{3}$$

where ψ is the potential for water pressure in the soil expressed as the height of an equivalent water column (positive in the saturated zone and negative in the unsaturated zone of soil); C(ψ) = dθ/dψ is called the specific capacity of soil moisture; θ = θ(ψ) is the water volume per volume unit of soil or water volume content and is a function of ψ, known as the moisture characteristic curve or water retention

curve; K = K(ψ) is the hydraulic conductivity, which in a partially saturated soil is a function of the potential of pressure; the gravitational potential is assimilated to the spatial coordinate z positively oriented downwards, x is a spatial coordinate and t is the time.

For the description of the water flow during an irrigation test, the hydrodynamic characterization of soil was necessary. As pointed out by Fuentes et al. [34] in experimental studies, it was more convenient to use the combination of the retention curve proposed by van Genuchten [35], considering the Burdine restriction [36], with the hydraulic conductivity curve proposed by Brooks and Corey [37], due to the fact that they satisfy the integral properties of infiltration and to the ease of identification of their parameters. The retention curve proposed by van Genuchten is written as:

$$\frac{\Theta(\boldsymbol{\psi}) - \Theta\_{\rm r}}{\Theta\_{\rm s} - \Theta\_{\rm r}} = \left[ 1 + \left( \frac{\boldsymbol{\psi}}{\boldsymbol{\psi}\_{\rm d}} \right)^{\rm n} \right]^{-\rm m} \tag{4}$$

where θ<sup>s</sup> is the volumetric water content at effective soil saturation, θ<sup>r</sup> is the volumetric content of residual water, ψ<sup>d</sup> is a characteristic value of water pressure in the soil, m and n are two parameters of empirical form related here by the Burdine restriction [36]: m = 1 − 2/n, with 0 < m < 1 and n > 2.

The hydraulic conductivity proposed by Brooks and Corey [37] is represented as:

$$\mathbf{K}(\boldsymbol{\Theta}) = \mathbf{K}\_{\mathbf{s}} \left( \frac{\boldsymbol{\Theta} - \boldsymbol{\Theta}\_{\mathbf{r}}}{\boldsymbol{\Theta}\_{\mathbf{s}} - \boldsymbol{\Theta}\_{\mathbf{r}}} \right)^{\eta} \tag{5}$$

where η is a parameter of positive form whose value can be estimated with η = 2s (2/mn + 1), being s a function of porosity (φ) defined implicitly by (1 − φ) <sup>s</sup> + φ2s = 1 [38].

#### *2.3. Hydraulic Resistance Law*

The phase of advance in gravity irrigation is represented by the following initial and boundary conditions in the Barré de Saint-Venant equations:

$$\mathbf{q}(\mathbf{x},0) = 0 \quad \text{and} \quad \mathbf{h}(\mathbf{x},0) = \mathbf{0} \tag{6}$$

$$\mathbf{q}(0,\mathbf{t}) = \mathbf{q}\_{0'} \quad \mathbf{q}(\mathbf{x}\_{\mathbf{t}}, \mathbf{t}) = 0 \quad \text{and} \quad \mathbf{h}(\mathbf{x}\_{\mathbf{t}}, \mathbf{t}) = \mathbf{0} \tag{7}$$

where x<sup>f</sup> (t) is the position of the wave front for the time t and q<sup>o</sup> the flow of supply at the entrance of the border.

An analysis of the singularity present in the advance phase in very short times has established that the law of hydraulic resistance, which makes compatible the coupling of Barré de Saint-Venant and Richards equations in this singularity, and it has the following structure [2]:

$$\mathbf{q} = \mathbf{k}\mathbf{v} \left(\frac{\mathbf{h}^3 \mathbf{g} \mathbf{J}}{\mathbf{v}^2}\right)^{\mathbf{d}} \tag{8}$$

where ν is the coefficient of kinematic viscosity, k is a dimensionless factor of friction, and d is an exponent such that 1/2 ≤ d ≤ 1 in a way that d = 1/2 corresponds to the Chézy turbulent regime and d = 1 to the Poiseuille depth regime.

#### **3. Results and Discussion**

#### *3.1. Numerical Relationship between the Optimal Flow and Length*

The solution of Barre de Saint-Venant and Richards equations to represent surface and subsurface movement, respectively, has been solved numerically based on a Eulerian–Lagrangian scheme that eliminates the traditional instabilities in short times and is available on Saucedo et al. [25].

#### 3.1.1. Optimal Flow in Border Irrigation

With the system of Barré de Saint-Venant and Richards equations, the irrigation design consists of determining the flow of optimal supply and irrigation time to achieve the greatest uniformity along the border, with high levels of application efficiency and irrigation requirement for an irrigation depth and soil hydrodynamic predetermined properties. The optimal flow should be determined for a border length, and its value should be updated in proportion to the new length, which was verified by Rendón et al. [39]. In fact, the optimal flow design follows a linear proportion with the border length that should be applied. The result is obtained using a model formed by the Lewis and Milne [40] equations to describe the water flow on the soil surface and by the Green and Ampt equation [31] to describe the water flow on the soil.

#### 3.1.2. Irrigation Efficiencies

In irrigation, it is essential to distinguish at least three related efficiencies: Application efficiency, irrigation requirement efficiency, and irrigation uniformity efficiency. Application efficiency (ηA) is defined as [1,23]:

$$\eta\_{\rm A} = \frac{\mathbf{V\_n}}{\mathbf{V\_b}} = \frac{\ell\_{\rm n}}{\ell\_{\rm b}} \tag{9}$$

where V<sup>n</sup> is the water volume required in the root zone of the crop or net volume and V<sup>b</sup> is the amount of applied irrigation water. The first is obtained w.ith the expression: V<sup>n</sup> = `nAr, where `<sup>n</sup> is the net irrigation depth, defined according to the crop irrigation requirements, and A<sup>r</sup> is the irrigated area considered. The second is obtained as V<sup>b</sup> = `bAr, where `<sup>b</sup> is the gross irrigation depth.

Irrigation requirement efficiency (ηR) is defined as [23]:

$$
\eta\_{\rm R} = \frac{\mathbf{V\_{d}}}{\mathbf{V\_{n}}} \tag{10}
$$

where V<sup>d</sup> is the available volume by the crop. This efficiency indicates how the water needs for the crop are met.

The ideal situation regarding uniformity occurs when all plants receive the same amount of water, a situation which is equivalent to applying the same water depth to the entire length of the border. To evaluate the uniformity in distribution of the infiltrated depth, the Christiansen uniformity coefficient is used. This coefficient (CUD) results from partitioning the length in N sections of size ∆x<sup>i</sup> , not necessarily equal, namely [23]:

$$\text{CU}\_{\text{D}}(\mathbf{t}) = 1 - \frac{1}{\bar{\mathbf{I}}(\mathbf{t})\mathbf{L}} \sum\_{i=1}^{N} \left| \mathbf{I}\_{\text{i}}(\mathbf{t}) - \mathbf{\tilde{I}}(\mathbf{t}) \right| \Delta \mathbf{x}\_{\text{i}} \ \mathbf{\tilde{I}}(\mathbf{t}) = \frac{1}{\mathbf{L}} \sum\_{i=1}^{N} \mathbf{I}\_{\text{i}}(\mathbf{t}) \Delta \mathbf{x}\_{\text{i}} \tag{11}$$

where I<sup>i</sup> is the infiltrated depth at any section i of the border strip or furrow, I is the average infiltrated depth, and N is the number of sections considered along the furrow or border strip.

Christiansen classic uniformity coefficient (CUC) results when sections are taken of the same size, L = N∆x.

#### 3.1.3. Optimal Flow-Length Relationship

The uniformity efficiency measured by the Christiansen uniformity coefficient can be obtained for different combinations of length and supply flow at the head of the border. Saucedo et al. [27] showed an example of four lengths of the border for the soil under study, where it was observed that the uniformity efficiency varied considerably with the irrigation flow.

For each border length, it was possible to determine the value of the supply flow that produced a maximum in the uniformity coefficient with high values of application efficiencies and irrigation requirements. When the supply flow was modified, application efficiencies and irrigation requirements did not vary significantly, not in the case of the uniformity efficiency, which varied substantially with the irrigation flow, i.e., application and irrigation requirement efficiencies can be considered that are not decision variables in defining the optimal flow and, therefore, the uniformity efficiency is that which allows defining the optimal irrigation flow. *Water* **2020**, *12*, x FOR PEER REVIEW 5 of 14

The numerical simulation of the irrigation with the system of Barre de Saint-Venant and Richards equations indicates that the relationship between optimal flow (q<sup>o</sup> ) and length (L) is approximately linear, for a soil type, topographic slope, friction coefficient, and irrigation depth, that is: The numerical simulation of the irrigation with the system of Barre de Saint-Venant and Richards equations indicates that the relationship between optimal flow (qo) and length (L) is approximately linear, for a soil type, topographic slope, friction coefficient, and irrigation depth, that is:

$$\mathbf{q}\_{\rm o} = \mathbf{q}\_{\rm u} \mathbf{L} \tag{12}$$

where q<sup>u</sup> has units of unitary flow per length unit. In addition, as q<sup>u</sup> is a constant, it follows that for the application of a specific irrigation depth, there is an irrigation time, unique and independent of the length, which allows obtaining a maximum value of the uniformity coefficient. where q<sup>u</sup> has units of unitary flow per length unit. In addition, as q<sup>u</sup> is a constant, it follows that for the application of a specific irrigation depth, there is an irrigation time, unique and independent of the length, which allows obtaining a maximum value of the uniformity coefficient.

o u

The relation (12) is illustrated by Saucedo et al. [27] in the loam soil of the experimental field of the Colegio de Postgraduados, Montecillo, State of Mexico, for water depths of 8, 10, and 12 cm, by making q<sup>u</sup> = αuKs, where α<sup>u</sup> is a dimensionless parameter (Figure 1). It is observed that there is monotony in the sense that the slope of the relationship between the border length and optimal flow decreases as the irrigation depth increases. The relation (12) is illustrated by Saucedo et al. [27] in the loam soil of the experimental field of the Colegio de Postgraduados, Montecillo, State of Mexico, for water depths of 8, 10, and 12 cm, by making q<sup>u</sup> = αuKs, where α<sup>u</sup> is a dimensionless parameter (Figure 1). It is observed that there is monotony in the sense that the slope of the relationship between the border length and optimal flow decreases as the irrigation depth increases.

**Figure 1.** Relationship between the border strip length and the optimum input discharge in the loam soil of Montecillo for three irrigation depths: 8, 10, and 12 cm. K<sup>s</sup> in m/s. **Figure 1.** Relationship between the border strip length and the optimum input discharge in the loam soil of Montecillo for three irrigation depths: 8, 10, and 12 cm. Ks in m/s.

#### 3.1.4. Irrigation Design Table 3.1.4. Irrigation Design Table

To generate the design table, it was necessary to obtain the relationships between optimal flow and border length for various soil types. The moisture content of residual water (θr) was assumed to be zero, according to Fuentes et al. [34]. The moisture content at saturation (θs) was assimilated to the total soil porosity (ϕ) determined, as the hydraulic conductivity at saturation (Ks), starting from the soil texture according to the relationships provided by Rawls and Brakensiek [41]. To generate the design table, it was necessary to obtain the relationships between optimal flow and border length for various soil types. The moisture content of residual water (θr) was assumed to be zero, according to Fuentes et al. [34]. The moisture content at saturation (θs) was assimilated to the total soil porosity (φ) determined, as the hydraulic conductivity at saturation (Ks), starting from the soil texture according to the relationships provided by Rawls and Brakensiek [41].

To estimate the value of the shape m parameter of the soil retention curve, one granulometric curve was reconstructed for each soil based on the percentages of sand, silt, and clay present in the To estimate the value of the shape m parameter of the soil retention curve, one granulometric curve was reconstructed for each soil based on the percentages of sand, silt, and clay present in the

triangle of textures [34]; we have followed the procedure suggested by Fuentes [38] to determine the values of m and η. The pressure scale (ψd) was determined from the suction of the wetting front of triangle of textures [34]; we have followed the procedure suggested by Fuentes [38] to determine the values of m and η. The pressure scale (ψd) was determined from the suction of the wetting front of Green and Ampt equation [31] according to the texture and porosity of soils [41], as this suction was identified with the Bouwer scale [42] defined by:

$$\lambda\_{\mathbf{c}} = \frac{1}{\mathbf{K}\_{\mathbf{s}} - \mathbf{K}\_{0}} \int\_{\boldsymbol{\Psi}\_{0}}^{\boldsymbol{0}} \mathbf{K}(\boldsymbol{\psi}) \mathbf{d}\boldsymbol{\psi} = \frac{1}{\mathbf{K}\_{\mathbf{s}} - \mathbf{K}\_{0}} \int\_{\boldsymbol{\Theta}\_{0}}^{\boldsymbol{\Theta}\_{\mathbf{s}}} \mathbf{D}(\boldsymbol{\theta}) \mathbf{d}\boldsymbol{\theta} \tag{13}$$

where K<sup>0</sup> = K(ψ0) is the hydraulic conductivity corresponding to the initial pressure ψ<sup>0</sup> = ψ(θ0) related to the initial moisture content θ0, D(θ) = K(θ) dψ/dθ is the hydraulic diffusivity.

The parameter ψ<sup>d</sup> is deduced by introducing the hydrodynamic characteristics defined by Equations (4) and (5) in Equation (13), considering the initial moisture content equal to the residual moisture content (ψ<sup>0</sup> → ∞), as follows:

$$
\lambda\_{\mathbf{c}} = \left| \psi\_{\mathbf{d}} \right| \frac{1}{\mathbf{n}} \mathbf{B} \left( \mathfrak{m} - \frac{1}{\mathbf{n}}, \frac{1}{\mathbf{n}} \right) \tag{14}
$$

where B(p,q) = Γ(p)Γ(q)/Γ(p + q) is the complete beta function, with p > 0 and q > 0, and Γ(x) the Euler complete gamma function.

The parameter values of the hydrodynamic characteristics are shown in Table 1 for different types of soil [27], where the residual moisture content is θ<sup>r</sup> = 0 cm<sup>3</sup> /cm<sup>3</sup> .


**Table 1.** Hydrodynamic characteristics of soils for irrigation design by borders.

The initial moisture content is considered as that which corresponds when the available moisture of each soil type has been consumed in a certain fraction. The available soil moisture is defined as the difference between the moisture contents at field capacity (θCC) and permanent wilting point (θPMP), whose values for each type of soil are estimated according to the soil texture triangle [41]. The initial moisture content is calculated as:

$$\Theta\_0 = \Theta\_{\text{PMP}} + \mathcal{F}\_{\text{ap}} (\theta\_{\text{CC}} - \theta\_{\text{PMP}}) \tag{15}$$

where Fap is the permissible depletion factor of the crop. The average value of 0.5 has been assumed. The values of initial moisture content are reported in Table 2.


**Table 2.** Moisture constants.

As for the parameters of the law of resistance defined by Equation (8), the values were taken from a loam soil border in the Montecillo experimental field. Considering the Reynolds number, the regime is depth, d = 1, the value k = 1/54 is obtained thus that the advance curve provided by the numerical solution describes the advance curve observed in an irrigation test; the coefficient of kinematic viscosity is taken as ν = 10−<sup>6</sup> m<sup>2</sup> /s. The longitudinal topographical slope of the border is of J<sup>0</sup> = 0.002 m/m, value that is used to simulate irrigation in other borders with different soil types. With the hydrodynamic characterization of soils and θ0, the constant involved in the relationship between the border length and the optimal flow for a given irrigation depth is calculated.

The value of the constant is expressed in terms of flow per unit area, i.e., per unit width, and per unit length of border, results are shown in Table 3. The same table shows the irrigation time (τb) obtained at the moment the flow supply is cut off when the volume of irrigation per width unit is already stored both on the surface as well as inside the soil. In Rendón et al. [39], a table of similar design to Table 3 is shown containing some inconsistencies of monotony between the relationship that the variables optimal flow, irrigation time, and applied depth since the used model has difficulties in reproducing gravity irrigation in relatively long times. The coupling of Saint-Venant and Richards equations allows obtaining results that keep the monotony in the design variables, as shown in Table 3.


**Table 3.** Table of the border irrigation design: flow in l/s/m<sup>2</sup> for the optimal application of the irrigation depth.

#### *3.2. Analytical Representation of the Relationship between the Optimal Flow and Length*

Equation (12) structure is deduced considering that the net water volume per width unit of the border is equal to the product of the border length (L) by the net irrigation depth (`n) and is also equal to the product of supply unitary flow (qo) by the time necessary to infiltrate the net depth (τn): qoτ<sup>n</sup> = L`n. The relationship is also attested by involving the irrigation time to obtain the volume per width unit provided by the gross depth: q<sup>o</sup> τ<sup>b</sup> = L`b. Comparing both results to Equation (12), we have:

$$\mathbf{q}\_{\mathbf{u}} = \frac{\ell\_{\mathbf{n}}}{\pi\_{\mathbf{n}}} = \frac{\ell\_{\mathbf{b}}}{\pi\_{\mathbf{b}}} \tag{16}$$

It should be noted that this relationship involves, considering Equation (9), the following expression for application efficiency: It should be noted that this relationship involves, considering Equation (9), the following expression for application efficiency:

$$
\eta\_{\rm A} = \frac{\ell\_{\rm n}}{\ell\_{\rm b}} = \frac{\tau\_{\rm n}}{\tau\_{\rm b}} \tag{17}
$$

From the continuity equation, it can be shown that the unitary flow of minimal irrigation thus that the water wave arrives at the end of the channel is given by q<sup>m</sup> = K<sup>s</sup> L; then it follows that the optimal flow must meet the inequality q<sup>o</sup> ≥ qm. If q<sup>o</sup> = α<sup>u</sup> q<sup>m</sup> is written, α<sup>u</sup> is a dimensionless parameter that must satisfy α<sup>u</sup> ≥ 1. Equations (12) and (16) should be written as follows: From the continuity equation, it can be shown that the unitary flow of minimal irrigation thus that the water wave arrives at the end of the channel is given by q<sup>m</sup> = K<sup>s</sup> L; then it follows that the optimal flow must meet the inequality q<sup>o</sup> ≥ qm. If q<sup>o</sup> = α<sup>u</sup> q<sup>m</sup> is written, α<sup>u</sup> is a dimensionless parameter that must satisfy α<sup>u</sup> ≥ 1. Equations (12) and (16) should be written as follows:

$$\mathbf{q}\_{\rm o} = \alpha\_{\rm u} \mathbf{K}\_{\rm s} \mathbf{L}, \; \alpha\_{\rm u} = \frac{\mathbf{q}\_{\rm u}}{\mathbf{K}\_{\rm s}} = \frac{\ell\_{\rm n}}{\mathbf{K}\_{\rm s} \tau\_{\rm n}} \tag{18}$$

in which the dependence of α<sup>u</sup> should be investigated regarding the irrigation depth and soil properties. in which the dependence of α<sup>u</sup> should be investigated regarding the irrigation depth and soil properties.

The extreme behavior of α<sup>u</sup> is deduced from the extreme behavior of the infiltrated depth. In very short times I = S √ t [43], where S is sorptivity, and, therefore, τ<sup>n</sup> = ` 2 n /S 2 , that is α<sup>u</sup> = S 2 /(K<sup>s</sup> `n). In long times I~I<sup>o</sup> + K<sup>s</sup> t, where I<sup>o</sup> is the ordinate at the origin depending on S and Ks, and on the Green and Ampt model on the time logarithm, and, therefore, αu~`n/(`<sup>n</sup> − Io). From the above, it follows that the limits: The extreme behavior of α<sup>u</sup> is deduced from the extreme behavior of the infiltrated depth. In very short times I = S√t [43], where S is sorptivity, and, therefore, τ<sup>n</sup> = ℓ<sup>n</sup> 2 /S<sup>2</sup> , that is α<sup>u</sup> = S2/(K<sup>s</sup> ℓn). In long times I~I<sup>o</sup> + K<sup>s</sup> t, where I<sup>o</sup> is the ordinate at the origin depending on S and Ks, and on the Green and Ampt model on the time logarithm, and, therefore, αu~ℓn/(ℓ<sup>n</sup> − Io). From the above, it follows that the limits:

$$\lim\_{\ell\_{\mathbf{u}} \to 0} \mathfrak{a}\_{\mathbf{u}} = \infty, \lim\_{\ell\_{\mathbf{u}} \to \infty} \mathfrak{a}\_{\mathbf{u}} = 1 \tag{19}$$

must be satisfied by the general function α<sup>u</sup> (`n).

Irrigation time (τb) shown in Table 3 corresponds to the gross irrigation depth and is greater than the infiltration time corresponding to the net depth calculated from Equation (16): τ<sup>n</sup> = `n/qu. Figure 2 shows the relationship between the two times; one has <sup>τ</sup><sup>n</sup> <sup>≈</sup> 0.83 <sup>τ</sup><sup>b</sup> with r<sup>2</sup> <sup>=</sup> 0.9995, which indicates that, according to Equation (17), the average application efficiency with the optimal flow is η<sup>A</sup> ≈ 83% for analyzed soils. must be satisfied by the general function α<sup>u</sup> (ℓn). Irrigation time (τb) shown in Table 3 corresponds to the gross irrigation depth and is greater than the infiltration time corresponding to the net depth calculated from Equation (16): τ<sup>n</sup> = ℓn/qu. Figure 2 shows the relationship between the two times; one has τ<sup>n</sup> 0.83 τ<sup>b</sup> with r <sup>2</sup> = 0.9995, which indicates that, according to Equation (17), the average application efficiency with the optimal flow is η<sup>A</sup> 83% for analyzed soils.

**Figure 2.** Relationship between the irrigation time of the gross depth (ℓb) and the infiltration time of the net depth (ℓn). **Figure 2.** Relationship between the irrigation time of the gross depth (`b) and the infiltration time of the net depth (`n).

The design table can be represented algebraically if the infiltration function is provided

*Water* **2020**, *12*, 2710

The design table can be represented algebraically if the infiltration function is provided analytically. This is obtained from the Parlange et al. model [44] deduced from the Richards equation, assuming that the hydraulic diffusivity tends to behave like a Dirac density and form a relationship between the hydraulic diffusivity and hydraulic conductivity. The model with the effect of water depth on the soil surface is as follows [30]:

$$\mathbf{I(t)} - \mathbf{K\_0t} = \frac{\mathbf{K\_s}\mathbf{h}\Delta\theta}{\mathbf{q\_s(t)} - \mathbf{K\_s}} + \frac{\mathbf{S^2}}{2\beta(\mathbf{K\_s} - \mathbf{K\_0})} \ln\left[1 + \beta \frac{\mathbf{K\_s} - \mathbf{K\_0}}{\mathbf{q\_s(t)} - \mathbf{K\_s}}\right], \mathbf{q\_s} = \frac{\mathbf{dI}}{\mathbf{d}t} \tag{20}$$

where ∆θ = θ<sup>s</sup> − θ<sup>0</sup> is the storage capacity of the soil, h is the water depth on the soil surface, S is the sorptivity, and β is a shape parameter thus that 0 < β < 1, the lower limit corresponds to the Green and Ampt model while the higher limit corresponds to the Talsma and Parlange model [45]. Time in Equation (20) should be interpreted as the contact time that water has at a given point of the border.

Sorptivity S can be calculated with the expression proposed by Parlange [46]:

$$\mathbf{S}^2 = \int\_{\psi\_0}^0 [\theta\_\mathbf{s} + \theta(\psi) - 2\theta\_0] \mathbf{K}(\psi) d\psi \tag{21}$$

and the parameter of β shape can be calculated with the expression proposed by Fuentes [47]:

$$1 - \frac{1}{2}\beta = \frac{\int\_{\Theta\_0}^{\Theta\_s} \left[\frac{\mathbf{K}(\Theta) - \mathbf{K}\_0}{\mathbf{K}\_s - \mathbf{K}\_0}\right] \left(\frac{\Theta\_s - \Theta\_0}{\Theta - \Theta\_0}\right) \mathbf{D}(\Theta) d\Theta}{\int\_{\Theta\_0}^{\Theta\_s} \mathbf{D}(\theta) d\Theta} \tag{22}$$

Variation in time of water depth on the soil is provided by the system of Saint-Venant and Richards equations, but their analytical representation is not known, the reason whereby it is assumed that it is represented by an average value; the mean depth of water can be estimated as a fraction of the normal depth: h = 4/5 hn. With the dimensionless variables:

$$\mathbf{t}^\* = \frac{2(\mathbf{K}\_\mathbf{s} - \mathbf{K}\_0)^2 \mathbf{t}}{\mathbf{S}^2 + 2\mathbf{K}\_\mathbf{s} \overline{\mathbf{h}} \Delta \theta}, \mathbf{I}^\*(\mathbf{t}^\*) = \frac{2(\mathbf{K}\_\mathbf{s} - \mathbf{K}\_0)[\mathbf{I}(\mathbf{t}) - \mathbf{K}\_0 \mathbf{t}]}{\mathbf{S}^2 + 2\mathbf{K}\_\mathbf{s} \overline{\mathbf{h}} \Delta \theta} \tag{23}$$

$$\mathbf{q\_s^\*(t^\*)} = \frac{d\mathbf{l^\*}}{d\mathbf{t^\*}} = \frac{\mathbf{q\_s(t)} - \mathbf{K\_0}}{\mathbf{K\_s - K\_0}}, \ \boldsymbol{\gamma} = \frac{2\mathbf{K\_s}\overline{\mathbf{h}}\Delta\boldsymbol{\theta}}{\mathbf{S^2} + 2\mathbf{K\_s}\overline{\mathbf{h}}\Delta\boldsymbol{\theta}}\tag{24}$$

Equation (20) is written as:

$$\Gamma = \frac{\gamma}{\mathbf{q}\_s^\*-1} + \frac{1-\gamma}{\beta} \ln\left[1 + \frac{\beta}{\mathbf{q}\_s^\*-1}\right] \tag{25}$$

The relationship dI\*/dq∗ <sup>s</sup> = q ∗ s dt\*/dq∗ s , considering constant the water depth on the surface, along with the initial condition t\* = 0, I\* = 0, q ∗ <sup>s</sup> → ∞, leads to find the time as a function of the infiltration flow:

$$\mathbf{t}^\* = \frac{\gamma}{\mathbf{q}\_\mathbf{s}^\* - 1} + \frac{1 - \gamma}{\beta (1 - \beta)} \ln \left[ 1 + \frac{\beta}{\mathbf{q}\_\mathbf{s}^\* - 1} \right] - \frac{1 - \beta \gamma}{1 - \beta} \ln \left[ 1 + \frac{1}{\mathbf{q}\_\mathbf{s}^\* - 1} \right] \tag{26}$$

Thus, the function defining the infiltrated depth in the function of time is of a parametric nature: I\* = I\* (q∗ s ) and t\* = t\* (q∗ s ), with the flow q∗ s as a parameter.

[44]:

(Equation (28)).

Due to the high nonlinearity of the function K(θ) it can be assumed that K<sup>0</sup> << Ks. The dimensionless version of α<sup>u</sup> defined in Equation (18) is as follows:

$$\alpha\_{\rm u} = \frac{\ell\_{\rm n}^\*(\mathbf{q}\_{\rm sn}^\*)}{\pi\_{\rm n}^\*(\mathbf{q}\_{\rm sn}^\*)} \tag{27}$$

where the dimensionless net irrigation time and net irrigation depth are defined by Equation (23) and the corresponding dimensionless infiltration flow by Equation (24). For a given border and initial medium depth, the irrigation depth in a dimensionless form (` ∗ n ) was calculated, and the dimensionless flow (q ∗ sn) was calculated with Equation (25) iteratively. Finally, the dimensionless net time (τ ∗ n ) was calculated with Equation (26), and α<sup>u</sup> was calculated with Equation (27). *Water* **2020**, *12*, x FOR PEER REVIEW 10 of 14 where the dimensionless net irrigation time and net irrigation depth are defined by Equation (23) and the corresponding dimensionless infiltration flow by Equation (24). For a given border and initial medium depth, the irrigation depth in a dimensionless form ( ℓ<sup>n</sup> ∗ ) was calculated, and the dimensionless flow (qsn ∗ ) was calculated with Equation (25) iteratively. Finally, the dimensionless net

It should be noted that the process of calculating the optimal flow, for a given irrigation length, was also iterative since the medium depth depended on the normal depth and this, in turn, of the optimal flow, through Equation (8): h<sup>n</sup> = [ν 2 (qo/kν) 1/d /gJo] 1/3 . time (τ<sup>n</sup> ∗ ) was calculated with Equation (26), and α<sup>u</sup> was calculated with Equation (27). It should be noted that the process of calculating the optimal flow, for a given irrigation length, was also iterative since the medium depth depended on the normal depth and this, in turn, of the 1/3

When the water depth was small (h << S 2 /2Ks∆θ) from Equation (25) it was deduced an explicit function of time with respect to the infiltrated water and corresponds to the Parlange et al. equation [44]: optimal flow, through Equation (8): h<sup>n</sup> = [ν<sup>2</sup> (qo/kν)1/d/gJo] . When the water depth was small (h̅ << S2/2Ks∆θ) from Equation (25) it was deduced an explicit function of time with respect to the infiltrated water and corresponds to the Parlange et al. equation

$$\mathbf{t}^\* = \mathbf{I}^\* - (1 - \beta)^{-1} \ln \left[ \beta^{-1} [1 - (1 - \beta) \exp(-\beta \mathbf{I}^\*)] \right] \tag{28}$$
 
$$\text{...r of constraints and the above numerator for different soil are connected.}$$

In Table 4, the values of sorptivity and the shape parameter for different soils are reported. In Figure 3, the optimal infiltration time obtained with Saint-Venant and Richards equations is compared with that obtained with Equation (28) of Parlange et al. [44]: r<sup>2</sup> = 0.9866. In Table 4, the values of sorptivity and the shape parameter for different soils are reported. In Figure 3, the optimal infiltration time obtained with Saint-Venant and Richards equations is compared with that obtained with Equation (28) of Parlange et al. [44]: r <sup>2</sup> = 0.9866.

**Table 4.** Parameters of the Parlange et al. infiltration equation [44]: sorptivity (S) and the shape parameter (β), calculated with Equations (21) and (29). **Table 4.** Parameters of the Parlange et al. infiltration equation [44]: sorptivity (S) and the shape parameter (β), calculated with Equations (21) and (29).


**Barré de Saint-Venant and Richards Infiltration time (s)**

The shape parameter β varies with respect to the initial moisture content; however, it does not

**Figure 3.** Comparison between the optimal irrigation time numerically obtained with the Barré de Saint-Venant/Richards equations and those calculated ones with the Parlange et al. equation **Figure 3.** Comparison between the optimal irrigation time numerically obtained with the Barré de Saint-Venant/Richards equations and those calculated ones with the Parlange et al. equation (Equation (28)).

The shape parameter β varies with respect to the initial moisture content; however, it does not vary significantly when this moisture content approaches the residual moisture content. In this case the introduction of the hydrodynamic characteristics defined by Equations (4) and (5) in Equation (22) provides: the introduction of the hydrodynamic characteristics defined by Equations (4) and (5) in Equation (22) provides: B 2η 1 m 1 n ,1 n β 2 1 (29)

*Water* **2020**, *12*, x FOR PEER REVIEW 11 of 14

$$\beta = 2 \left\{ 1 - \frac{\mathbb{B}[(2\eta - 1)\mathbf{m} - 1/\mathbf{n}, 1/\mathbf{n}]}{\mathbb{B}(\eta \mathbf{m} - 1/\mathbf{n}, 1/\mathbf{n})} \right\} \tag{29}$$

where B(p,q) is the complete beta function. where B(p,q) is the complete beta function. Introduction of Equation (28) into Equation (27) gives the approximate formula for calculating

Introduction of Equation (28) into Equation (27) gives the approximate formula for calculating the optimal unitary flow in function of the border length, the net depth, and the characteristic parameters of infiltration representing the capillary forces (sorptivity) and gravitational forces (hydraulic conductivity at saturation) and a shape parameter of the hydrodynamic characteristics, namely: the optimal unitary flow in function of the border length, the net depth, and the characteristic parameters of infiltration representing the capillary forces (sorptivity) and gravitational forces (hydraulic conductivity at saturation) and a shape parameter of the hydrodynamic characteristics, namely:

$$\alpha\_{\rm ul} = \frac{\ell\_{\rm n}^{\star}}{\ell\_{\rm n}^{\star} - (1 - \beta)^{-1} \ln \left[ \beta^{-1} [1 - (1 - \beta) \exp(-\beta \ell\_{\rm n}^{\star})] \right]}, \; \ell\_{\rm n}^{\star} = \frac{2 \mathbf{K}\_{\rm s} \ell\_{\rm n}}{\mathbf{S}^{2}} \tag{30}$$

where ` ∗ n is the net irrigation depth in dimensionless writing. where ℓ<sup>n</sup> ∗ is the net irrigation depth in dimensionless writing.

As can be seen, Equation (30) contains a shape parameter in the function of the soil type; however, it does not have a great variation in the range of soils reported in Table 4, the mean value β ≈ 3/4 can be taken. With this value in Figure 4, the graph of Equation (30) and its comparison with the numerical results are presented. A good agreement is clearly demonstrated. As can be seen, Equation (30) contains a shape parameter in the function of the soil type; however, it does not have a great variation in the range of soils reported in Table 4, the mean value β 3/4 can be taken. With this value in Figure 4, the graph of Equation (30) and its comparison with the numerical results are presented. A good agreement is clearly demonstrated.

**Figure 4.** The multiple factors of the minimum discharge as a function of the irrigation depth numerically obtained with the Saint-Venant and Richards equations, and those calculated with the Parlange et al. model [44], Equation (30). **Figure 4.** The multiple factors of the minimum discharge as a function of the irrigation depth numerically obtained with the Saint-Venant and Richards equations, and those calculated with the Parlange et al. model [44], Equation (30).

#### *3.3. Application of the Analytical Formula in Furrow Irrigation Systems 3.3. Application of the Analytical Formula in Furrow Irrigation Systems*

This analytical formula has been applied in field experiments realized by Chávez and Fuentes [1] with good results. Irrigation tests (250) were performed in 1010 ha with the next crops: *Zea mayz, Sorghum vulgare, Medicago sativa, Phaseolus vulgaris, Pachyrhizus erosus, Hordeum vulgare, Triticum aestivum,* and *Allium cepa*. This analytical formula has been applied in field experiments realized by Chávez and Fuentes [1] with good results. Irrigation tests (250) were performed in 1010 ha with the next crops: *Zea mayz*, *Sorghum vulgare*, *Medicago sativa*, *Phaseolus vulgaris*, *Pachyrhizus erosus*, *Hordeum vulgare*, *Triticum aestivum*, and *Allium cepa*.

The characteristics and properties that were measured in the plots were: Length, slope, texture, bulk density, initial, and saturation moisture content. The results that were obtained from the evaluation of the irrigation tests were: Discharge at the entrance of the plot, number of furrows by The characteristics and properties that were measured in the plots were: Length, slope, texture, bulk density, initial, and saturation moisture content. The results that were obtained from the evaluation of the irrigation tests were: Discharge at the entrance of the plot, number of furrows by irrigation set,

saturated hydraulic conductivity, the efficiency of the application, and the slope of the direction of irrigation. The irrigation tests were performed in plots.

With the advanced and recession data and the characteristics of the soils from the locations where the irrigation tests were carried out, the calibration of the test was performed using the kinematic wave model. With the parameters found (K<sup>s</sup> and h<sup>f</sup> ) for the evaluation of the irrigation tests and the net irrigation depth that is intended to be applied on the plot (water depth depending on each of the crops established in the plots), Equations (18) and (30) were used to make the design. The obtained result is the optimal flow that must go into each furrow; for this value, the discharge for the entrance of the plot between that obtained with Equation (18) was divided, and then approached the nearest whole value. As a result, the equation gives us the number of furrows that the user has to open by set and the time that must pass before cutting off the water.

In this study, the evaluation of irrigation tests, the data of the plot, and the net irrigation depth to be applied demonstrated that the optimal flow expense that can be put in each furrow during an irrigation event can be calculated under the hypothesis that with this expense, the historical water depths applied in the evaluated plots can be reduced. The average water depths decreased by 19 cm, irrigation time decreased 12 h ha−<sup>1</sup> on average, and the average volume saved was 2150 m<sup>3</sup> ha−<sup>1</sup> , which represented a total of 49% of the total volume used. In addition, the average efficiency rose from 51% to 86%.

#### **4. Conclusions**

A linear relationship has been validated between the length of the border and the optimal irrigation flow, defined as the inflow rate that has to be applied to obtain a maximum value of Christiansen's uniformity coefficient with high values of application efficiency and irrigation requirement efficiency. The linear form of the proportion between both variables was corroborated by [39] using a hydrological model for the flow of surface water and the Green and Ampt infiltration equation.

The proportion between optimal flow and border length has been established by numerical simulation with the hydrodynamic model, formed by the coupled equations of Barré de Saint-Venant and Richards. In the numerical simulation, 10 types of soils of contrasting texture have been used and three water depths applied, which has allowed us to form an irrigation design table, that is, 10 linear relationships for each irrigation depth.

To establish an analytical representation of the optimal flow in function of the border length, the Parlange infiltration theory proposed for integrating the differential equation of the vertical one-dimensional infiltration has been used. In the formula obtained to relate the optimal irrigation flow and length, the irrigation depth intervenes and as soil parameters sorptivity that comes from the capillary forces and the saturated hydraulic conductivity that comes from gravitational forces. The good accordance between the results of numerical simulation and analytical representation allows us to obtain the formula for the design of gravity irrigation.

**Author Contributions:** Conceptualization, C.F.; methodology, C.F. and C.C.; software, C.F.; validation, C.C. and C.F.; formal analysis, C.C.; investigation, C.F.; resources, C.F. and C.C.; data curation, C.F. and C.C.; writing—original draft preparation, C.F.; writing—review and editing, C.C.; visualization, C.F. and C.C.; supervision, C.F.; project administration, C.F. and C.C.; funding acquisition, C.F. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).
