Using root water uptake estimated by a hydrological model to evaluate the least limiting water range

The least limiting water range (LLWR) has been extensively determined, but evaluating if LLWR can indeed indicate soil physical stress on plant growth is still a controversial issue. In this study, we used the Hydrus-1D hydrological model to simulate root water uptake (RWU) to analyze if RWU and LLWR are correlated under stress conditions. The LLWR was determined in a sandy-loam Ultisol and a clayey Oxisol. In both soils, RWU extracted by plants (leaf area index set as 3) from a rooted layer of 0.4 m was simulated over 20 days under a potential evapotranspiration rate of 6 mm day. For each soil, RWU was simulated over the same range of soil compaction in which LLWR was determined. The cumulative RWU over the 20 days varied between 23 to 58 mm in the Ultisol and 20 to 48 mm in the Oxisol, indicating that plants were able to take up only a small part of the cumulative potential transpiration (93 mm) and experienced severe water stress in some soil conditions. However, RWU under water stress was poorly correlated with both bulk density and LLWR. The correlation between RWU and LLWR was 0.5 (p<0.01) for the Ultisol and 0.22 (p<0.19) for the Oxisol, suggesting that LLRW has little (for Ultisol) or almost no (for Oxisol) ability to indicate soil quality related to plant water availability. Our simulations suggest that RWU in the water availability range (between field capacity and wilting point) may be little affected or even improved by light soil compaction. Studies to elucidate this phenomenon would contribute to the understanding of the compaction effect on RWU and the weak correlation between RWU and LLWR.


INTRODUCTION
An extensive review on least limiting water range (LLWR) by Gubiani et al. (2013a) showed that this soil physical quality indicator for plants was evaluated using specific plant parameters in only 9 of 38 publications. Considering studies in which plant parameters and LLWR were correlated, Gubiani et al. (2013b) stated that LLWR is generally weakly correlated to a single measurement of plant variables (e.g., leaf area, yield, and root length). Taking the plant process into account, a detailed theoretical analysis of LLWR by van Lier and Gubiani (2015) showed that LLWR is unable to mimic any mechanism of soil physical stress that plants may experience over time. However, it is possible to expect LLWR to be statistically correlated with some biological processes, and the correlation degree can be used to evaluate the usefulness of LLWR.
In this study, we used a physically-based model to simulate root water uptake (RWU) under stress conditions to analyze its correlation with LLWR. One of our hypotheses is that water flow from soil to roots may be improved by soil compaction. There is evidence showing that soil compaction can increase unsaturated hydraulic conductivity (K) within the range of water available to plants (Aravena et al., 2011;Durigon and van Lier, 2011). This is because soil compaction increases the volume of micropores that remain filled with water within the range of water available to plants (Richard et al., 2001). As K is defined by the water flow density that can be transported by a unit of gradient of energy, an increase in K caused by soil compaction might also increase RWU. As a result of the first hypothesis, the second one is that LLWR and RWU are poorly correlated.
Plant water uptake is a complex process influenced by soil, plant, and climatic factors (van Lier et al., 2013). Measuring water uptake in growing plants experiments is the best way to study how soil physical factors control this process. However, modeling plant water uptake using mechanistic models is also a useful exploratory strategy, avoiding laborious experiments. However, a few mathematical models describe water uptake by mimicking the physical and physiological mechanisms governing the water flow throughout the soil-plant-atmosphere continuum. One example is the Hydrus-1D model (Šimůnek and van Genuchten, 2008). Hydrus-1D numerically solves the Richards equation for variably-saturated water flow and also incorporates a sink term that accounts for RWU. Hydrus-1D mimics the process of plant water uptake by relating the potential transpiration, which is a function of atmospheric and plant parameters, with soil restriction to water flow, the so-called root-water uptake water stress response function. For example, Hydrus-1D has been used for evaluating the effect of salinity (Shouse et al., 2011) and the dynamics of root distribution (Cai et al., 2017)

on RWU.
It is possible to estimate the effects of soil compaction on water flow using Hydrus-1D because soil compaction influences soil water-retention curve parameters and saturated hydraulic conductivity, which are input parameters in Hydrus-1D. Therefore, the objective of this study was to evaluate our hypotheses by using Hydrus-1D to simulate RWU in two soils in which LLWR was also determined.

