Modeling Two-Dimensional Infiltration with Constant and Time-Variable Water Depth

Abstract

Water infiltration is simulated by obtaining the time infiltrated depth evolution and humidity profiles with the numerical solution of the two-dimensional Richards’ equation. The contact time hypothesis is accepted in this study and used to apply a unique form on time of the water depth evolution in the solution domain (furrow), as boundary condition. The specific form of such evolution in time was obtained from results reported in the literature based on the internal numerical full coupling of the Saint-Venant and Richards’ equations in border irrigation. Moreover, the equivalent hydraulic area between the border and the furrow was achieved by scaling the values of water depth. The analysis was made for three contrasting soil textures, and the comparison was done by computing the root mean square error (RMSE) indicator. The comparison was performed from the selection of five finite element meshes with different densities to discretize the solution domain of the two-dimensional Richards’ equation, combined with several time steps. Finally, a comparison was made between infiltrated depth evolution calculated with a constant water depth in the furrow to the one proposed in this work, finding important differences between both approaches. To expand the scope of this study and for a fuller exploration of the subject, the results were compared with results obtained by applying the HYDRUS-2D software. The results confirm that it is important to consider an internal full coupling of the Saint-Venant and Richards´ equations to improve furrow irrigation simulations.

1. Introduction

Richards’ equation [1] is used in most hydrological models to describe groundwater flow in porous media. Richards’ equation is derived from the combination of a continuity equation and velocity ranges obtained with Darcy’s law; the three-dimensional form without water extraction by plants is:

$$\mathrm{C}\left(\mathsf{\psi }\right)\frac{\partial \mathsf{\psi }}{\partial \mathrm{t}}=\nabla ·\left[\mathrm{K}\left(\mathsf{\psi }\right)\nabla \mathsf{\psi }\right]-\frac{\mathrm{dK}}{\mathrm{d}\mathsf{\psi }}\frac{\partial \mathsf{\psi }}{\partial \mathrm{z}}$$

(1)

where $\mathsf{\psi }$ is the water pressure potential into soil, expressed as the height of an equivalent water column (L) (positive in saturated zones and negative in unsaturated zones of soil); $\mathrm{C}\left(\mathsf{\psi }\right)=\mathrm{d}\mathsf{\theta }/\mathrm{d}\mathsf{\psi }$ is the humidity specific capacity of soil; $\mathsf{\theta }=\mathsf{\theta }\left(\mathsf{\psi }\right)$ is the water volume per unit volume of soil or water content; (L³·L⁻³) is a $\mathsf{\psi }$ function, known as the humidity characteristic curve or water retention curve; $\mathrm{K}=\mathrm{K}\left(\mathsf{\psi }\right)$ is hydraulic conductivity (L·T⁻¹), in a saturated soil as a function of pressure potential; gravitational potential is assigned to spatial coordinate $\mathrm{z}$ positively oriented down (L); $\nabla =\left(\partial /\partial \mathrm{x},\partial /\partial \mathrm{y},\partial /\partial \mathrm{z}\right)$ is the gradient operator; $\mathrm{x}$ and $\mathrm{y}$ are spatial coordinates (L), and $\mathrm{t}$ is time (T).

