Modeling Soil Water Redistribution under Gravity Irrigation with the Richards Equation

Soil water movement is important in fields such as soil mechanics, irrigation, drainage, hydrology, and agriculture. The Richards equation describes the flow of water in an unsaturated porous medium, and analytical solutions exist only for simplified cases. However, numerous practical situations require a numerical solution (1D, 2D, or 3D) depending on the complexity of the studied problem. In this paper, numerical solution of the equation describing water infiltration into soil using the finite difference method is studied. The finite difference solution is made via iterative schemes of local balance, including explicit, implicit, and intermediate methods; as a special case, the Laasonen method is shown. The found solution is applied to water transfer problems during and after gravity irrigation to observe phenomena of infiltration, evaporation, transpiration, and percolation.


Introduction
Surface gravity irrigation consists of the supply of water at the head of a channel or inclined channel built on a plot, such as a border or a furrow, to take advantage of the gravitational field and provide the necessary amount of water for the development of cultivated plants [1,2].
In surface gravity irrigation, there are three phases-the advanced phase, the storage phase, and the recession phase-which together are studied using a variety of models to understand the phenomenon [3], ranging from empirical to physical-based models that use the Barré de Saint-Venant [4] and Richards equations [5] to model the surface movement and underground movement, respectively [2].
The differential equation that describes the water transfer in soils [5] is difficult to integrate in an analytical way for most initial and boundary conditions of practical interest, due to its highly nonlinear nature [2][3][4][5][6]. This situation leads to the use of numerical methods for their integration and application in simulations of the moisture profile evolution and of the infiltrated depth during gravity irrigation [2,[7][8][9], and also in simulations of water redistribution after irrigation due to the gravitational field, water evaporation through the soil surface, and crop transpiration [7][8][9].
The general equation of water transfer in soil, 2 of 13 results from the combination of the continuity equation which results from the principle of conservation of mass, ∂θ ∂t + ∇ · q + Υ = 0, (2) and from Darcy's Law [10], argued experimentally, which represents the dynamic equation resulting from the principle of the amount of momentum: where θ is the volumetric water content, also known as the moisture content; q = (q x , q y , q z ) is the flow of water per unit area of soil or Darcy flow, with its components in the directions of space (x, y, z); t is the time; Υ is a term of source or sink as the volume of extracted water by plants per unit volume of soil and unit of time; H is the hydraulic potential or hydraulic load, being equal to the sum of pressure potential (ψ), which is negative in the unsaturated area of soil and positive in the saturated area, and the gravitational potential assimilated to the vertical coordinate positively oriented downward (z); and K is the hydraulic conductivity of the partially saturated soil and depends on pressure or moisture content. The relationships θ(ψ) and K(ψ) are called soil hydrodynamic characteristics. The differential equation contains two apparently dependent variables, namely, the moisture content and water pressure, but these are related through the retention curve. With the introduction in the equation of the specific capacity, defined as the slope of the curve retention C(ψ) = dθ/dψ, that is, the pressure only appears to be a dependent variable; this is known as the Richards equation [5].
As for the analytical representation of the hydrodynamic characteristics, according to [11], the combination of the hydrodynamic characteristics of Fujita [12] and Parlange [13] is used for theoretical studies, such as the construction of exact analytical solutions; in experimental studies, a combination of the retention curve proposed by van Genuchten [14] is used, considering the restriction of Burdine [15], with the curve of hydraulic conductivity proposed by Brooks and Corey [16], due to the fact that they satisfy the integral properties of infiltration and easy identification of its parameters.
In this study, we focus on solving the Richards vertical one-dimensional differential equation: which is a combination of the continuity equation, and Darcy's Law,

Solution in Finite Difference by Local Balance of Mass
To build the numerical solution, a scheme proposed by Zataráin et al. [17] was used. According to the algorithm, the soil profile is divided into layers that may be non-constant, as shown in Figure 1. The interpolation parameters are introduced and defined by The pressure ψ in an intermediate node i + γ z for all k is estimated via The pressure ψ k i−1+γ z is deduced by replacing i with i−1 in Equation (9), and the pressure at the intermediate time k + ω for all i is estimated as The time derivative of Equation (6) for the ith node is calculated as follows: where C k+ω i = C ψ k+ω i . The discretization of the spatial derivative in Equation (6) is According to Equation (7), Thus, the divergence defined by Equation (12) is as follows: The dependent variable involved in the hydraulic conductivity of Equation (13) is defined as The pressure ψ k+ω i−1+γ z is deduced replacing i with i−1 in Equation (15). Substituting Equation (10) into Equation (14) and this, in turn, into Equation (6) together with Equation (11) leads to the following system of algebraic equations: where It is important to note that the system is not linear since the coefficients (17)-(20) depend on the solution itself. Therefore, the resolution for each time step is iterative. The first estimator can be obtained by the Laasonen scheme, which is the result of considering that the coefficients, specific capacity, and hydraulic conductivity are calculated in the previous time step while the hydraulic gradient is calculated at the current time. The method is partially implicit, and the resulting algorithm is no longer iterative in a time step. The coefficients of the system of equations defined by Equation (16) are defined by the following expressions deduced from Equations (17)- (20): Mathematics 2020, 8, 1581 5 of 13 The system (16) forms a tridiagonal matrix and can be solved efficiently by an algorithm of Gaussian elimination known as the Thomas algorithm or double scan. The solution of system (16) has the form ψ k+1 can be deduced, which, substituted into Equation (16), leads to Equation (25) provided that Equations (26) and (27) require knowledge of F 1 and G 1 , which are deduced from the boundary condition in z = 0; this is the forward scan, i = 2, 3, . . . , n−1. Subsequently, ψ k+1 n is defined from the boundary condition in z = L, and ψ k+1 i is calculated for i = n−1, n−2, . . . , 1; this is the backward scan. If the coefficients in system (16) depend on the pressure itself, they are updated and the algorithm reapplied until the next estimator of the solution t k+1 is equal, given an error criterion, to the previous estimator. The convergence can be accelerated if second-order algorithms like Newton's method are used to find the pressure in a given time step.

Boundary Conditions
The boundary conditions can be grouped into three types: Dirichlet (first type), Neumann (second type), and Robbins (third type). In gravity irrigation (basin irrigation, border irrigation, or furrow irrigation), the Dirichlet condition on the upper boundary (water depth over the soil surface) is one of the main components, and its characterization is fundamental for irrigation design [2,[6][7][8][9]. In this condition, the value of the dependent variable of the differential equation for all time on the soil surface (pressure) is known. It is introduced into the numerical solution by defining F 1 and G 1 from Equation (25). For the upper boundary, I = 1 is used in Equation (25), resulting in the following expression: Due to the fact that ψ 1 is the water pressure on the soil surface and is independent of ψ 2 , we have For the lower boundary, the length of the column must be taken into consideration. In a semi-infinite column, pressure corresponds to the initial condition. In a finite column, pressure is introduced into the base position of the column by using i = n − 1 in Equation (25): With this expression, the backward scan is performed.

The Analytical Solution of Infiltration
Equation (1) can be written with the moisture content as the dependent variable by introducing hydraulic diffusivity, defined as The structure of the resulting equation corresponds to a nonlinear Fokker-Planck equation. The equation for one-dimensional vertical infiltration without a sink term is This equation defines the function θ = θ(z,t) for each type of initial and boundary conditions. Solutions can be built more or less simply considering that there is an implicit function z = z(θ,t), as in the case of infiltration in a semi-infinite column with a constant initial moisture content (θ 0 ) and constant humidity condition in the upper part of the column, generally equal to the moisture content at saturation (θ s ). From these considerations, the total differentials dθ = (∂θ/∂z) t dz + (∂θ/∂t) z dt and dz = (∂z/∂θ) t dθ + (∂z/∂t) θ are deduced, where the subscripts are the variables that remain constant in the differentiation. Substituting the first into the second, we have [(∂z/∂θ) t (∂θ/∂t) z + ((∂z/∂t)] dt = 0. The identity (∂z/∂θ) t (∂θ/∂t) z = −(∂z/∂t) θ is then deduced, which allows us to write Equation (32) as follows: A first integration with the initial condition provides For the second integration, the relationship of the concentration of flows is introduced and defined as where q s (t) = q z (θ s , t) is the Darcy flux on the surface of the column or infiltration velocity. From Equations (34) and (35), the moisture profile is deduced: in which the boundary condition z(θ s , t) = 0 has been used. The moisture profile integration provides water storage in the soil column: where I(t) is the water depth accumulated over time, with q s (t) = dI/dt. Equations (36) and (37) can be integrated in an approximated manner using the hydrodynamic characteristics of Fujita [12] and Parlange [13]. Fujita diffusivity is provided by where α is a shape parameter; K s the hydraulic conductivity at saturation, and λ c is the Bouwer scale [18] defined by and the degree of saturation, θ * , is defined by The introduction of Fujita diffusivity into the relationship between hydraulic conductivity and diffusivity proposed by Parlange et al. [13] gives the following expression for the hydraulic conductivity: where β is a shape parameter.
It should be noted that the Fujita diffusivity contains extreme behaviors, since with α = 0, we have constant diffusivity, and when α→1, the diffusivity is assumed to be a Dirac density; the behavior of soils is nearer the latter than the former. In the case of a Dirac density, the relationship of the concentration of flows is equal to the degree of saturation defined by Equation (40).
The introduction of Equations (38) and (41) into Equations (36) and (37), considering that F(θ,t)→θ* when α→1, leads to the following solutions of the infiltration phenomenon: where the dimensionless variables are defined by The parameter S is the sorptivity that for a Dirac medium (α→1) is defined by S 2 = 2 (K s − K 0 ) λ c (θ s − θ 0 ). Equation (43) is a differential equation since q * s = dI*/dt*. Integration with the condition I = 0 in t = 0 leads to the infiltration formula implicitly describing the evolution of the infiltrated water depth as a function of time, firstly obtained by Parlange et al. [13]: The retention curve of Fujita and Parlange is obtained from the definition of hydraulic diffusivity, Equation (31), considering the condition ψ(θ s ) = ψ s , where ψ s is a bubble pressure in saturation that can be taken equal to zero: For the case in which θ 0 = θ r , where θ r is the residual moisture content defined in such a way that K (θ r ) = 0 and ψ(θ r )→−∞, that is, k = 0, the next curve is deduced: where is the effective degree of saturation.

Comparison between Numerical and Analytical Solutions
A comparison between the numerical solution and analytical solution was carried out in the soil of Montecillo, State of Mexico, the American textural classification of which is a sandy loam soil; the soil was characterized hydrodynamically via the internal drainage method [19]. According to Fuentes et al. [20], the values of the parameters of the Fujita and Parlange characteristics are θ s = 0.520 cm 3 /cm 3 , θ 0 = 0.185 cm 3 /cm 3 , ψ c = 13.500 cm, ψ s = 0.000 cm, K s = 2.5 cm/h, K 0 = 0.000 cm/h, α = 0.969, and β = 0.998.
Of the behavior in short time periods of the infiltration solution under the conditions studied, a relationship between space and time steps can be deduced. In short time periods, I = S √ t, so in the first time step, I ≈ ∆z∆θ and t = ∆t, where ∆θ = θ s − θ 0 ; since S 2 ≈ 2 (K s − K 0 ) λ c ∆θ, the time and space steps must be such that the dimentionless parameter is of the order of unity [17].
The comparison between the numerical solution of finite differences of the Richards equation and the analytical solution of Parlange [13] is shown in Figure 2. The results are similar if the parameter interpolator in time (ω) is close to unity, that is, if the scheme is fully implicit. The convergence of the numerical scheme and its stability are discussed in depth elsewhere [21].

Simulations under Different Boundary Conditions
Numerical simulations of infiltration, redistribution, and evaporation of water in a soil column in Montecillo, State of Mexico were performed using the hydrodynamic characteristics of van Genuchten [14] for the retention curve and those of Brooks and Corey [16] for the hydraulic

Simulations under Different Boundary Conditions
Numerical simulations of infiltration, redistribution, and evaporation of water in a soil column in Montecillo, State of Mexico were performed using the hydrodynamic characteristics of van Genuchten [14] for the retention curve and those of Brooks and Corey [16] for the hydraulic conductivity curve. The retention curve proposed by van Genuchten is written as where θ s is the volumetric water content at saturation, θ r is the volumetric content of residual water, ψ d is a characteristic value of the water pressure in soil, and m and n are two parameters of empirical shape related here by the Burdine restriction [15]: m = 1 − 2/n, with 0 < m < 1 and n > 2. The hydraulic conductivity proposed by Brooks and Corey is where η is a parameter of positive form.
According to the methodology proposed by Fuentes [22], the parameters θ s , θ r , m, and η of Equations (50) and (51) were estimated from the porosity and granulometric curve and parameters Ks and ψ d by applying an inverse parameter estimation method using advance phase data in an irrigation test and the Barré de Saint-Venant [4] and Richards [5] equations, in a border with a sandy loam soil in Montecillo, Mexico. The values were θ s = 0.4865 cm 3 /cm 3 , θ r = 0.0 cm 3 /cm 3 , m = 0.125, n = 2.2857, η = 11.0, ψ d = −32.75 cm, and K s = 1.84 cm/h.
The water content at the permanent wilting point (θ PMP ) and field capacity (θ CC ) were estimated as corresponding to water pressure ψ PMP = −15,300 cm and ψ CC = −340 cm, respectively, as follows: θ PMP = 0.0840 cm 3 /cm 3 and θ CC = 0.2492 cm 3 /cm 3 . The initial moisture content at which soil irrigation must be carried out was calculated as θ 0 = θ PMP + F ap (θ CC − θ PMP ), where F ap is the permissible depletion factor of the crop, for example, F ap = 1/3; using the above values, we obtained θ 0 = 0.1391 cm 3 /cm 3 . The water pressure on the surface was taken to be equal to the average, h = 1.5 cm. The net water depth of irrigation was calculated as (θ CC − θ 0 )/P r , where P r is the thickness of the considered soil profile; for a depth of 70 cm, the net water depth is 7.71 cm, which corresponds to an application efficiency of 83.3% to a gross water depth of 9.25 cm, which is applied in an irrigation time of 86 min. Figure 3 shows the results of the numerical simulations. water pressure on the surface was taken to be equal to the average, h = 1.5 cm. The net water depth of irrigation was calculated as (θCC − θ0)/Pr, where Pr is the thickness of the considered soil profile; for a depth of 70 cm, the net water depth is 7.71 cm, which corresponds to an application efficiency of 83.3% to a gross water depth of 9.25 cm, which is applied in an irrigation time of 86 min. Figure 3 shows the results of the numerical simulations.

Evaporation and Transpiration Processes
Immediately after irrigation, processes of water redistribution due to the gravitational field, water consumption by plants, or transpiration and evaporation through the soil surface are started. For purposes of illustration, the evaporation process was simplified and represented by the following harmonic describing its variation in time:

Evaporation and Transpiration Processes
Immediately after irrigation, processes of water redistribution due to the gravitational field, water consumption by plants, or transpiration and evaporation through the soil surface are started. For purposes of illustration, the evaporation process was simplified and represented by the following harmonic describing its variation in time: where ω is the angular frequency, P = 24 h is the period, e v is the evaporated water depth hourly mean, φ = −ωt M is the phase angle, t M is the time at which the hourly evaporation reaches its maximum value e vM = e v (t M ), and c v = e vM − e v . The minimum value is e vm = e v − c v ; therefore, e vm + e vM = 2e v is satisfied.
To calculate the parameters of the harmonic, it was considered that the average annual water depth evaporated is Ev a = 1,168 mm, which provides average daily water depth Ev d = 3. = −e v (t k + 1). For transpiration, there are different possibilities, such as the models that include root density or laws of type Ohm. For purposes of illustration, the sink term studied by Zataráin et al. [17] was taken, namely, where γ is a dimensionless coefficient reflecting the demand for water by the sink; γ = 0.015 was taken as an example. The Bouwer scale [15] is defined by With the hydrodynamic characteristics defined by Equations (50) and (51), it is calculated as where B(p, q) = Γ(p) Γ(q)/Γ(p + q) is the complete beta function with p > 0 and q > 0, and Γ(x) is the complete Euler gamma function. For the sandy loam soil under study, λ c = 33.95 cm. Within the local balance scheme, the discrete form of Equation (53) is where K k+ω j = K (ψ k+ω j ). The interpolation defined by Equation (15) should be used and similar terms should be grouped as in system (16). Figure 4 shows the evolution of the moisture profile after irrigation, due to redistribution caused by the gravitational field, to evaporation through the soil surface, to transpiration by crops, and to percolation in the base of the control column of the sandy loam soil in Montecillo, Mexico. Figure 5 shows the evolution of the corresponding accumulated water depth as a function of time in the profile of the soil control. One month after irrigation, the average moisture content was θ = 0.2060 cm 3 /cm 3 , which should continue decreasing over time to achieve the average initial moisture content θ 0 = 0.1391 cm 3 /cm 3 in order to apply the second irrigation. surface, to transpiration by crops, and to percolation in the base of the control column of the sandy loam soil in Montecillo, Mexico. Figure 5 shows the evolution of the corresponding accumulated water depth as a function of time in the profile of the soil control. One month after irrigation, the average moisture content was θ = 0.2060 cm 3 /cm 3 , which should continue decreasing over time to achieve the average initial moisture content θ0 = 0.1391 cm 3 /cm 3 in order to apply the second irrigation.

Conclusions
A differential equation describing the transfer of water in soil was integrated analytically for a few initial and boundary conditions of practical interest, such as in a semi-infinite soil column with constant initial moisture content and constant moisture conditions at its upper boundary. In the general case, the integration is done in a numerical manner, which is why the numerical scheme was presented via the finite difference method.
Using the finite difference method, a scheme of local balance in one dimension of space was given. This scheme has the peculiarity that the coefficients involving specific capacity and hydraulic conductivity depend on the desired pressure in a given time step, making the resolution of the system of discrete equations occur via an iterative method. The presented solution was compared with the analytical solution, which permitted us to establish a relationship between the steps of time and space. The solutions are similar if the scheme is taken to be close to or equal to a fully implicit scheme.
Finally, with the numerical solution, we have a tool allowing us to know more about the evolution of the moisture profile and infiltrated water depth during irrigation and the evolution of the moisture profile after irrigation: the redistribution of moisture due to the gravitational field, evaporation of water through the soil surface, transpiration from plants with a simplified model of the sink term, and percolation. The use of the numerical solution helps us to understand the movement of water in an irrigation event, which allows us to give pertinent recommendations for the efficient application of water to a plot.

Conclusions
A differential equation describing the transfer of water in soil was integrated analytically for a few initial and boundary conditions of practical interest, such as in a semi-infinite soil column with constant initial moisture content and constant moisture conditions at its upper boundary. In the general case, the integration is done in a numerical manner, which is why the numerical scheme was presented via the finite difference method.
Using the finite difference method, a scheme of local balance in one dimension of space was given. This scheme has the peculiarity that the coefficients involving specific capacity and hydraulic conductivity depend on the desired pressure in a given time step, making the resolution of the system of discrete equations occur via an iterative method. The presented solution was compared with the analytical solution, which permitted us to establish a relationship between the steps of time and space. The solutions are similar if the scheme is taken to be close to or equal to a fully implicit scheme.
Finally, with the numerical solution, we have a tool allowing us to know more about the evolution of the moisture profile and infiltrated water depth during irrigation and the evolution of the moisture profile after irrigation: the redistribution of moisture due to the gravitational field, evaporation of water through the soil surface, transpiration from plants with a simplified model of the sink term, and percolation. The use of the numerical solution helps us to understand the movement of water in an irrigation event, which allows us to give pertinent recommendations for the efficient application of water to a plot.