Soil sampling and measurements
Undisturbed soil samples were collected in a sandy-loam Ultisol (Argissolo) and a clayey Oxisol (Latossolo), both cultivated under no-tillage (the first and the second name are according to Soil Taxonomy (Soil Survey Staff, 2014) and Brazilian soil classification system (Santos et al., 2011), respectively). Soil texture in the 0.00-0.20 m layer is sandy loam (580 g kg -1 of sand, 300 g kg -1 of silt, and 120 g kg -1 of clay) in the Ultisol and very clayey (130 g kg -1 sand, 230 g kg -1 of silt, and 640 g kg -1 of clay) in the Oxisol. The soil samples were collected with metal rings (0.05 m long and 0.06 m in diameter) at a layer of 0.05-0.10 m. The collection points were randomly distributed over the sites to compose a set of samples in each soil with variation in bulk density. To do this, we collected a total of 108 and 93 samples in the Ultisol and Oxisol, respectively.
Excess soil from the rings was removed in the laboratory. A permeable screen was attached to one end of the rings and the samples were kept on trays with the water level near the upper edge of the rings for 48 hours for saturation. Afterward, we determined water retention (θ, m 3 m -3 ) and soil resistance to penetration (R, MPa) for several values of water tension (h, m). Bulk density (ρ, Mg m -3 ) and saturated hydraulic conductivity (K s , m day -1 ) were also determined.
Water retention was determined at h of 0.1, 0.6, and 1.0 m using a sand column (Reinert and Reichert, 2006) and at 3.3 and 10 m using a pressure-plate extractor (Klute, 1986). In this range of h, 32 samples of the Ultisol and 28 samples of the Oxisol were used to measure R. In addition, other samples drier than h =10 m were used to determine R. To do this, 18 samples of the Ultisol and Oxisol were removed from the pressureplate extractor after equilibrating to h of 10 m, and more water loss was allowed by evaporation. After the water loss and before measuring R, the samples were kept in plastic packaging for three days for redistribution of the water content. Soil resistance to penetration was measured using an electronic penetrometer with a metal rod (cone of 30° angle and 4 mm diameter at the base), inserted in the soil sample at a constant rate of 10 mm min -1 . Soil resistance to penetration was measured at two or three values of h in some samples. Each penetration creates a hole of 4 mm of diameter. The area of three penetration is just 1 % of the area of a sample with 60 mm of diameter. Thus, we assumed that penetrations did not affect each other. With this strategy, the number of measurements went beyond the number of samples, totaling 92 R measurements in the Oxisol and 67 in the Ultisol. All samples were then dried at 105 °C until constant weight.
The K s was determined with a constant head permeameter (Libardi, 2005) in 40 samples of the Ultisol and in 37 samples of the Oxisol, in which R was not measured. These samples were then dried at 105 °C until a constant weight was reached. In all samples, we calculated ρ by dividing the dry weight of soil by the volume of the ring.
The water retention curve was determined only in the 40 Ultisol samples and in 37 Oxisol samples in which K s was determined. A value of θ for h of 150 m was included in the set of θ-h obtained in these samples in a sand column and pressure-plate extractor. This value for each soil was obtained from measurements with a dew point psychrometer by Gubiani et al. (2012), as described below by equations 12, 13, and 14. Equation 2 described in the next section was fitted to the total set of θ-h of each soil.

