Parameter estimation of soil hydraulic characteristics by inverse modeling of the analytical equation for unsaturated subsurface water flow

Localized irrigation has increased in popularity in recent decades due to the low emitter flow rates, high water application efficiency, and possibility for fertirrigation and systems automation (Li et al. 2004; Arbat et al. 2013; Moncef & Khemaies 2016; Arraes et al. 2019). This irrigation method has been used for crops with high added value and in locations with water scarcity, and regions in conflict due to the use of water resources.

Success in the design of this irrigation modality depends on the design of system stages and on water application management (Samadianfard et al. 2012; Subbaiah 2013). Due to its high implementation cost, it is necessary to acquire information on the soil to be irrigated and on the hydraulic characteristics of the emitter to be used in an attempt to reduce failures in the design phase.

The soil physical and hydraulic characteristics, combined with the emitter flow rate and the frequency of irrigation, define the shape and size of the moist soil volume below the water application point. Thus, these characteristics are the most influential variables in the design and management of drip irrigation systems (Cook et al. 2006; Hao et al. 2007; Kandelous & Sǐmůnek 2010a, 2010b; Arbat et al. 2013; Tolentino Júnior et al. 2014).

The wet part of the wetting bulb where the root zone of plants is concentrated is called the wetting pattern (Elnesr & Alazba 2017). Soil wetting patterns under a single emitter can be used to suggest planning factors for a drip irrigation system.

The design of drip irrigation systems involves the selection of an appropriate combination of emitter flow rate and spacing between emitters. This combination allows the understanding of the wetting pattern (Subbaiah 2013; Abid 2015), that is the area from where plant roots absorb water and nutrients (Kunz et al. 2014). The spatial distribution of roots in a drip-irrigated crop is closely related to the dimensions of the wetting bulb (Souza & Matsura 2004; Li et al. 2017). This provides important information for defining the percentage of the area to be effectively irrigated (Li et al. 2004), spacing between emitters (Levien et al. 2012; Al-Ogaidi et al. 2016), and the depth of installation of the emitters in the case of subsurface irrigation (Kandelous & Sǐmůnek 2010a, 2010b; Santos et al. 2016; Ren et al. 2017; Vásquez et al. 2017).

The wetting pattern can be obtained by directly measuring the soil moistened in the field, which is site-specific, or by simulation using models (Souza et al. 2007; Abid et al. 2012). Field tests consist of excavating a trench in the soil, perpendicular to the lateral drip line, to measure the dimensions of the wetting bulb and to facilitate the installation of sensors for determining soil moisture (Barros et al. 2009; Bizari et al. 2016). However, field tests are expensive, and it is sometimes difficult to visually distinguish the wet soil from dry soil (Arbat et al. 2013).

The use of mathematical models is an important alternative for estimating the water distribution in the area delimited by the wetting bulb (Coelho et al. 1999; Abid et al. 2012; Abid 2015). The use of computer simulations with mathematical models validated to optimize management practices is a quick and inexpensive approach that allows the understanding of water distribution in the soil in drip-irrigated fields (Li et al. 2015).

Numerous empirical, analytical, and numerical models have been developed to simulate the wetting pattern and the dimensions of wetting fronts in drip-irrigated soil, both in permanent and transient flow regimes (Kandelous & Sǐmůnek 2010a, 2010b; Levien et al. 2012; Subbaiah 2013; Abid 2015).

Analytical models quickly allow determination of the position of the wetting pattern (Cook et al. 2006) and are based on the hypothesis of a point source of water and the use of analytical relationships to characterize some soil physical and hydraulic properties (Cook et al. 2003; Arbat et al. 2013). However, these models are limited by the number of variables involved and therefore require a smaller number of input parameters, and they are represented by more accessible equations than those used in numerical models (Coelho et al. 1999). Numerical models, in turn, require fewer assumptions (Cook et al. 2006), yield more accurate results, and are computationally efficient (Ku et al. 2017). Nonetheless, they require more input parameters in comparison to analytical models. For this reason, analytical models demonstrate widespread use (Moncef & Khemaies 2016; Elnesr & Alazba 2017).