In the literature, studies on infiltration modeling in surface irrigation can be found; for example, Liu et al. [2] proposed a coupled model in which surface water flow and solute transport are described using the zero-inertia equation and the average cross-sectional convection–dispersion equation, respectively, while the two-dimensional Richards’ equation and the convection–dispersion equation are used to simulate water flow and solute transport in soils, respectively. Dong et al. [3] developed a new numerical methodology based on the physical-based model of surface irrigation and the Monte Carlo simulation method to improve the modeling accuracy of surface irrigation performance at a field scale. Border irrigation was simulated by coupling the Saint-Venant equations and a one-dimensional Richards’ equation. Duncan et al. [4] presented a mathematical model that describes the movement of water and solutes in a ridge and furrow geometry. They focus on the effects of two physical processes: root water uptake and pond formation in the furrows. Special attention was paid to the physical description of ponding as an effect of transient rain events. Fan et al. [5] used vertical line source irrigation as a water-saving irrigation method in order to enhance water and nutrient delivery to the root zone. Therefore, reduction on soil evaporation and improvement on water and nutrient efficiency was achieved. To identify influencing factors, they performed computer simulations using the HYDRUS-2D software. Brunetti et al. [6] studied a hybrid Finite Volume–Finite Element (FV–FE) model that describes the coupled surface–subsurface flow processes occurring during furrow irrigation and fertigation. Furman [7] discussed some research gaps, led by the need to include vertical momentum transfer and to expand the use of fully coupled models toward surface irrigation applications. Banti et al. [8] developed a model for the simulation of furrow irrigation advance based on the Saint-Venant equations for one-dimensional surface flow and the two-dimensional Richards’ equation for porous media flow. Bautista et al. [9] proposed a methodology for estimating furrow infiltration under time-variable ponding depth. The methodology approximated the solution to the two-dimensional Richards’ equation, and was a modification of a procedure that was originally proposed for computing infiltration under constant ponding depth. Bautista et al. [10] expanded the analysis of a proposed furrow infiltration formulation based on an approximate solution to the two-dimensional Richards’ equation. The approach calculates two-dimensional infiltration flux as the sum of one-dimensional infiltration and a second term labeled the edge effect. Dong et al. [11] presented a hybrid coupling strategy by distinguishing the mass conservation and momentum conservation of surface flow and subsurface flow systems to reduce the iteration times of the momentum conservation equation and increase the simulation accuracy of coupled models. They proposed a coupling strategy, where the mass conservation component of surface flow model and subsurface flow model are iteratively coupled at the same time step to obtain the convergence value of surface flow depth, and then the momentum conservation component of surface flow model was externally coupled based on the convergence value of both the surface flow depth and infiltration rate to update the surface flow velocity. In the hybrid coupling procedure, infiltration was used as the common internal boundary condition. Soroush et al. [12] developed a furrow irrigation model based on the slow-change/slow-flow routing equation, which is an approximate reduced form of the Saint-Venant equations to a single equation with a single variable, the upstream volume of water. Downstream-propagating disturbances show that the rate of change of upstream inflow is small, with no limit on the Froude number, so it can be used for all slopes and with all common end conditions. To calculate resistance to flow, a composite model in terms of almost any boundary roughness was proposed. Infiltration is assumed to follow the Kostiakov formula. Wöhling et al. [13] presented a process based on a seasonal furrow irrigation model which describes the interacting one-dimensional surface–two-dimensional subsurface flow and crop growth during a whole growing period. The irrigation advance model was based on an analytical solution of the zero-inertia surface flow equations and was iteratively coupled with the two-dimensional subsurface flow model HYDRUS-2D. Wöhling et al. [14] showed a proposed furrow advance phase model (FAPS) which further develops the concepts of a previous study. An analytical zero-inertia surface flow model was iteratively coupled with the two-dimensional (2D) subsurface water transport model HYDRUS-2D. Wöhling et al. [15] proposed a model that overcomes the restrictions of traditional furrow irrigation models by rigorously describing the subsurface flow at computational nodes using the physically based two-dimensional (2D) model HYDRUS-2D. Surface flow is portrayed by an analytical zero-inertia model. Saucedo et al. [16] coupled one-dimensional Saint-Venant equations and two-dimensional Richards’ equation in border irrigation modeling.

Most of the above cited studies apply simplified models of either surface or subsurface water flow. The aim for this study is to note the disadvantage of using water depth average in two-dimensional infiltration simulations, and moreover to compare the infiltrated depth evolution in two-dimensional infiltration, in two cases: (a) considering a constant water depth in the furrow, and (b) considering a time–water depth evolution adapted as reported in literature in relation to a numerical full coupling of Saint-Venant equations and Richards’ equation in border irrigation. Otherwise, the importance of possible differences are analyzed comparing them to those calculated only by the spatial-temporal discretization of the two-dimensional Richards’ equation. This study focuses on phenomena taking place in the furrow cross-section, not along the furrow.

2. Materials and Methods

Irrigation is a three-dimensional phenomenon, therefore, water flux into the soil should be described as shown in Equation (1). This study accepts the contact time hypothesis, which implies that water has a predominantly vertical movement in the soil, i.e., the infiltration law is unique along the furrow. For the above it is possible to use the two-dimensional Richards’ equation form:

$$\mathrm{C}\left(\mathsf{\psi }\right)\frac{\partial \mathsf{\psi }}{\partial \mathrm{t}}=\frac{\partial }{\partial \mathrm{x}}\left[\mathrm{K}\left(\mathsf{\psi }\right)\frac{\partial \mathsf{\psi }}{\partial \mathrm{x}}\right]+\frac{\partial }{\partial \mathrm{z}}\left[\mathrm{K}\left(\mathsf{\psi }\right)\left(\frac{\partial \mathsf{\psi }}{\partial \mathrm{z}}-1\right)\right]$$

(2)