Water uptake by plant roots simulated by HYDRUS-1D
Hydrus-1D numerically solves the Richards equation for variably-saturated water flow, and incorporates a sink term to account for water uptake by plant roots (Šimůnek et al., 2013): in which h is water pressure head (m), θ is volumetric water content (m 3 m -3 ), t is time (day), x is the vertical coordinate (m) (positive upward), S is the sink term (m 3 m -3 day -1 ), and K is the unsaturated hydraulic conductivity function (m day -1 ). Note that h is defined as a positive value and will hereon be referred to as water tension.
The θ (h) and K(h) relations in equation 1 were defined according to the analytical functions of van Genuchten-Mualem: Rev Bras Cienc Solo 2020;44:e0190096 in which S e is effective saturation (dimensionless), θ s and θ r are saturated and residual water contents, respectively (m 3 m -3 ), K s is saturated hydraulic conductivity (m day -1 ), α (m -1 ), and n are semi-empirical shape parameters, m = 1-1/n. The values of θ r , α, and n were determined by fitting equation 2 to each sample for both soils. This allows considering the effect of texture and ρ on parameters and on the simulations of RWU. The value of γ, a pore-connectivity parameter (dimensionless), was fixed in HYDRUS-1D simulations at 0.5. The pore-connectivity is affected by texture and aggregation (Pinheiro et al., 2019), and the value of γ should be different in Ultisol and Oxisol, and over the range of ρ for each soil as well. However, the θ-h-K relation needs to be determined in each sample, under an evaporation experiment for example, to calibrate γ. Due to the lack of appropriate equipment, it was not possible to calibrate γ in this study. Thus, the differences in θ (h) parameters were responsible for transmitting the effect of texture and ρ on RWU simulations.
The sink term S, in equation 1, is defined as the volume of water removed per unit of time from a unit volume of soil as a result of plant water uptake. The simulations of this study were carried out using the Feddes et al. (1978) approach available in HYDRUS-1D, according to which: in which the root-water uptake water stress response function ω(h) is a prescribed dimensionless function of the soil water tension, 0 < ω ≤ 1, and S p the potential water uptake rate (day -1 ).
The root-water uptake water stress response function, ω(h), is assumed to be zero close to saturation (h < h 1 ) and after water tension close to wilting point (h > h 4 ); ω(h) is optimal (=1) between water tensions h 2 and h 3 , whereas for water tensions between h 3 and h 4 , and between h 1 and h 2 , ω(h) decreases (or increases) linearly with h. In this study, we set h 1 = 0.1 m, h 2 = 0.25 m, h 3 = 2 m, and h 4 = 80 m, which are values close to those suggested by HYDRUS-1D for corn in the vegetative period.
The potential water uptake rate S p was distributed over the root zone using a non-uniform normalized distribution function: in which T p is the potential transpiration rate (m day -1 ) and b(x) is the normalized water uptake distribution in depth (m -1 ) available in Hydrus-1D: in which z is the x-coordinate of the soil surface [L] and z r is the root depth [L].
The application of equation 6 considers that the bottom of the soil profile is located at x = 0 and the soil surface at x = z. In this study, z r was set at 0.4 m and maintained Rev Bras Cienc Solo 2020;44:e0190096 constant over the simulation period of 20 days. Using this configuration, b(x) was 0.042 m -1 in the range of 0 ≤ z ≤ 0.08 m, from which it decreased linearly to zero at z = 0.4 m.
The potential transpiration rate is estimated by HYDRUS-1D using the following equation: in which ET p is potential evapotranspiration (m day -1 ), LAI is the leaf area index (dimensionless), and k is a constant governing the radiation extinction by the canopy (dimensionless), varying between 0.5 and 0.75. In this study, k was set at 0.5, LAI was set at 3, and ET p was set at 0.006 m day -1 , all maintained constant over the simulation period of 20 days.
Total root water uptake (RWU, m) at the end of 20 days was calculated by substituting equation 5 in equation 4, and integrating it over depth and time as follows: HYDRUS-1D simulations were run using a spatial discretization involving 101 nodes with the same density over the 0.4 m profile depth. The upper boundary was atmospheric boundary conditions (no precipitation, default value of 10 3 m for h crit related to evaporation, LAI of 3, and ET of 0.006 m day -1 ), while the bottom boundary was free drainage. For the initial condition, the water content was considered to be closer to field capacity by setting an equilibrium profile with a water tension of 1 m. It is important to mention that no plant growth and constant atmospheric conditions were assumed during simulations because our hypotheses are related to the variation of soil compaction degree.

Determining the least limiting water range
LLWR is defined by (van Lier and Gubiani, 2015): in which "max" is the "maximum" function [max(i,j) = i if i ≥ j; max(i,j) = j if i < j] and "min" is the "minimum" function [min(i,j) = i if i ≤ j; min(i,j) = j if i > j]. The θ cc is the water content at field capacity, θ pa is the limiting water content for adequate aeration, θ pwp is the water content at the permanent wilting point, and θ pr is the water content corresponding to the limit of the resistance to penetration.
The θ pa was considered θ corresponding to 10 % air-filled porosity: Particle density (ρ p , Mg m -3 ) was determined in previous studies in both soils (2.65 Mg m -3 and 2.7 Mg m -3 for the Ultisol and Oxisol, respectively).
The θ cc was estimated by relating θ at h = 1 m (θ 100 ) to ρ by means of a linear equation: in which a 1 and b 1 are fitting parameters.
The θ pwp was considered as the θ at h = 150 m, which was related to ρ as follows: in which W 150 is gravimetric soil water content (Mg Mg -1 ) at h = 150 m and ρ w is water density (1 Mg m -3 ).
W 150 was calculated using h = 150 m in equations of Gubiani et al. (2012), which were generated using a large set of W measurements in a 50 < h < 2500 range for both soils used in this study: For the Ultisol: W 150 = 0.3148h -0.3753 (R² = 0.95) Eq. 13 For the Oxisol: W 150 = 0.347h -0.1309 (R² = 0.85) Eq. 14 The θ pr was considered the θ corresponding to an R of 2 MPa. The nonlinear model proposed by Busscher (1990) was fit to the dataset of R, ρ and θ of each soil: in which a 2 , b 2 , and c 2 are fitting parameters; a 2 in MPa m 3n kg -n , and b 2 and c 2 are dimensionless.
Rewriting (15), θ pr for R of 2 MPa was calculated as a function of ρ: The LLWR was calculated by equation 1 using θ cc , θ pa , θ pwp , and θ pr of each soil sample.
The fitting of equations 2, 11, and 15 was done by minimizing the sum of squared residuals. The RWU was associated with ρ by means of regression and with LLWR by the Pearson correlation coefficient. The K estimated with equation 3 at a tension of 1, 10, and 100 m was correlated with ρ by means of the Pearson correlation coefficient. All statistical calculations were made using SAS software (SAS, 1999).

