Modeling of Border Irrigation in Soils with the Presence of a Shallow Water Table. I: The Advance Phase

: The overelevation of the water table in surface irrigation plots is one of the main factors affecting salinization in agricultural soils. Therefore, it is necessary to develop simulation models that consider the effect of a shallow water table in the process of advance-inﬁltration of the water in an irrigation event. This paper, the ﬁrst in a series of three, develops a simple mathematical model for the advance phase of border irrigation in soils with the presence of a shallow water table. In this study, the hydrodynamic model of the Barr é de Saint-Venant equations is used for the water surface ﬂow, and the equations are solved using a Lagrangian ﬁnite-differences scheme, while in the subsurface ﬂow, an analytical solution for inﬁltration in soils with a shallow water table is found using the bisection method to search for roots. In addition, a hydraulic resistance law is used that eliminates the numerical instabilities presented by the Manning–Strickler law. The model results for difference irrigation tests show adjustments with an R 2 > 0.98 for the cases presented. It is also revealed that, when increasing the time step, the precision is maintained, and it is possible to reduce the computation time by up to 99.45%. Finally, the model proposed here is recommended for studying the advance process during surface irrigation in soils with shallow water tables.


Introduction
Surface irrigation is one of the main factors in the fluctuation of shallow water tables in agricultural plots. Around the world, surface irrigation is the most frequently used method for water application in agricultural plots, where the water is distributed over the field and through the soil [1].
The presence of a shallow water table changes the water content of the soil profile ( Figure 1), due principally to the capillary process, which is present in porous media [2]. It has been found that in shallow water conditions, irrigation can be reduced by up to 80% without affecting yield and without increasing soil salinity [3].
Poor design of surface irrigation, for example losses from tailwater due to the selection of an inappropriate irrigation flow or due to poor operation in the water distribution network as a result of poor leveling in the ground by percolation, can cause overelevation of the water table and the progressive salinization of the soil [4].
It is estimated that 20% of the arable land and 33% of irrigated agricultural plots in the world are affected by high salinity [5]. Water-table levels play an essential role in the salt distribution in the soil profile and could be controlled by subsurface drainage [6]. However, drainage leads to an increase in irrigation costs. Therefore, it is necessary to model surface irrigation both for the surface water movement and infiltration into the soil, using an equation that considers soil with a shallow water table for the latter. The border irrigation water flow can be characterized using the hydrodynamic model of the Barré de Saint-Venant equations, which accurately describes surface irrigation [7]. The relationship between the width and the water depth in a border makes it possible to consider the equations corresponding to runoff for an infinite surface width [8]. The continuity equation in the hydrodynamic model is written as follows: x t (1) and the momentum equation is as follows [9]: where q(x,t) = U(x,t)h(x,t) is the discharge per unit width of the border or the unitary discharge, x is the spatial coordinate in the main direction of the water movement in the border, t is the time, U is the mean velocity, h is the water depth, VI = ∂I(x,t)/∂t is the infiltration flow, that is, the water volume infiltrated per unit width per unit length of the border, I is the infiltrated depth, g is gravitational acceleration, β = UIX/U is a dimensionless parameter where UIX is the projection in the direction of the output velocity of the water mass due to the infiltration, Jo is the topographic slope, and J is the friction slope that can be determined by the fractal law of hydraulic resistance [10]: where ν is the kinematic viscosity coefficient, k is a dimensionless factor that includes the effects of soil roughness, and the exponent d has a fractal interpretation. From this law, the Chezy formula is deduced with d = 1/2 and the Poiseuille law with d = 1. The initial and boundary conditions for a closed border, avoiding loss by tailwater outside the irrigation domain, are as follows: where xf (t) is the position of the wave front at time t and qo is the unitary discharge at the entrance of the border. The border irrigation water flow can be characterized using the hydrodynamic model of the Barré de Saint-Venant equations, which accurately describes surface irrigation [7]. The relationship between the width and the water depth in a border makes it possible to consider the equations corresponding to runoff for an infinite surface width [8]. The continuity equation in the hydrodynamic model is written as follows: and the momentum equation is as follows [9]: where q(x,t) = U(x,t)h(x,t) is the discharge per unit width of the border or the unitary discharge, x is the spatial coordinate in the main direction of the water movement in the border, t is the time, U is the mean velocity, h is the water depth, V I = ∂I(x,t)/∂t is the infiltration flow, that is, the water volume infiltrated per unit width per unit length of the border, I is the infiltrated depth, g is gravitational acceleration, β = U IX /U is a dimensionless parameter where U IX is the projection in the direction of the output velocity of the water mass due to the infiltration, J o is the topographic slope, and J is the friction slope that can be determined by the fractal law of hydraulic resistance [10]: where ν is the kinematic viscosity coefficient, k is a dimensionless factor that includes the effects of soil roughness, and the exponent d has a fractal interpretation. From this law, the Chezy formula is deduced with d = 1/2 and the Poiseuille law with d = 1. The initial and boundary conditions for a closed border, avoiding loss by tailwater outside the irrigation domain, are as follows: where x f (t) is the position of the wave front at time t and q o is the unitary discharge at the entrance of the border.
Numerous practical situations require a numerical solution of the Richards equation, with dimensions depending on the complexity of the studied problem. In surface irrigation, a one-dimensional equation is sufficient to represent the infiltration phenomenon [11]. However, the equation lacks general analytical solutions, and therefore complex numerical methods are often required in order to solve it [12][13][14]. There are other physics-based models resulting from the simplification of the initial conditions, in particular the Green and Ampt equation [15]. This equation has been used in surface irrigation [16][17][18]. However, with this equation, only a homogeneous moisture profile can be represented in the soil profile.
Fuentes et al. [19] developed an analytical solution of the Richards equation using the Green and Ampt hypotheses to describe the infiltration of water into the soil with a shallow water table (P f ). The solution considers a hydrostatic initial moisture distribution (Figure 2), where the initial moisture is calculated with the expression θ i (z) = θ o + (θ s − θ o )(z/P f ). Thus, the moisture content at the soil surface is θ i (0) = θ o and at the water-table surface is θ i (P f ) = θ s . The suction in the wetting front is a linear function of the moisture content at the front, h f (θ i ,θ s ) = h f (θ s − θ i )/(θ s − θ o ), that is: such that h f [θ i (0),θ s ] = h f and h f [θ i (P f ),θ s ] = 0. The infiltrated depth is defined by: that is: where I M is the maximum infiltrated depth. Numerous practical situations require a numerical solution of the Richards equation, with dimensions depending on the complexity of the studied problem. In surface irrigation, a one-dimensional equation is sufficient to represent the infiltration phenomenon [11]. However, the equation lacks general analytical solutions, and therefore complex numerical methods are often required in order to solve it [12][13][14]. There are other physicsbased models resulting from the simplification of the initial conditions, in particular the Green and Ampt equation [15]. This equation has been used in surface irrigation [16][17][18]. However, with this equation, only a homogeneous moisture profile can be represented in the soil profile. Fuentes et al. [19] developed an analytical solution of the Richards equation using the Green and Ampt hypotheses to describe the infiltration of water into the soil with a shallow water table (Pf). The solution considers a hydrostatic initial moisture distribution (Figure 2), where the initial moisture is calculated with the expression θi(z) = θo + (θs − θo)(z/Pf). Thus, the moisture content at the soil surface is θi(0) = θo and at the water-table surface is θi(Pf) = θs. The suction in the wetting front is a linear function of the moisture content at the front, hf (θi,θs) = hf (θs − θi)/(θs − θo), that is: such that hf[θi(0),θs] = hf and hf[θi(Pf),θs] = 0. The infiltrated depth is defined by: that is: where IM is the maximum infiltrated depth. Introducing Equation (6) into the Green and Ampt equation yields the differential equation for the analytical solution that considers a shallow water table [19]: Introducing Equation (6) into the Green and Ampt equation yields the differential equation for the analytical solution that considers a shallow water table [19]: When the water depth is independent of time, h = h, and integration of Equation (9) with the condition I = 0 at t = 0 leads to the following infiltration equation, where h f = P f : and when h f = P f : where h ∼ = [4d/(5d + 1)]h o is the mean water depth and h o = (ν 2 /gJ o ) 1/3 (q o /kν) 1/3d is the normal depth [20]. Equation (10) reduces to the Green and Ampt infiltration equation when P f → ∞, considering that I M = 1/2∆θP f . In this limit, 1 -(1 -I/I M ) 1/2 ∼ = I/2I M = I/∆θP f holds. The third term is of the order of 1/P f and tends to zero. In the second term, the argument of the logarithm tends to 1 + I/∆θ(h + h f ) and its coefficient to ∆θ(h + h f ). Finally, in the first term, the coefficient of I → 1. Using the definition λ = ∆θ (h + h f ), the Green and Ampt equation is deduced.
In recent years, many different software packages have been developed to model surface irrigation (e.g., [21,22]); however, they have some limitations due to the infiltration equations used. For example, they only represent constant initial moisture conditions along a homogeneous soil column, and they are empirical equations for a specific irrigation event. Furthermore, they are not representative for the soils with shallow water tables (where the moisture profile is not constant) existing in some irrigated agricultural areas.
The objectives of this study, the first in a series of three, are: (a) to model the advance phase of surface irrigation in a soil with a shallow water table by coupling an analytical solution to the Barré de Saint-Venant equations and (b) to validate the model obtained with data from an irrigation test reported in the literature. In the second paper, the three phases of irrigation will be included: advance, storage, and recession. Finally, in the third article, this solution will be compared with the Barré de Saint-Venant equations coupled internally with the Richards equation.