This equation can be solved on a solution domain as shown in Figure 1.

The function to define solution domain geometry on its upper border is [17]:

$$\mathrm{z}=\frac{{\mathrm{P}}_{\mathrm{s}}}{2}\left[1-\cos \left(\mathsf{\pi }\frac{2\mathrm{x}}{\mathrm{L}}\right)\right]$$

(3)

Spatial pressure distribution as the initial condition necessary in obtaining two-dimensional Richards’ equation solution is:

$${\mathsf{\psi }=\mathsf{\psi }}_{\mathrm{o}}\left(\mathrm{x},\mathrm{z}\right)$$

(4)

Boundary conditions are: $\overline{\mathrm{EF}}$ Dirichlet frontier, prescribed potential according to the water depth in furrow center, $\overline{\mathrm{AB}}$, $\overline{\mathrm{CD}}$, $\overline{\mathrm{DE}}$ and $\overline{\mathrm{FA}}$ Neumann frontiers with zero flux, and $\overline{\mathrm{BC}}$ frontier with unitary gradient:

$$\mathsf{\psi }=\mathrm{h}-\left({\mathrm{z}}_{\mathrm{L}/2} -\mathrm{z}\right),\mathrm{x}\in \overline{\mathrm{EF}},\mathrm{z}\in \overline{\mathrm{AB}},\mathrm{t}>0$$

(5)

$$-\mathrm{K}\left(\mathsf{\psi }\right)\frac{\partial \left(\mathsf{\psi }-\mathrm{z}\right)}{\partial \mathrm{z}}=0,\mathrm{x}\in \overline{\mathrm{FA}},\mathrm{z}\in \overline{\mathrm{AB}},\mathrm{t}>0$$

(6)

$$-\mathrm{K}\left(\mathsf{\psi }\right)\frac{\partial \left(\mathsf{\psi }-\mathrm{z}\right)}{\partial \mathrm{z}}=0,\mathrm{x}\in \overline{\mathrm{DE}},\mathrm{z}\in \overline{\mathrm{CD}},\mathrm{t}>0$$

(7)

$$-\mathrm{K}\left(\mathsf{\psi }\right)\frac{\partial \left(\mathsf{\psi }\right)}{\partial \mathrm{x}}=0,\mathrm{x}=0,\mathrm{z}\in \overline{\mathrm{AB}},\mathrm{t}>0$$

(8)

$$-\mathrm{K}\left(\mathsf{\psi }\right)\frac{\partial \left(\mathsf{\psi }\right)}{\partial \mathrm{x}}=0,\mathrm{x}=\mathrm{L},\mathrm{z}\in \overline{\mathrm{CD}},\mathrm{t}>0$$

(9)

$$\frac{\partial \left(\mathsf{\psi }-\mathrm{z}\right)}{\partial \mathrm{z}}=-1,\mathrm{x}\in \overline{\mathrm{BC}},\mathrm{z}=\mathrm{P},\mathrm{t}>0$$

(10)

Richards’ equation solution is required to represent the hydrodynamic properties of soil, expressing the pressure potential $\left(\mathsf{\psi }\right)$ as a function of water content $\left(\mathsf{\theta }\right)$ and hydraulic conductivity $\left(\mathrm{K}\right)$ as a function of $\mathsf{\theta }$.

Simulation of water flow requires a soil hydraulic characterization, as noted by Fuentes et al. [18], with a combination of characteristics from Fujita [19] and Parlange et al. [20] used in theoretical studies. Experimental studies may require the combination of the retention curve proposed by van Genuchten [21], with the restriction by Burdine [22], and the hydraulic conductivity curve proposed by Brooks and Corey [23] as they meet the integral properties of infiltration, and facilitate the identification of their parameters.

The retention curve is [21]:

$$\frac{\mathsf{\theta }\left(\mathsf{\psi }\right){-\mathsf{\theta }}_{\mathrm{r}}}{{\mathsf{\theta }}_{\mathrm{s}}{-\mathsf{\theta }}_{\mathrm{r}}}= {\left[1+ {\left(\frac{\mathsf{\psi }}{{\mathsf{\psi }}_{\mathrm{d}}}\right)}^{\mathrm{n}}\right]}^{-\mathrm{m}}$$

(11)

where ${\mathsf{\psi }}_{\mathrm{d}}$ is the water pressure characteristic value into the soil; m and n are empirical parameters related by Burdine’s restriction [22]: m = 1 – 2/n, with $0<\mathrm{m}<1$ and $\mathrm{n}>2$; ${\mathsf{\theta }}_{\mathrm{s}}$ is the water content at effective saturation of soil, and ${\mathsf{\theta }}_{\mathrm{r}}$ is the residual water content.