Accurate knowledge of soil hydraulic properties is crucial to solve soil water flow problems under unsaturated conditions (Jhorar et al. 2002; Iden & Durner 2007). The hydraulic soil properties used in simulation models involve the relationships between hydraulic conductivity, soil matric potential and moisture (Communar & Friedman 2015; Ameli et al. 2016). The determination of soil hydraulic properties is difficult, requiring methods that are time-consuming and expensive. Meanwhile, the parameter estimation method requires the use of complex mathematical analysis methods and numerical techniques (Silva Junior 2015). An alternative is the use of inverse modeling methods, in which the modeling results are compared with data measured experimentally (Nakhaei & Šimůnek 2014; Vogelerb et al. 2019). Optimization techniques are used in the determination of soil hydraulic properties by the inverse method, which aims to minimize the sum of squares of the deviations between estimated and experimental data on the study variables (Hopmans et al. 2002; Olyphant 2003). In the studies performed by Jhorar et al. (2002), Ritter et al. (2003), Antonino et al. (2004), Moret-Fernández & Latorre (2017), Costa (2019), and Latorre & Moret-Fernández (2019), the authors obtained accurate estimates of soil hydraulic parameters by using the inverse modeling technique for determining the water movement in unsaturated soil. Thus, this alternative method allows the acquisition of soil hydraulics parameters quickly and at a reduced cost.

Given the above findings, the aim of this study was to determine the hydraulic parameters of a clayey Yellow Oxisol by the inverse modeling of an analytical equation used in wetting bulb simulation based on soil moisture measurements obtained in the field.

In this study, the analytical solution of Equation (4) proposed by Warrick & Lomen (1976) was used (Equation (5)). This equation is applied for point sources of water application, as is the case of drip irrigation, and subject to the following initial and boundary conditions (Figure 1):

• for T = 0 and 0 ≤ z ≤ ∞: Φ = 0 and h = –∞ (extremely dry soil);

• for T > 0 and R and Z → ∞: Φ = 0 and h = –∞ (points distant from the source and not reached by the wetting pattern);

• for T > 0 and Z = 0 and for any value of R, the water flow at the soil surface is null, i.e. q_z = 0, except at the source where there is water application.

  

  (5)

This study was performed at the experimental field of the Federal University of Goiás (Universidade Federal de Goias – UFG). The coffee cultivars Catimor, Catucaí, Iapar59, Obatã, and Oeiras were cultivated and drip-line irrigated at a flow rate of 3.0 L h^–1 and average water application time of 2 h. The soil of the experimental field was classified as clayey Yellow Oxisol (42.5% clay, 13.1% silt, and 44.4% sand), and its hydraulic characteristics were determined in the laboratory. For this purpose, undeformed soil samples were collected to determine the hydraulic conductivity of saturated soil and soil water retention curves.

The saturated hydraulic conductivity (K_o) was determined using a constant-head permeameter. For the retention curve points at pressures of 1, 2, 4, 6, 8, and 10 kPa, a suction unit with a Buchner funnel was used; at pressures of 33, 60, 100, 500, and 1,500 kPa, a Richards' pressure chamber was used. The K_o was 10.5 cm h^–1 and classified as moderately fast according to Ferreira (1999). The parameters of the van Genuchten model for prediction of the moisture retention curve, adjusted using the RETC application (van Genuchten et al. 1991), were θ_s = 0.443 m³ m^–3, θ_r = 0.162 m³ m^–3, α = 0.0497, n = 0.7188, and m = 3.556.

A trench was opened perpendicularly to the lateral line to identify the wetting bulb, which was visible after one irrigation cycle. After identification of the formed wetting bulb, the soil was sampled with the aid of a hand screw auger, in a 30 × 60 cm grid, sampling points 5 cm apart. For this purpose, the grid of sampling points was traced from the flow emission point (dripper), considering a symmetric distribution of the wetting bulb. Soil samples were collected during the dry season to facilitate visualization of the wetting bulb formed in the soil from the water application source.

The following options were considered for parameters estimates:

• estimating parameter β while considering the other parameters, K_o, θ_r, θ_s, α, n, and m, as known and obtained experimentally;

• estimating parameters K_o and β while considering the experimental retention curve as known;

• estimating parameters K_o, β, α, n, and m while considering the values of the volumetric moistures θ_r and θ_s as known;