Numerical Solution
The numerical solution of the Barré de Saint-Venant equations uses a Lagrangian scheme [9,23]. The discrete form of the continuity equation is as follows: Due to the numerical instability presented in the calculation cells at the beginning, a discrete form of the momentum equation was developed, considering the following Agriculture 2022, 12, 426 5 of 12 assumptions: (a) the derivatives in time were calculated as a rectangular cell (Eulerian) and (b) average flow and water depth coefficients were considered for the calculation of an average friction slope. The discrete form of the momentum equation is as follows: where Rc and Rm represent the residuals of the continuity equation and momentum, respectively, the coefficient q , taking into account the extreme values of each calculation cell, and consequently the coefficient J = ν 2 (q/kν) 1/d /gh 3 .
The weight factors for time and space are denoted ω and ϕ, respectively. In the discrete forms, a weight factor in space of ϕ = 1/2 was considered for interior cells [23][24][25][26]. For the last cell and the first time step, ϕ = π/4 was used, deduced from the analysis of the short-time coupling of the Saint-Venant and Richards equations [27]. The weight factor for time was taken as ω = 0.60 [24][25][26]28]. Figure 3 shows the Lagrangian cells during the advance phase. The subscripts L and R denote the values of each variable to the left and right of the next time t i+1 , and J and M represent the values at the current time t i .
where Rc and Rm represent the residuals of the continuity equation and momentum, re- , and the coef- In the discrete forms, a weight factor in space of φ = 1/2 was considered for interior cells [23][24][25][26]. For the last cell and the first time step, φ = π/4 was used, deduced from the analysis of the short-time coupling of the Saint-Venant and Richards equations [27]. The weight factor for time was taken as ω = 0.60 [24][25][26]28]. Figure 3 shows the Lagrangian cells during the advance phase. The subscripts L and R denote the values of each variable to the left and right of the next time ti+1, and J and M represent the values at the current time ti. The solution developed is implicit, where Equations (12) and (13) are solved simultaneously for all cells at each time step. By linearizing the discrete equations, a system of 2N + 2 nonlinear algebraic equation with 2N + 4 unknown variables is produced, where N is the number of computational cells. The solution of the system is obtained by applying the double sweep algorithm, which assumes a linear relation between the variation of the water depth δh and the discharge δq in a given time step [29]: and for the last cell the equation is expressed in terms of the advance distance, considering that δhN and δqN are equal to zero: The solution developed is implicit, where Equations (12) and (13) are solved simultaneously for all cells at each time step. By linearizing the discrete equations, a system of 2N + 2 nonlinear algebraic equation with 2N + 4 unknown variables is produced, where N is the number of computational cells. The solution of the system is obtained by applying the double sweep algorithm, which assumes a linear relation between the variation of the water depth δh and the discharge δq in a given time step [29]: Agriculture 2022, 12, 426 6 of 12 and for the last cell the equation is expressed in terms of the advance distance, considering that δh N and δq N are equal to zero: where the values of E and F are calculated by the partial derivatives of the residual equations of continuity and momentum for each cell with respect to h L , q L , h R , and q R .
With the coefficient values E and F, the improved solution for the new values of flow, water depth, and position of the wave front is determined by: The convergence criterion selected for advancing in time is that the residual values of the continuity equation (Rc) and the momentum equation (Rm) in the current iteration must be less than 1 × 10 −5 .
The analytical solution for infiltration in a soil with a shallow water table is obtained using the bisection method, which is widely used to find the roots of a function.