The hydraulic conductivity curve is [23]:

$$\mathrm{K}\left(\mathsf{\theta }\right){=\mathrm{K}}_{\mathrm{s}}{\left(\frac{{\mathsf{\theta }-\mathsf{\theta }}_{\mathrm{r}}}{{\mathsf{\theta }}_{\mathrm{s}}{-\mathsf{\theta }}_{\mathrm{r}}}\right)}^{\mathsf{\eta }}$$

(12)

where $\mathsf{\eta }$ is an empirically positive parameter. Cumulative infiltration per unit furrow length is:

$${\mathrm{A}}_{\mathrm{I}}\left(\mathrm{x},\mathrm{t}\right)={{\displaystyle \int }}_{{\mathrm{P}}_{\mathrm{m}}}{\mathrm{q}}_{\mathrm{I}}\left(\mathrm{x},\mathrm{y},\mathrm{z},\mathrm{t}\right){\mathrm{dP}}_{\mathrm{m}}$$

(13)

where ${\mathrm{q}}_{\mathrm{I}}\left(\mathrm{x},\mathrm{y},\mathrm{z},\mathrm{t}\right)$ is the infiltration rate per unit at the furrow’s surface and ${\mathrm{P}}_{\mathrm{m}}$ is the wet perimeter of the furrow.

2.1. Numerical Solution with Finite Element Method

The numerical two-dimensional Richards’ equation solution used is made of finite elements in space, and implicit finite differences in time. The solution of the algebraic equations system that results from the application of the finite element method is effected using the preconditioned conjugate gradient method [24], which has been adapted to use a free from zeros computational storage matrix.

Multiplying Equation (2) by a weight function $\left(\mathsf{\nu }\right)$, the weak expression form of the two-dimensional Richards’ equation is obtained through integration using Green’s theorem on solution domain $\left(\mathrm{R}\right)$ limited by the border $(\Re )$:

$${{\displaystyle \int }}_{\mathrm{R}}\mathrm{C}\left(\mathsf{\psi }\right)\frac{\partial \mathsf{\psi }}{\partial \mathrm{t}}\mathsf{\nu }\mathrm{dR}+{{\displaystyle \int }}_{\mathrm{R}}\mathrm{K}\left(\mathsf{\psi }\right)\left[\frac{\partial \mathsf{\psi }}{\partial \mathrm{z}}\frac{\partial \mathsf{\nu }}{\partial \mathrm{z}}+\frac{\partial \mathsf{\psi }}{\partial \mathrm{x}}\frac{\partial \mathsf{\nu }}{\partial \mathrm{x}}\right]\mathrm{dR}={{\displaystyle \int }}_{\mathrm{R}}\mathrm{K}\left(\mathsf{\psi }\right)\frac{\partial \mathsf{\nu }}{\partial \mathrm{z}}\mathrm{dR}+{{\displaystyle \int }}_{\mathsf{\Gamma }}\mathrm{q}\mathsf{\nu }\mathrm{d}\mathsf{\Gamma }$$

(14)

where $\mathsf{\Gamma }$ is a fraction of $\Re$ under the Neumann’s condition with prescribed flux $\mathrm{q}$. The Equation (2) solution is assumed as a linear combination of base functions ${\mathsf{\phi }}_{\mathrm{i}}\left(\mathrm{x},\mathrm{z}\right)$ defined in relation to the Kronecker’s delta function and applied to each particular node:

$${\hat{\mathsf{\psi }}}_{\mathrm{n}}(\mathrm{x},\mathrm{z},\mathrm{t})={\displaystyle \sum }_{\mathrm{j}=1}^{\mathrm{n}}{\mathrm{a}}_{\mathrm{j}}(\mathrm{t}){\mathsf{\phi }}_{\mathrm{j}}(\mathrm{x},\mathrm{z})$$

(15)

where ${\mathrm{a}}_{\mathrm{j}}$(t) are coefficients to be determined and n is the number of nodes where the finite element solution is obtained. It is replaced in the first weak form of Richards’ equation (Equation (14)), considering the following: (i) Weight functions are considered equal to the base functions $\left(\mathsf{\phi }\right)$ corresponding to the interior nodes. (ii) A linear variation of the hydrodynamic characteristics is assumed on each element expressing it through the functions of form, i.e., $\hat{\mathrm{C}}={\mathsf{\phi }}_{\mathrm{g}}{\mathrm{C}}_{\mathrm{g}}$ and $\hat{\mathrm{K}}={\mathsf{\phi }}_{\mathrm{g}}{\mathrm{K}}_{\mathrm{g}}$. (iii) A concentrated mass system is used in order to obtain a diagonal matrix and to improve the stability of the scheme [25,26]. Temporal derivative is approximated by an implicit scheme in finite differences and is obtained:

$$\left[\frac{{\mathbf{M}}^{\mathrm{k}+1}}{\Delta \mathrm{t}}{+\mathbf{K}}^{\mathrm{k}+1}\right]{\mathbf{a}}^{\mathrm{k}+1}=\left[\frac{{\mathbf{M}}^{\mathrm{k}+1}}{\Delta \mathrm{t}}\right]{\mathrm{a}}^{\mathrm{k}}{+\mathbf{B}}^{\mathrm{k}+1}{+\mathbf{Q}}^{\mathrm{k}+1}$$

(16)

where the matrices are calculated as indicated below when using linear base functions:

$${\mathrm{M}}_{\mathrm{kj}}={\displaystyle \sum }_{\mathrm{j}=1}^{\mathrm{n}}\left[{\mathrm{C}}_{\mathrm{g}}{{\displaystyle \int }}_{\mathrm{R}}{\overline{\mathsf{\phi }}}_{\mathrm{g}}{\overline{\mathsf{\phi }}}_{\mathrm{j}}{\overline{\mathsf{\phi }}}_{\mathrm{k}}\mathrm{dR}\right]={\displaystyle \sum }_{\mathrm{e}}{\mathrm{C}}_{\mathrm{j}}\frac{\Delta }{3}$$

(17)

$${\mathrm{K}}_{\mathrm{kj}}={\displaystyle \sum }_{\mathrm{j}=1}^{\mathrm{n}}\left[{\mathrm{K}}_{\mathrm{g}}\underset{\mathrm{R}}{{\displaystyle \int }}{\mathsf{\phi }}_{\mathrm{g}}\left(\frac{\partial {\mathsf{\phi }}_{\mathrm{j}}}{\partial \mathrm{x}}\frac{\partial {\mathsf{\phi }}_{\mathrm{k}}}{\partial \mathrm{x}}+\frac{\partial {\mathsf{\phi }}_{\mathrm{j}}}{\partial \mathrm{z}}\frac{\partial {\mathsf{\phi }}_{\mathrm{k}}}{\partial \mathrm{z}}\mathrm{dR}\right)\right]={\displaystyle \sum }_{\mathrm{e}}\frac{\overline{\mathrm{K}}}{4\Delta }({\mathrm{m}}_{\mathrm{j}}{\mathrm{m}}_{\mathrm{k}}+{\mathrm{p}}_{\mathrm{j}}{\mathrm{p}}_{\mathrm{k}})$$

(18)

$${\mathrm{B}}_{\mathrm{k}}{=\mathrm{K}}_{\mathrm{g}}\underset{\mathrm{R}}{{\displaystyle \int }}{\mathsf{\phi }}_{\mathrm{g}}\frac{\partial {\mathsf{\phi }}_{\mathrm{k}}}{\partial \mathrm{z}}\mathrm{dR}={\displaystyle \sum }_{\mathrm{e}}\frac{\overline{\mathrm{K}}}{2}{\mathrm{p}}_{\mathrm{k}}$$

(19)

$${\mathrm{Q}}_{\mathrm{k}}=\underset{\Gamma }{{\displaystyle \int }}\mathrm{q}{\overline{\mathsf{\phi }}}_{\mathrm{k}}\mathrm{d}\Gamma ={\displaystyle \sum }_{\mathrm{e}}\frac{\mathrm{q}{\mathrm{L}}_{\mathrm{j}}}{2}$$

(20)

In previous equations, $\overline{\mathsf{\phi }}$ represents the so-called concentrated mass functions defined as unitary functions on a barycentric region corresponding to a specific node and zero in the rest of the domain. $\mathsf{\Delta }$ is the element area, $\overline{\mathrm{k}}$ is the element conductivity calculated as the arithmetic average of the conductivities obtained in each of its corners (as a consequence of the form adopted for the base functions), ${\mathrm{C}}_{\mathrm{j}}$ is the specific capacity estimated at the node j, ${\mathrm{L}}_{\mathrm{j}}$ is the border length corresponding to each node under Neumann’s condition, m and p are geometric factors defined according to the base functions: ${\mathrm{m}}_{\mathrm{i}}{=\mathrm{z}}_{\mathrm{j}}{-\mathrm{z}}_{\mathrm{k}}$ and ${\mathrm{p}}_{\mathrm{i}}{=\mathrm{x}}_{\mathrm{j}}{-\mathrm{x}}_{\mathrm{k}}$ where the subscripts i, j and k correspond to the corners of the triangular element and run on their three possible sequenced permutations.

