text in 
English (pdf)
Article in xml format
Article references
Automatic translation
Send this article by e-mail
Cited by SciELO 
Similars in
SciELO 
Permalink
Introduction
Soil water flow is important in areas such as soil mechanics, irrigation, drainage, hydrology and agriculture (Fuentes et al., 2020). This phenomenon can be described physically with the Richards equation (1931), which produces accurate results if the hydrodynamic properties and boundary conditions are known, but it lacks general analytical solutions, and therefore highly complex numerical methods are usually required for its solution (Damodhara-Rao et al., 2006; Ma et al., 2010; Malek & Peters, 2011). Alternatively, there are physically based approximate models, which are the result of simplifying the initial conditions. In particular, the equation of Green and Ampt (1911) is an alternative to simulate the process of vertical infiltration of water into the soil, which has been used in surface irrigation. (Chávez et al., 2018, 2020; Naghedifar et al., 2020; Saucedo et al., 2015).
The application of the Green and Ampt equation involves soil physical property parameters (Ali & Islam, 2018; Damodhara-Rao et al., 2006); however, there are two key parameters that are not viable to measure experimentally: saturated hydraulic conductivity (Ks) and wetting front suction (hf). The value of Ks can be measured directly in the field using rainfall simulators, infiltration meters and constant or variable load permeameters (Gómez-Tagle et al., 2008), as well as Guelph permeameter. However, these methods need long periods of time to achieve a stable value, which leads to the use of large volumes of water for estimation (Castiglion et al., 2018). This uncertainty in the values of these soil hydrodynamic parameters can be considered an obstacle to optimal irrigation design due to the high degree of difficulty in calculating them, although most of the time it depends on the experience of the modeler. The parameters of the Green and Ampt equation have been optimized by least-squares fitting under rainfall conditions (Chen et al., 2015; Xiang et al., 2016). Chen et al. (2015) showed that these parameters are affected by the duration of precipitation used in their modeling. Currently, different theoretical and empirical approaches such as pedotransfer functions (Saxton & Rawls, 2006; Trejo-Alonso et al., 2020) and artificial neural networks are used to estimate these parameters (Trejo-Alonso et al., 2021). Nevertheless, pedotransfer functions are only applicable in the area where the model was calibrated, and neural networks represent a high computational cost (higher computational time).
Therefore, the objectives of this work were to optimize the parameters Ks and hf of the Green and Ampt equation by using a nonlinear optimization algorithm and to validate the solution to achieve the optimization of parameters as a function of soil textures.
Materials and methods
The Green and Ampt equation
The Green and Ampt equation is established from the continuity equation and Darcy's law (1856) with the following assumptions: a) the initial moisture profile in a soil column is uniform: θ = θo, b) water pressure at the soil surface is hydrostatic: ψ = h ≥ 0 (where h is the water depth on the soil surface), c) there is a well-defined wetting front characterized by a negative pressure: ψ = -hf () 0 (where hf is the wetting front suction) and d) the region between the soil surface and the wetting front is completely saturated (piston flow): θ = θs and K = Ks (where Ks is the saturated hydraulic conductivity; that is, the value of the hydraulic conductivity of Darcy's law corresponding to the volumetric water content at saturation). The resulting ordinary differential equation is as follows:
where I the cumulative infiltration, t is the time, and θs and θo are the initial and saturated moisture contents, respectively.
If the water depth over the surface is considered independent of time, Equation (1) is analytically integrated with the initial condition I = 0 at t = 0, which results in the following:
When the amount of water flowing at infinity is neglected, the infiltrated water depth as assumed by Green and Ampt is equal to the volume per unit area stored in the piston: I(t) = zf(t)Δθ with Δθ = θs - θo, and zf(t) is the piston front (Fuentes et al., 2012).
Parameter optimization
In nonlinear regression models, each observation yi is written as a response function f(xi; β). An important difference in nonlinear regression models is that the number of parameters β of the regression is not directly related to the number of variables xi in the model (Cornejo-Zuniga & Rebolledo-Vega, 2016). To estimate the parameters Ks and hf using an infiltration test, the Levenberg-Marquardt algorithm (Moré, 1978) is applied, which has been a standard technique for nonlinear least squares problems (Chávez et al., 2022; Fuentes, et al., 2022; Šimůnek & Hopmans, 2018). The parameter update is performed iteratively with the following expression:
where J is the Jacobian matrix related to the variations of the infiltration function for each parameter to be optimized, Id is the identity matrix, r is the vector of differences between the measured infiltrated water depth and the one calculated with the optimization algorithm, and μ is the damping parameter that is fitted in each iteration. The Jacobian matrix is calculated as follows:
where Im is the equation of Green and Ampt and m is the number of measured data. In this algorithm the values of the Jacobian matrix are approximated numerically by centered derivatives to improve computational time.
The optimization of the parameters was coded in a program with functions and subroutines, which start with a pair of initial values, to subsequently calculate the Jacobian matrix and the error vector. This allows to solve Equation (3), which results in a new pair of values. This allows to calculate the error vector and check if it is smaller than the previous iteration due to a tolerance. If it is higher than the tolerance, the algorithm starts again with the pair found in the last iteration, and if this error vector is lower than the tolerance, the program stops and prints the results. The diagram is shown in Figure 1.
Input parameters
The input data required are the initial and saturation moisture contents, as well as data derived from an infiltration test (time and infiltrated water depth). To determine the initial moisture content (θo), the input variables are bulk density of soil (ρt), gravimetric moisture content (ωo) and density of water (ρw):
Porosity (φ) is the void space volume of the porous media, which is calculated from bulk density (ρt) and quartz particle density (ρs):
where ρs = 2.65 g∙cm-3.
Saturation moisture content (θs) is the volume of water retained in the pore space, usually assimilated to the volumetric porosity (φ) by the following inequality: 0 ≤ θs ≤ φ; nevertheless, it is important to note that, for a soil apparently saturated with water, a certain amount of air is usually trapped. Therefore, the saturated water content can be taken θs = 0.9φ (Haverkamp et al., 2016; Rogowski, 1971). In this study, θs will be taken as the total porosity: θs ≈ φ.
Laboratory experiment
Soil samples were taken from agricultural plots in the region of San Juan del Río, Querétaro, Mexico, following the methodology proposed by Reynolds and Topp (2007). The soil was sieved with a number 10 mesh (2 mm) and was outdoor dried for one week. Before placing the soil in the column, a sample was taken and sent to the laboratory to determine the initial moisture content using the gravimetric method.
The infiltration test was carried out using two acrylic columns of different lengths and circular cross section. Both columns had a porous plate covered by a filter and placed at the bottom to retain the soil and allow water and air to escape (Figure 2). The inside of the columns were partially coated with wax to create roughness between the soil and the acrylic column. The soil was placed in the column in 5 cm thick layers at a density similar to that obtained in the field. A constant water depth was kept during the infiltration test.
Table 1, adapted by Saucedo et al. (2013), shows the mean values of some soil parameters (θo, θs, Ks and hf) based on soil texture, using the relationships proposed by Rawls et al. (1991). To determine these values, the soils used were analyzed at the laboratory with the Bouyoucos hydrometer using the methods of the Mexican standard NOM-021-SEMARNAT-2000 (Secretaría del Medio Ambiente y Recursos Naturales [SEMARNAT], 2002) and the texture triangle of the United States Department of Agriculture (USDA). Within the operation of the optimization algorithm, the following condition was included: if the difference of parameters in the current iteration with the previous one is lower than 0.00001 %, or the number of iterations is equal to 20, end the algorithm sequence and display the results.
Table 1. Mean parameters of the Green and Ampt infiltration equation according to soil texture.
| Textura del suelo / Textura del suelo |
Parameters / Parámetros | |||
|---|---|---|---|---|
| θo (cm3∙cm-3) | θs (cm3∙cm-3) | hf (cm) | Ks (cm∙h-1) | |
| Clay / Arcilla | 0.36 | 0.49 | 140.26 | 0.05 |
| Silty clay / Arcillo-limoso |
0.32 | 0.48 | 100.16 | 0.05 |
| Silty clay loam / Franco-arcillo-limoso |
0.26 | 0.49 | 60.12 | 0.15 |
| Clay loam / Franco-arcilloso |
0.25 | 0.48 | 36.00 | 0.4 |
| Sandy clay / Arcilla-arenosa |
0.25 | 0.42 | 25.72 | 0.5 |
| Silt / Limo |
0.14 | 0.50 | 30.52 | 0.8 |
| Loam / Franco |
0.20 | 0.46 | 20.04 | 1.5 |
| Silt loam / Franco-limoso |
0.17 | 0.55 | 30.07 | 1.0 |
| Sandy clay loam / Franco-arcillo-arenoso |
0.18 | 0.42 | 35.61 | 2.0 |
| Sandy loam / Franco-arenoso |
0.16 | 0.46 | 10.00 | 2.9 |
Results and discussion
Decoding the algorithm
The optimization algorithm was programmed using the Visual Basic.net language in the Microsoft Visual Studio 2019 integrated development environment. This program can be installed on any Windows 10 platform, and the speed in RAM required is minimal, so any computer with a mid-range processor and 4 GB of RAM supports it. The program consists of four main sections: a) main menu, b) displayed data from an infiltration test, c) graphical results of measured vs. optimized irrigation water depth and d) optimization results (Figure 3).
Table 2 shows the routines used for the proper operation of the computer program we developed. The task performed by each routine is also described.
Table 2. Routines used in the computer program.
| Titel / Título | Routine / Rutina | Task / Acción |
|---|---|---|
| File / Archivo | Import Excel data / Importar datos de Excel |
It opens data from an infiltration test in *.xls format (time and cumulative infiltration). / Abre datos de una prueba de infiltración en formato* .xls (tiempo y lámina infiltrada acumulada). |
| Open examples / Abrir ejemplos |
Data from an infiltration test is loaded, including input parameters. / Se cargan datos de una prueba de infiltración, donde también se incluyen sus parámetros de entrada. |
|
| Save optimized data / Guardar datos optimizados |
It saves the simulated data with the optimized parameters in *.xls format (time and cumulative infiltration). / Guarda los datos simulados con los parámetros optimizados en formato *.xls (tiempo y lámina infiltrada acumulada). |
|
| Exit program / Salir del programa |
Closes the program. / Cierra el programa. | |
| Soil data / Datos del suelo |
Input data / Datos de entrada |
It displays a window in which the user must enter the required input parameters (bulk density, initial moisture, water depth, infiltration test column length and soil texture). / Muestra una ventana en la cual el usuario debe introducir los parámetros de entrada necesarios (densidad aparente, humedad inicial, carga superficial, longitud de columna de la prueba de infiltración y textura del suelo). |
| Run / Ejecutar |
Run / Ejecutar |
It starts the parameter optimization process, while charting the iterations in real time. / Inicia el proceso de optimización de parámetros, al mismo tiempo que grafica en tiempo real las iteraciones. |
| Calculation memory / Memoria de cálculo |
Calculation memory / Memoria de cálculo |
It displays a window with the results for each iteration performed by the Levenberg-Marquardt algorithm. / Muestra una ventana con los resultados en cada iteración realizada por el algoritmo de Levenberg-Marquardt. |
| Help / Ayuda | Help / Ayuda | It opens a window displaying the authors' contact information. / Abre una ventana donde se muestran los datos de contacto de los autores. |
Infiltration test
The initial conditions of the experiment are described in Table 3, showing the results of the analyses conducted at the laboratory, as well as the lengths of the columns (L) filled with soil, diameters (D) and surface water depths (h).
Table 3. Physical characteristics of the study sites.
| Soil sample / Muestra de suelo |
D (cm) | ρt (g∙cm-3) | θs (cm3∙cm-3) | θo (cm3∙cm-3) | L (cm) | h (cm) | Texture / Textura |
|---|---|---|---|---|---|---|---|
| S1 | 9.2 | 1.0810 | 0.5921 | 0.1462 | 70.0000 | 6.0000 | Sandy loam / Franco-arenosa |
| S2 | 8.8 | 1.1588 | 0.5627 | 0.1259 | 85.0000 | 6.0000 | Clay / Arcillosa |
| S3 | 8.8 | 1.1713 | 0.5580 | 0.0280 | 85.0000 | 5.0000 | Sandy loam / Franco-arenosa |
With the data from the initial conditions of the experiments and the infiltrated water depth over a period of time, the parameters Ks and hf of the Green and Ampt equation were optimized according to the soil texture. The results of the optimization for each site are displayed in Figure 4, showing the fit of the Green and Ampt model to the measured values.
Figure 4. Comparison between infiltrated water depth measured and calculated with the Green and Ampt equation.
Table 4 shows the optimization achieved, and the root mean square error (RMSE), resulting from the infiltrated water depth measured and the theoretical infiltrated water depth. Errors lower than 0.45 cm are observed, showing that the optimization of parameters with initial values related to soil texture works adequately. It is important to highlight that, from the previous graph, it is possible to observe the good fit of the optimized model with the data measured at the laboratory.
Numerical validation
The validation of the numerical algorithm was performed to remove programming errors and check the consistency of correct solutions. Data reported in the literature were used, and the results were obtained using the one-dimensional solution of the Richards equation (Chávez et al., 2016; Fuentes et al., 2020). The data used correspond to a sandy loam soil from Montecillo, Mexico (Zataráin et al., 2003). Parameter values for this soil are: ρt = 1.3607 g∙cm-3,θo = 0.1391 cm3∙cm-3, h = 1.50 cm, φ = 0.4865 cm3∙cm-3, θs = φ and L= 70 cm. The parameters Ks and hf of the Green and Ampt equation were optimized using the computational program according to soil texture (Figure 5).
Figure 5. Comparison of the infiltrated water depth derived from the Richards equation and the one modeled withthe Green and Ampt equation.
When comparing the result and the one reported in the literature (Table 5), it is observed that there is no considerable difference, since RMSE (0.1953 cm) is considered the same, and this indicates that the result measured using the Green and Ampt model optimizing Ks and hf varies very little.
The parameter hf is characteristic of the Green and Ampt equation; therefore, it cannot be compared with any other parameter of other equations. However, the value is adapted depending on the proportions of sand, silt and clay in soil, so it can be considered inversely proportional to Ks, because when Ks decreases, the value of hf increases.
If only soil texture is considered for the surface irrigation design, the average Ks value for the sandy loam soil analyzed would be 2.9 cm·h-1 (Table 1). However, this type of soil is dominated by sand characteristics, which leads to Ks values higher than 1.5 cm·h-1. The range of variation between both values is very wide, which causes a substantial change in the characteristic values of infiltration.
Therefore, it is recommended to perform infiltration tests in soil columns to adequately represent the flow in a one-dimensional manner, since other tests, such as the double ring infiltrometer shows a two-dimensional flow and, in addition, there is a preferential flow along the cylinder walls, which in most cases is attributed to the infiltration velocity. The advantage of using columns with soil samples altered at the laboratory allows control of all variables, such as bulk density, soil depth, water depth, initial moisture, among others.
Temporal sensitivity of optimized infiltration test parameters
A sensitivity analysis of the Green and Ampt parameters was performed to verify how the estimation errors of the two parameters propagate to the infiltration estimates. Optimizations were made with the experimental results at different times for each site. For S1 a cumulative increment of Δt = 1 h, for S2 a Δt = 4 h and for S3 a Δt = 1 h. Figure 6 shows that the parameters have the same trend, while Ks gradually decreases, the value of hf increases. The inverse process occurs with the same trend; this indicates that the algorithm optimizes the values according to the soil texture and the more infiltration data one has, the optimized parameters can vary considerably, either decreasing or increasing their values gradually and inversely proportional to each other. It also shows the importance of performing an infiltration test at long times, because, if it is performed at short times, the pores are not filled with water, and one would only have a value of hydraulic conductivity and not the saturated hydraulic conductivity.
The program is delimited so that each pair of new values proposed are within the limits of the parameters, i.e., Ks )( 0 and 0 () hf ≤ 200. Parameters are adapted according to soil texture, with emphasis on the Equation (7):
where S is the sorptivity of the porous medium (Philip, 1957). Without the use of initial parameters related to soil texture, there can be an infinite number of parameters to fit the data measured in an infiltration test.
In this study, a computer with the following characteristics was used: Intel® CoreTM i7-4710MQ CPU @ 2.50 GHz and 32 GB RAM memory. The computation time required to find the optimal values of Ks and hf in the experiments was 2 s. This is because the Green and Ampt solution requires no knowledge of values at the previous time level; therefore, it is possible to solve it exclusively in the required time, as in the exact times where it is compared with the measured data. On the other hand, the Richards equation needs the pressure values at the previous time level, which delays, in at least one iteration, the calculation of the infiltrated water depth.
Conclusions
The estimation of hydrodynamic parameters of saturated hydraulic conductivity and wetting front suction was performed with the Levenberg-Marquardt algorithm. The main advantage of the presented optimization is the short computational time required to optimize the parameters of the Green and Ampt equation, compared to the Richards equation reported in the literature. It was shown that, by using infiltration test data from the literature, the values derived for Ks with the proposed nonlinear optimization algorithm are very similar to those found with the Richards equation, which provides certainty and reliability when calculating the parameters.
The sensitivity analysis suggests that the duration time of the infiltration test significantly affects the Ks and hf values, so it is recommended to perform an infiltration test up to the maximum infiltration, since variations in the infiltrated water depth with time may cause the curve to be adapted to a different pair of parameters..