• estimating all the parameters of the Gardner and van Genuchten models (K_o, β, θ_r, θ_s, α, n, m).

Table 1 shows the results of the parameters estimation of the Gardner and van Genuchten models by the inverse modeling of the analytical equation proposed by Warrick (1974) and Warrick & Lomen (1976) to evaluate the distribution of the soil matric flow potential from a point source of water application. These models describe the dependence of unsaturated hydraulic conductivity and volumetric moisture as a function of matric potential. The results show that the concordance coefficients (c) were greater than 0.9 for all the options employed to estimate the parameters. Thus, according to Camargo & Sentelhas (1997), the performance of the estimates are classified as optimum.

The best concordance coefficients (c) between the soil moistures observed in the field and estimated by the inverse modeling were obtained by option D, in which all the parameters of the Gardner and van Genuchten models were adjusted. Option C showed the second best performance, in which the θ_r and θ_s moisture values were fixed and the other parameters were adjusted. Option A, in which the experimentally known parameters (K_o, θ_r, θ_s, α, n, and m) were fixed, adjusting only parameter β of the Gardner model showed the lowest concordance coefficient among the evaluated options. Regarding option B, the concordance coefficient was close to that obtained for option A due to the proximity of the values of the adjusted β and K_o parameters of the Gardner model. This result demonstrates that the benefits of the decrease of the deviations between the moisture observed in the field and the estimated values were reduced when opting for fitting the moisture retention curve.

These variations can be explained in part by the methodological issues employed in obtaining the saturated hydraulic conductivity (constant-head permeameter) and the moisture retention curve (suction units with the Buchner funnel and Richards pressure chamber) in comparison to the results obtained by modeling the equations that describe the movement of water in the soil. The soil control volume evaluated in the modeling exceeded the size of the samples used in the characterization of hydraulic soil properties in the laboratory. Another point is that, in the interior of the soil, the water is subjected to the action of potential gradients between the points of the sampling grid at different energy potentials and moisture, whereas in the laboratory, the samples are subjected to the forcing action of the pre-established pressures to extract water from its pores. The different values of the matric potentials at the grid points were reflected in the unsaturated hydraulic conductivity values, which were mainly indicated in the value of parameter β.

In the determination of K_o in the laboratory, maintaining the water head constant at the top of the soil column results in an increase in the total potential that forces the movement of water through the soil column in the vertical direction. In contrast, in the field without the presence of a hydrostatic head, a thin zone of saturated soil is formed below the dripper emission point. From this saturated zone, the wetting bulb spread in the soil due to the action of the potential gradients. Because of that, the estimated values by options C and D were lower than the measured values in laboratory. Agreeing with this, Mesquita & Moraes (2004) indicated that K_o is highly influenced by sample size, flow geometry, sample collection process, and other soil attributes. Therefore, the fit of the parameters when using option D was more adequate, avoiding expenditures of time and resources to characterize the properties of the evaluated soil.

Figure 2 shows the dispersion of the volumetric moisture values around the 1:1 line for the grid points for the four options used to adjust the model parameters that characterize the soil hydraulic properties, confirming the concordance between the measured and calculated values (Al-Ogaidi et al. 2016; Elnesr & Alazba 2017). Options A and B showed a greater dispersion of the data relative to the 1:1 line for the peripheral points of the wetting bulb where the moistures were lower. When using options C and D, the dispersion from the line 1:1 were lower for points with lower moisture, thus better representing the distribution of moisture in the soil from the point source of water application.

Figure 3 shows the observed and estimated moisture profiles of the advancing wetting front by the analytical solution developed by Warrick & Lomen (1976) using the four options for estimating the parameters of the Gardner and van Genuchten models. The moisture isolines were obtained using the kriging technique, a nonbiased interpolator that provides the lowest variance for the estimated values (Yamamoto & Landim 2013). For all moisture profiles shown in Figure 3, the theoretical semivariogram was fitted to a linear model, not reaching the plateau, the point at which there was no dependence on the moisture values with the distance from the point of water emission (Figure 4). This behavior of the theoretical semivariograms allowed us to conclude that the soil moisture levels within the grid of points were dependent on the values of their neighboring points, promoting the gradients of potentials responsible for the movement of water in the soil. At points beyond the perimeter of the wetting bulb, the soil moisture did not depend on the application of water by the point source but rather, it relies on the initial moisture conditions that depends on other sources of input water to the soil such as rainfall prior to the soil sampling period, as well as on evaporation between the coffee planting rows and on internal soil drainage.