2.2. Solution Domain Characteristics

The cross section of the furrow had the following dimensions: a width of 100 cm, and a depth of 150 cm; in the upper zone of the domain, the shape of the furrow obtained with Equation (3) was introduced.

To solve two-dimensional Richards’ equation numerically five finite element meshes with different densities were built (Figure 2). Characteristics of each mesh are shown in

Table 1. Characteristics of finite element meshes.

Mesh	Elements	Nodes	$\mathbf{Maximum}\mathbf{Increase}\mathsf{\Delta }\mathbf{x}\left(\mathbf{cm}\right)$	$\mathbf{Maximum}\mathbf{Increase}\mathsf{\Delta }\mathbf{z}\left(\mathbf{cm}\right)$
1	575	334	10.0	10.0
2	1401	764	5.0	5.0
3	8122	4197	2.0	2.0
4	32,014	16,272	1.0	1.0
5	127,732	64,387	0.5	0.5

.

Time steps discretization proposed are $\Delta \mathrm{t}=1.0$, $\Delta \mathrm{t}=5.0$, $\Delta \mathrm{t}=10.0$, $\Delta \mathrm{t}=30.0$ and $\Delta \mathrm{t}=60.0$ s.

Selection of three different soil textures was made to cover a spectrum that involved different conditions. The contrasting chosen soil textures ensure that other soil textures have intermediate behaviors among the analyzed ones. Hydrodynamic characteristics of the three contrasting soil types are shown in

Table 2. Hydrodynamic characteristics for three soil textures.

Soil	${\mathsf{\theta }}_{\mathbf{r}}({\mathbf{cm}}^3/{\mathbf{cm}}^3)$	${\mathsf{\theta }}_{\mathbf{s}}({\mathbf{cm}}^3/{\mathbf{cm}}^3)$	${\mathsf{\psi }}_{\mathbf{d}}\left(\mathbf{cm}\right)$	${\mathbf{K}}_{\mathbf{s}}(\mathbf{m}/\mathbf{s})$	$\mathbf{m}=1-2/\mathbf{n}$	$\mathsf{\eta }$
Sandy loam	0.0	0.450	9.52	1.390 × 10−4	0.1004	13.62
Silt loam	0.0	0.525	29.35	0.0167 × 10−4	0.1165	12.01
Clay loam	0.0	0.475	34.15	0.0042 × 10−4	0.0714	19.30

. Residual saturation was assumed equal to zero in accordance with Fuentes et al. [18]. Total soil porosity was assumed as saturation water content $\left(\mathsf{\varphi }\right)$, determined based on values provided by Rawls et al. [27]. For determined values of m and $\mathsf{\eta }$ parameters, a granulometric curve was built for each soil, based in sand, silt and clay percentages of triangle textures of USDA (U.S. Department of Agriculture), following Fuentes’ procedure [18]. Water content necessary to assign Richards’ equation initial condition was determined considering usable humidity of each soil type, assuming 50% humidity before irrigation. Usable humidity was determined by subtracting water content corresponding to field capacity and permanent wilting point.

2.3. Border Condition in the Furrow

To solve two-dimensional Richards’ equation numerically, it is necessary to know time–water depth evolution in furrow. Accepting the contact time hypothesis [28] allows proposing a unique form of time–water depth evolution in furrow throughout its length as a first approximation. The internal numerical coupling of the Saint-Venant and Richards’ equations for border irrigation was done by using a Eulerian–Lagrangian method based on [29], where the water depth over the soil was obtained by the numerical solution of the full Saint-Venant equations for mass and momentum conservation, while the infiltrated water depth was calculated by the finite element solution of the two-dimensional Richards’ equation. For a given border length, the Christiansen Uniformity Coefficient (CUC) was evaluated for a varied water intake flow. It was found that the CUC was very sensitive to the water intake flow; the water intake flow that corresponded to a maximum CUC was denominated optimal intake flow. The process was repeated for several lengths of a border and was found that a basically lineal relationship exist between the optimal intake flow and the border length. This last process was done for ten soil types of the texture triangle of USDA.

Table 3. Indicators to obtain time–water depth evolution form, reported by Saucedo et al. [16].

	Soil
	Sandy Loam	Silt Loam	Clay Loam