RESULTS
Among samples, ρ varied from 1.46 to 1.88 Mg m -3 for the Ultisol and from 1.27 to 1.62 Mg m -3 for the Oxisol. In these ranges of ρ, θ cc , θ pa , θ pwp , and θ pr (Figures 1a and 1b) were consistently determined by equations 10, 11, 12, and 16, respectively. The effect of ρ was significant on θ cc for the Oxisol and on R for both soils (Table 1). As the effect of ρ was not significant on θ cc for the Ultisol (Table 1), the average value of θ cc was used. The θ fc , θ pwp , and θ pr increased with higher ρ, while θ pa decreased as ρ increases (Figures 1a and 1b). The parameters θ pr and θ pa were the most sensitive to the increase of ρ and caused a remarked decrease on LLWR with increasing ρ, which was zero at a ρ of 1.71 Mg m -3 for the Ultisol (Figure 1c) and 1.42 Mg m -3 for the Oxisol (Figure 1d).
The cumulative RWU over the 20 days varied between 23 to 58 mm in the Ultisol and 20 to 48 mm in the Oxisol. These values are between 21 and 62 % of the cumulative Table 1. Fitted equations for water content at field capacity (θ fc ) as a function of bulk density (ρ) (Equation 11) and the equation for soil resistance to penetration (R) as a function of ρ and soil volumetric moisture (θ) (Equation 15)

Soil
Equation 11 Equation 15 Ultisol θ fc = -0.1347ρ + 0.  Rev Bras Cienc Solo 2020;44:e0190096 potential transpiration (93 mm), and they show that plants would experience severe water stress. The stress-factor for root-water uptake (ω) indicated that RWU ceased at 0.1 m depth, and it was lower than potential RWU at 0.3 m depth (Figure 3c). There was little association between the simulated RWU and ρ (Figures 1c and 1d). Although the regression fitting coefficients were significant (p<0.05), R 2 shows that only 43 and 14 % of the variance was explained by the models for the Ultisol and Oxisol, respectively.
The correlation between RWU and LLWR by Pearson's analysis was 0.5 (p<0.01) for the Ultisol and 0.22 (p<0.19) for the Oxisol, suggesting little association for Ultisol and almost none for Oxisol. Also, the positive correlations between Log 10 K and ρ in both soils were at most 0.53 (Figure 2), but it was stronger in the Oxisol.