Soil Characterization
Simulations were performed with the information reported by Pacheco [30] for a clay soil from La Chontalpa, Tabasco, Mexico, using the coupling of the Barré de Saint-Venant equations and the analytical solution to model infiltration in a soil with a shallow water table. The measured data from the irrigation tests were: border width B = 10.5 m, border length L = 100 m, border slope J o = 0.00085 m/m, moisture content at saturation θ s = 0.5245 cm 3 /cm 3 , and the dimensionless parameter β = 0. The values of the hydraulic conductivity at saturation (K s ) and the suction in the wetting front (h f ) were optimized using the Levenberg-Marquardt algorithm [31]. In the hydraulic resistance law, Equation (3), d = 1 was used. The values of q o , P f , I M , initial moisture (θ o ), h, K s , and h f are reported in Table 1 for each irrigation test. This table shows the correlation coefficient R 2 between the measured data and those calculated with the optimized parameters of K s and h f . Good adjustments were observed for the three irrigation tests.

Applications
Simulations corresponding to each irrigation test were performed using the values reported in Table 1 and the optimized values of K s and h f . Figure 4 shows the good fit of the measured data for the advance phase and those obtained by the proposed model when optimizing the values of K s and h f , which is corroborated by the high R 2 values obtained.

Comparison of the Theorical Advance Curve with the Measured Data
The information from the second irrigation test (Pf = 50 cm) was used for a more exhaustive analysis. Figure 5 shows the advance curve obtained by the model, which shows a good fit between the measured data obtained in the field and those obtained with the model.  Figure 6 shows that when the infiltrated depth reaches the IM calculated with Equation (8), the soil is no longer capable of storing more water. Therefore, the soil was considered completely saturated 36 min after irrigation had started.