Qopt (lps/m)	0.02476	0.00446	0.00088
Irrigation time (h)	1.2	6.2	31.4
Normal water depth (border irrigation) (cm)	2.42	0.98	0.50
Normal water depth (furrow irrigation) (cm)	8.92	4.99	3.24

presents these values relative to three soil types, including optimal intake flow and the corresponding irrigation time. If use is made of a hydraulic resistance flow law of potential form [30], then the normal water depth for a border is obtained, which is transformed to an equivalent normal water depth for a furrow. Figure 3a–c exposes the evolution of equivalent furrow water depth at the beginning of the furrow, scaled from the evolution of border water depth at the beginning of the border obtained by the application of the described numerical model.

Through the internal numerical coupling of the Saint-Venant and Richards’ equations in border irrigation, Saucedo et al. [16] report values of optimal flow and corresponding irrigation time, for a slope of 0.002 m/m, indicated in Table 3. Considering normal water depth (perpendicular to the soil surface), a per unitary border width is obtained for a furrow corresponding to water depth.

Stability criteria was assumed per Banti et al. [8], and therefore a discretization was made in a space smaller than that used in the referred validation process, and different time steps mostly less than or equal to 30.0 s, so that the relative error maximum level was similar to that estimated in the referenced work.

3. Results and Discussion

The two-dimensional Richards’ equation solution was obtained applying finite element method combined with a temporal approximation based on implicit finite differences. The above was done by computational programming in C++ language. The model was applied in the solution domain to calculate the infiltrated depth, considering three contrasting soil textures. Results obtained considering different spatial-temporal discretizations were compared. Five meshes of finite element and five time steps were used. Results were evaluated using the root mean squared error (RMSE) [9]:

$$\mathrm{RMSE}=\sqrt{\frac{1}{\mathrm{n}}{\displaystyle \sum }_{\mathrm{i}=1}^{\mathrm{n}}{\left({\mathrm{M}}_{\mathrm{i}}{-\mathrm{m}}_{\mathrm{i}}\right)}^2}$$

(21)

where n is the number of infiltrated depth data pairs, ${\mathrm{M}}_{\mathrm{i}}$ is the infiltrated depth obtained using the denser finite element mesh for a given time and ${\mathrm{m}}_{\mathrm{i}}$ is the infiltrated depth obtained with a different mesh from ${\mathrm{M}}_{\mathrm{i}}$ for the same time step.

Figure 4, Figure 5 and Figure 6 show the results of infiltrated depth for three different soil textures, with two denser finite element meshes and a time step of 1.0 s.

Figure 7, Figure 8 and Figure 9 present the volumetric soil water content (%). Typical irrigation times were used, and results presented consider a time step of 1.0 s, with spatial discretization according to a denser finite element mesh.

Figure 10a–c shows the values of the RMSE. These correspond to the results of infiltrated depth obtained with mesh 5 combined with meshes 1, 2, 3 and 4. With that combination the maximum error results; any other combination of meshes to determine the RMSE had a lower value than those presented below.

As presented above, results were compared considering five time steps as reference. Figure 11 shows the values of the RMSE indicator that compare the infiltrated depth evolution results, considering as standard the densest mesh combined with the shortest time step.

Instability of infiltrated depth evolution results for sandy loam soil were evident when time steps were greater than 10.0 s. This result is similar to studies reported in the literature [31,32].

Another aim was to compare results shown in Figure 3, Figure 4 and Figure 5 that were obtained when considering the infiltration model relative to the normal water depth indicated in Table 3 (constant water depth in time). The result is shown in Figure 12, Figure 13 and Figure 14. Mesh 5 and a time step of 1.0 s were used.

Relevant results are shown in Figure 12, Figure 13 and Figure 14. They make evident the convenience of considering time–water depth evolution in infiltration simulations. In fact, values of RMSE indicator shown in Figure 12, Figure 13 and Figure 14 are greater than values showed in Figure 10a–c.

To expand the scope of this study and for a fuller exploration of the subject, the results were compared with results obtained by applying the HYDRUS-2D program [33] (Version 2.0, PC-Progress, Prague, Czech Republic). The HYDRUS program numerically solves the two-dimensional Richards’ equation for saturated–unsaturated water flow. HYDRUS implements the soil-hydraulic functions of [21], using the statistical pore-size distribution model of [34] to obtain a predictive equation for the unsaturated hydraulic conductivity function in terms of soil water retention parameters. The expressions of [21] are given as follows:

$$\mathsf{\theta }\left(\mathsf{\psi }\right)=\{\begin{array}{cc}{\mathsf{\theta }}_{\mathrm{r}}+\frac{{\mathsf{\theta }}_{\mathrm{s}}-{\mathsf{\theta }}_{\mathrm{r}}}{{\left[1+{\left|\mathsf{\alpha }\mathsf{\psi }\right|}^{\mathrm{n}}\right]}^{\mathrm{m}}}& \mathsf{\psi }<0\\ {\mathsf{\theta }}_{\mathrm{s}}& \mathsf{\psi }\ge 0\end{array}$$

(22)

$$\mathrm{K}\left(\mathsf{\psi }\right){=\mathrm{K}}_{\mathrm{s}}{\mathrm{S}}_{\mathrm{e}}^{\mathrm{l}}{\left[1-{\left(1-{\mathrm{S}}_{\mathrm{e}}^{1/\mathrm{m}}\right)}^{\mathrm{m}}\right]}^2$$

(23)

where

$$\mathrm{m}=1-1/\mathrm{n},\mathrm{n}>1$$

(24)

The above equations contain six independent parameters: ${\mathsf{\theta }}_{\mathrm{r}}$, ${\mathsf{\theta }}_{\mathrm{s}}$, $\mathsf{\alpha }$, $\mathrm{n}$, ${\mathrm{K}}_{\mathrm{s}}$ and $\mathrm{l}$. The pore-connectivity parameter l in the hydraulic conductivity function was estimated [34] to be about 0.5 as an average for many soils. ${\mathrm{S}}_{\mathrm{e}}$ is the degree of saturation.

Results shown in Figure 12, Figure 13 and Figure 14 were compared with those obtained applying the normal water depth indicated in Table 3 (constant water depth in time) in the HYDRUS-2D program. Parameters used to implement the soil-hydraulic functions of van Genuchten–Mualem to obtain a predictive equation for unsaturated hydraulic conductivity function are shown in

Table 4. Hydrodynamic characteristics for three soil textures used in HYDRUS-2D.

Soil	${\mathsf{\theta }}_{\mathbf{s}}({\mathbf{cm}}^3/{\mathbf{cm}}^3)$	${\mathbf{K}}_{\mathbf{s}}(\mathbf{m}/\mathbf{s})$	$\mathbf{m}=1-1/\mathbf{n}$	n (n > 1)	l
Sandy loam	0.450	1.390 × 10−4	0.1004	1.1116	0.5
Silt loam	0.525	0.0167 × 10−4	0.1165	1.1318	0.5
Clay loam	0.475	0.0042 × 10−4	0.0714	1.0769	0.5

.

Results are shown in Figure 15, Figure 16 and Figure 17. A similar mesh to mesh 5 is used in the HYDRUS-2D program, as well as a time step of 1.0 s.

From the results shown in Figure 15, Figure 16 and Figure 17 it can be observed that the difference obtained with the HYDRUS-2D program and the model presented in this study were due to the different models used to obtain the humidity retention curve. It is possible to match the humidity retention curves of both models by adjusting the “m” parameter, but this is not enough to match the hydraulic conductivity functions. In addition, the difference between both models was due to spatial-temporal discretization, and the enforcement of a variable and constant water depth in each case.

4. Conclusions

By obtaining the time infiltrated depth evolution and humidity profiles, with a numerical solution of the two-dimensional Richards’ equation, the water infiltration during two-dimensional infiltration was simulated. The contact time hypothesis has been used to apply a unique form on time evolution of the water depth in the furrow throughout its length, as boundary condition. The specific form of such evolution in time was obtained from results reported in the literature based on the internal numerical coupling of the Saint-Venant and Richards’ equations in border irrigation. Moreover, the equivalent hydraulic area between the border and the furrow was achieved by scaling the values of water depth. The model was evaluated using three contrasting soil types and the comparison was done by computing the RMSE indicator (root mean square error). The comparison was performed from the selection of five finite element meshes with different densities, to discretize the solution domain of the two-dimensional Richards’ equation, combined with several time steps. Finally, a comparison was made between infiltrated depth evolution calculated with a constant water depth in the furrow to the one proposed in this work, finding important differences between both approaches. To expand the scope of this study and for a fuller exploration of the subject, the results are compared with results obtained by applying the HYDRUS-2D program. The results confirm that it is important to consider an internal full coupling of the Saint-Venant and Richards’ equations to improve furrow irrigation simulations. In other words, future efforts should be directed to the internal numerical full coupling of the Saint-Venant and two-dimensional Richards’ equations to adequately describe furrow irrigation.