DISCUSSION
The relationship between the LLWR parameters and ρ (Figures 1a and 1b) is in agreement with observations in the literature (Silva et al., 1994;Tormena et al., 1998;Junior et al., 2012;Kaiser et al., 2009). The values of ρ for which LLWR was zero (1.71 Mg m -3 for the Ultisol and 1.42 Mg m -3 for the Oxisol) are similar to those found previously in these same soils, 1.75 Mg m -3 for the Ultisol (Kaiser, 2010) and 1.48 Mg m -3 for the Oxisol (Gubiani et al., 2013b), indicating LLWR was correctly determined.
The weak statistical association between RWU and ρ suggests a decrease in RWU with increasing ρ, except before the maximum value of the parabolic function relating RWU with ρ for the Ultisol (Figure 1c). Except in this range for the Ultisol, the RWU simulated with Hydrus-1D refutes our hypothesis that water flow from soil to roots may be improved by soil compaction. In this study, the differences in θ(h) over the range of ρ were responsible for transmitting the effect of ρ on RWU simulations. Some possible effect of differences in the pore-connectivity parameter γ on RWU over the range of ρ was not assessed in our study since γ was set as 0.5. Moreover, the same imposed initial conditions by setting an equilibrium profile with a water tension of 1 m may not be a proper field capacity for both soils. Further investigation would be useful for reevaluating this hypothesis by calibrating the pore-connectivity parameter γ and soil field capacity.
The weak correlation between RWU and LLWR (of 0.5 for the Ultisol and 0.22 for the Oxisol) corroborates our hypothesis that LLWR and RWU are poorly correlated. The weak association between RWU and LLWR (Figures 1c and 1d) means that LLWR is unable to indicate if the soil is impairing the flow of water towards the root. Thus, the soil quality regarding its ability to move water towards the roots cannot be assessed by the LLWR.
Due to the poor correlation between RWU and LLWR, the possibility of LLWR being associated with crop performance is highly dependent on its ability to indicate the mechanical constraint that roots would face during growth, which is estimated by equation (16). However, equation (16) is highly uncertain, since what it expresses is the soil resistance to the penetration of a metal cone and not the resistance that the roots face during growth (van Lier and Gubiani, 2015). Plant physiological and morphological mechanisms (e.g., mucilage production and change growth direction) and soil factors (e.g., presence of biopores and cracks) allow the roots to perceive the mechanical resistance less than that detected by the metal cone (Bengough, 1997;Clark et al., 2003;Bengough et al., 2006Bengough et al., , 2011. Furthermore, LLWR also considers an abrupt biological change in θ pr and does not measure mechanical stress. This limitation has been tackled by using a continuous mechanical stress function which relates the root elongation rate to soil resistance to penetration of a metal cone, also introducing the bioporosity effect (Moraes et al., 2018). However, the LLWR abrupt approach does not deal with continuous stress functions. Moreover, for only assessing mechanical constraints, an easier soil resistance to penetration function determination has the same capability of the laborious LLWR.
Experimental findings support the positive correlations between Log10K and ρ in both soils. Durigon et al. (2011) observed that K measured with polymer tensiometers in an evaporation experiment remained higher in the sample with higher ρ, except at the wet end of the water retention curve. From a certain tension, a higher ρ would allow a larger volume of pores filled with water and effectively contributing to the water flow (Richard et al., 2001). It is a fact that aggregates and particles approach each other with increasing ρ. This improves water connectivity and its flow in situations of low soil water content (Carminati et al., 2008;Aravena et al., 2011). Both phenomena (conversion of large to small pores and increased contact between particles) are less intense when a large part of the soil that moves with the compaction is composed of large, non-porous particles such as sand. This would explain the lower change in K in the Ultisol (580 g kg -1 of sand) in comparison to that of the Oxisol (130 g kg -1 of sand).
The increase of K with increasing ρ shows that a unit of gradient of energy can transport more water if there is an increase in compaction. Although the direct relationship of K with ρ did not reflect in a direct relationship of RWU and ρ, it must have reduced the curvature of the parabolic (Oxisol) and the slope of the linear (Ultisol) functions relating RWU with ρ ( Figure 1). However, as the positive effect of ρ on water transport is unclear, its negative effect can not be supported by simulations output as well. At depths of 0.1 and 0.3 m at 20 th day of simulation, volumetric water content, water tension and stress-factor for root-water uptake were not correlated to ρ for both soils (Figure 3). This means that predicted soil hydraulic conditions were independent of the soil bulk density, that is, RWU was not sensitive to soil compaction.
If an increase in ρ leads to some compensation in RWU this can temporarily relieve some restriction in water uptake due to mechanical impediment to root penetration. Depending on how often this compensation occurs, it can relieve a lot of stress during the plant growth cycle. The plant is less forced to reduce biomass accumulation rate or to eliminate leaf area, both undesirable agronomic consequences of the physiological mechanisms acting to reduce water stress (Taiz and Zeiger, 2016). Our simulations suggest that water transport in the water available range (between field capacity and wilting point) could be little affected or even be improved by light soil compaction. Studies to elucidate this phenomenon would contribute to the understanding of the compaction effect on RWU.

CONCLUSION
The Least Limiting Water Range (LLWR) has little indicative ability for the soil water suppling capacity to roots, and crop water stress to root water uptake cannot be assessed by the LLWR.
Correlations between LLWR and plant performance are highly dependent on the ability of LLWR to indicate the mechanical constraint that roots would face during growth. However, a simple to perform soil resistance to penetration function determination has the same capability of the laborious LLWR for assessing mechanical constraints.
Studies to investigate how compaction affects water transport in soils between field capacity and wilting point would contribute significantly to understanding the compaction effect on RWU and the limitations of LLWR.