Comparison of the Theorical Advance Curve with the Measured Data
The information from the second irrigation test (P f = 50 cm) was used for a more exhaustive analysis. Figure 5 shows the advance curve obtained by the model, which shows a good fit between the measured data obtained in the field and those obtained with the model.

Comparison of the Theorical Advance Curve with the Measured Data
The information from the second irrigation test (Pf = 50 cm) was used for a more exhaustive analysis. Figure 5 shows the advance curve obtained by the model, which shows a good fit between the measured data obtained in the field and those obtained with the model.  Figure 6 shows that when the infiltrated depth reaches the IM calculated with Equation (8), the soil is no longer capable of storing more water. Therefore, the soil was considered completely saturated 36 min after irrigation had started.

Distance (m)
Measured data Model development Figure 5. Evolution of the advance front in the second irrigation test. Figure 6 shows that when the infiltrated depth reaches the I M calculated with Equation (8), the soil is no longer capable of storing more water. Therefore, the soil was considered completely saturated 36 min after irrigation had started.
The surface and subsurface flow profiles are shown in Figure 7, where it can be observed that the infiltration over all points of the border is limited by the maximum infiltration shown in Figure 6. It is also evident from the third curve (39 min), that a determined region of the soil was already at saturation, because the soil could no longer store more water. Due to the complexity of obtaining the roots in the analytical solution, the time used for this process was 55 min; however, the time discretization can be increased to obtain reliable results in a shorter time. The surface and subsurface flow profiles are shown in Figure 7, where it can be observed that the infiltration over all points of the border is limited by the maximum infiltration shown in Figure 6. It is also evident from the third curve (39 min), that a determined region of the soil was already at saturation, because the soil could no longer store more water. Due to the complexity of obtaining the roots in the analytical solution, the time used for this process was 55 min; however, the time discretization can be increased to obtain reliable results in a shorter time.