Figure 3 shows good concordance when comparing the experimental positions measured at the advancing wetting front with the estimated results. The soil wetting pattern was more prominent in the vertical direction. The soil was classified as of clayey texture, which would imply a more prominent wetting pattern in the horizontal direction because of the matric potential gradient. However, this soil type had a granular structure that endowed it with high porosity, conferring a saturated hydraulic conductivity near the emission point classified as moderately fast that resulted in vertical movement of water in the soil in response not only to the matric potential gradient but also to the gravitational potential. This behavior has also been observed by Cook et al. (2006), Abid (2015), and Elnesr & Alazba (2017). Irrigation time had a greater effect on vertical water advance in soil compared with radial advance, and according to Elnesr & Alazba (2017), increasing the irrigation time tends to increase the depth of the wetting pattern, which is not preferable in most cases to avoid water loss by percolation to below the limits of the roots.

Similarly to Arbat et al. (2013), we evaluated the absolute mean errors (AMEs) for all points of the sampling grid and for grid point intersecting vectors forming angles of 0, 30, 45, 60, and 90° relative to the soil surface to verify the horizontal and vertical advance of the wetting bulb. Table 2 shows the results of the AMEs, in which one can observe that the values were low, confirming that the inverse modeling technique using the analytical solution developed by Warrick & Lomen (1976) is a plausible alternative for fitting the parameters that describe soil hydraulic properties. Singh et al. (2006), Kandelous & Sǐmůnek (2010a, 2010b), Arbat et al. (2013), and Arraes et al. (2019) evaluated the performance of analytical and numerical models in the simulation of the wetting pattern of point sources of water application. These authors verified the suitability of the evaluated models, with AME values ranging from 0.0002 to 0.045. Therefore, the AME values obtained in this study for all options used to estimate the parameters of the Gardner and van Genuchten models were within the range found in the cited studies.

When evaluating the options used for fitting the parameters, the AME value determined using all grid points was lower when Option D was used, followed by options C, B, and A. In turn, when evaluating the advances of the wetting front in the horizontal direction (0°), vertical direction (90°), and for the 30° angle, the difference between the AME values was small, a behavior similar to that observed by Elnesr & Alazba (2017). For the 45 and 60° angles, the same trend in AME values was observed when evaluating the four fitting options employing all grid points. This suggests that C and D options were the most adequate in attenuating the differences between the values of the estimated moisture contents relative to those observed at the points inside the wetting bulb, which contains matric and gravitational potential activities. The variations between the observed and estimated moisture values were lower in the sampling grid points closer to the emission point and increased with distance, as observed in Figure 5 for the vertical and radial distances of 15 and 30 cm, confirming the previous discussion. Considering the spatial heterogeneity of the soil properties in the field, it can be concluded that the correspondence between the simulations and observations was very good (Kandelous & Sǐmůnek 2010a, 2010b).

In this work, the parameters of Gardner and van Genuchten models were estimated by the inverse modeling of the analytical equation for subsurface water flow in unsaturated soil. The inverse modeling used field observations of the soil moistures measured in a sampling grid around the wetting bulb formed in the soil below the emission point of a dripper. From the obtained results, it was found that estimating all the parameters of Gardner and van Genuchten models (K_o, β, θ_r, θ_s, α, n, and m) showed the higher value of the concordance index c (0.9437) and the lower value of the absolute mean error (0.008) between the soil moisture estimated by the analytical model and observed in the field. Estimating parameters of soil hydraulic characteristics by this technique requires only a few measurements of soil moisture in the field. Therefore, the inverse modeling of the analytical equation of wetting bulb simulation represents a gain in time and resources spent in laboratory experiments. Consequently, this is an important tool for irrigation design and management, since the understanding of the wetting pattern of the soil is essential in defining the area effectively irrigated, the spacing between emitters and their installation depth, and the irrigation rate.

All relevant data are included in the paper or its Supplementary Information.