Time Step Analysis
An analysis of the computation time was performed at different values of δt in the numerical coupling of the Barré de Saint-Venant equations and the analytical solution of infiltration with a shallow water table. Table 2 shows that there was no difference for the different values of δt studied here. The computer equipment used to perform the simulations has an Intel ® Core TM i7-4710 CPU @ 2.50 GHz and 32 Gb of RAM.  The surface and subsurface flow profiles are shown in Figure 7, where it can be observed that the infiltration over all points of the border is limited by the maximum infiltration shown in Figure 6. It is also evident from the third curve (39 min), that a determined region of the soil was already at saturation, because the soil could no longer store more water. Due to the complexity of obtaining the roots in the analytical solution, the time used for this process was 55 min; however, the time discretization can be increased to obtain reliable results in a shorter time.

Time Step Analysis
An analysis of the computation time was performed at different values of δt in the numerical coupling of the Barré de Saint-Venant equations and the analytical solution of infiltration with a shallow water table. Table 2 shows that there was no difference for the different values of δt studied here. The computer equipment used to perform the simulations has an Intel ® Core TM i7-4710 CPU @ 2.50 GHz and 32 Gb of RAM.

Time Step Analysis
An analysis of the computation time was performed at different values of δt in the numerical coupling of the Barré de Saint-Venant equations and the analytical solution of infiltration with a shallow water table. Table 2 shows that there was no difference for the different values of δt studied here. The computer equipment used to perform the simulations has an Intel ® Core TM i7-4710 CPU @ 2.50 GHz and 32 Gb of RAM.  According to Figure 8, for design purposes it is recommended to use δt = 10 s to increase the processing speed, and this give a fast response for the advance/infiltration process in each plot where irrigation was required, with presence of the water table. According to Figure 8, for design purposes it is recommended to use δt = 10 s to increase the processing speed, and this give a fast response for the advance/infiltration process in each plot where irrigation was required, with presence of the water table. For the purpose of observing whether, with a temporal increment, the accuracy of the solution was lost, simulations were performed for different values of δt. However, the coupling method used in this work yielded excellent results and when performing the temporal increments the errors were minimal (Figure 9). Nevertheless, it is recommended to use a small δt for research purposes, since it is necessary to observe the phenomenon in detail in order to provide solutions for particular cases that may occur in irrigation areas with these characteristics. For the purpose of observing whether, with a temporal increment, the accuracy of the solution was lost, simulations were performed for different values of δt. However, the coupling method used in this work yielded excellent results and when performing the temporal increments the errors were minimal (Figure 9). Nevertheless, it is recommended to use a small δt for research purposes, since it is necessary to observe the phenomenon in detail in order to provide solutions for particular cases that may occur in irrigation areas with these characteristics. According to Figure 8, for design purposes it is recommended to use δt = 10 s to increase the processing speed, and this give a fast response for the advance/infiltration process in each plot where irrigation was required, with presence of the water table. For the purpose of observing whether, with a temporal increment, the accuracy of the solution was lost, simulations were performed for different values of δt. However, the coupling method used in this work yielded excellent results and when performing the temporal increments the errors were minimal (Figure 9). Nevertheless, it is recommended to use a small δt for research purposes, since it is necessary to observe the phenomenon in detail in order to provide solutions for particular cases that may occur in irrigation areas with these characteristics.

Discussion
The numerical coupling of the Barré de Saint-Venant equations with the analytical solution for a soil with a shallow water table allows the effect of this on the infiltration

Discussion
The numerical coupling of the Barré de Saint-Venant equations with the analytical solution for a soil with a shallow water table allows the effect of this on the infiltration process to be analyzed and related to the maximum infiltrated depth of the soil and the most suitable crop for soils with this characteristic, as well as showing how the advance wave process is affected at the surface.
By optimizing the parameters K s and h f , it was possible to represent surface irrigation in soils with shallow water tables, obtaining excellent results for the coefficient of determination for each irrigation test, which confirms that the proposed model adequately reproduces the data measured in the field. In addition, with the values obtained it is possible to calculate, for the next irrigation, the optimal irrigation flow that should be applied in each border, by means of an analytical representation that takes into account the border length, the net irrigation depth, and the characteristic parameters of infiltration obtained through the inverse process shown here [32].
Solving the numerical coupling of the model developed by the explicit method has a high computational cost. However, by having δt > 1 s, the computational time decreases by up to 99.45% compared to δt = 1 s, and the results provided by the model still fit correctly. This method makes it possible to increase the processing speed and to ideally reproduce the advance/infiltration process in the border irrigation in soils with shallow water tables. It is important to note that for research purposes it is advisable to use a small δt to observe in detail what happens in the advance phase of irrigation.
By using an implicit scheme for the discretization of the equations, a higher degree of accuracy is provided, because the solution is iterative and considers what is occurring in the current state of the flow profiles, as well as in the previous state. In contrast, the explicit scheme only solves the system once, without considering the current state of the flow profiles, which can generate accuracy problems when using different values of δt.
The developed model has the advantage of not requiring as many physical parameters associated with the plots and irrigation, without losing the representativeness of the soils, while the use of the Richards equation requires a hydrodynamic characterization of the soil including the characteristic moisture curve and the hydraulic conductivity curve. This represents more exhaustive characterization work and more time to obtain solutions in the process of border irrigation, because the equation requires a complex numerical solution [11].
Our method is therefore a robust and efficient model for the advance phase in soils with shallow water table, compared to others reported in the literature that use infiltration equations without a physical basis, such as the Kostiakov (e.g., [21,23]) and Kostiakov-Lewis equations (e.g., [22,33]).

Conclusions
A numerical solution of the Barré de Saint-Venant equations coupled to the analytical solution of the infiltration for a soil with a shallow water table was implemented to describe the advance phase of border irrigation. The proposed solution of the model used a Lagrangian finite-difference scheme for the Barré de Saint-Venant equations, while the infiltration equation was solved by the bisection method. The model evaluation in the simulation of the advance phase revealed that the model fits a set of experimental data.
Finally, the solution shown here can be used with excellent results to design and model surface irrigation, since the infiltration equation used for a soil with a shallow water table requires only a few characteristic soil parameters (K s and h f ) that are easy to obtain.
Author Contributions: S.F. and C.C. contributed equally to this work. S.F. and C.C. All authors have read and agreed to the published version of the manuscript.