Modelling soil water retention to identify management inﬂ uence on soil pore characteristics

: The water retention curve is widely used in studies involving soil. The management systems directly infl uence soil structure by altering water retention dynamics. Several equations are used to adjust the retention of water in the soil, but most of the time, the choice of these models occurs in an arbitrary way. From this problem, it was proposed with the present study to relate the best mathematical model to water retention, taking into account the different management systems adopted, based on previously established adjustment criteria. For the accomplishment the study, a soil of caulinitic mineralogy and average texture was utilized. The treatments were area under native Cerrado (Savanna), eucalyptus plantation with six and twelve years of cultivation, pasture with two and six years of cultivation, conventional plantation with two and eight years of cultivation, no-till with three and six years of cultivation. From the adjustment criteria for non-linear models, it was found that the management infl uences the choice of the water retention model. It is possible to observe that the greatest divergences between the models occurred close to the soil saturation zone, and that the Fredlung-Xing model is more effi cient in adjusting the water retention curve under conservationist management systems.


INTRODUCTION
Soil water retention curve (WRC) is of fundamental importance (Botula et al. 2014) since the water content infl uence many physical and hydraulic properties of soil (Dias Junior & Estanislau 1999, Alaoui et al. 2011. WRC is used to determine pore size distribution, fi eld capacity (FC), permanent wilting point (PWP), and available water (AW), among others.
Water retention can be considered an important indicator of the proper management of soil and water resources, as well as environmental quality. This attribute derives from diverse interactions between soil properties (Rieu & Sposito 1991, Dias Junior & Estanislau 1999, Liyanage & Leelamanie 2016, Silva et al. 2017) that are susceptible to management (Alaoui et al. 2011).
Management systems are closely related to the dynamics of water in the soil (Kutílek 2004, Alaoui et al. 2011, Silva et al. 2018, considering the direct infl uence on soil structure (Liu et al. 2005, Alaoui et al. 2011, Pires et al. 2012. In this sense, it is important to emphasize that soil water retention is governed by two types of forces: adsorption and capillary (Marshall et al. 1996). Capillarity is related to the WRC head and depends on the porous geometry of the medium (shape, size, orientation and pore distribution), which is directly affected by the management system adopted (Rasiah & Aylmore 1998). Adsorption, another phenomenon linked to water retention, is strongly related to the texture of the soil (Hillel 1998, Durner 1994, which is not significantly affected by different management systems in the same soil.
Several mathematical models for WRC have been established (Too et al. 2014). The most popular was developed empirically (van Genuchten 1980) and alternative models were explored from the principles of soil physics (Fredlund & Xing 1994, Omuto 2009). The van Genuchten's model is arbitrarily applied in very diverse circumstances in Brazil, without considering any adjustment or comparison with other WRC models (Lucas et al. 2011, Barros et al. 2013. There are several factors to consider when evaluating an adjustment method for WRC. Statistical methods can be used with generalization and simplicity to assess the descriptive adequacy of the model. Nevertheless, the quality of a model is closely related to its interpretability (Navarro et al. 2004).
Consequently, the WRC model fitting without a well-defined criterion might result in the selection of models that are not representative of that soil conditions. It is fundamental to accurately adopt a WRC that reflects the structural changes that result in different management systems. In this context, it is hypothesized that alterations caused due to soil management may influence the adjustment of the WRC model.
Obtaining water retention curves that best fit specific conditions of the soil, under different management systems, is fundamental for evaluating and understanding water dynamics. Thus, this study aims to investigate the best mathematical model suitable for WRC related to the management systems adopted, based on previously established criteria.

Location of the study area
The study area is located in Uruçuí (UTM: 526919.52, 9009171.95, SAD 1969), southwest region of the state of Piauí, Brazil ( Figure 1). The soil is classified as Ferralsol (or Oxisol according to Soil Survey Staff 2014 and Latossolo Amarelo Distrófico according to Embrapa 2018) (FAO 2015, presenting a caulinitic mineralogy and medium texture. According to the Köppen climate classification, the climate is Tropical Aw: hot and humid, with average precipitation of 1,100 mm year -1 and an annual mean temperature of 29 ºC (Peel et al. 2007).

Soil sampling and treatments
The water retention models "WRC" were fitted to nine areas with different management practices in a medium-textured caulinitic soil. The treatments were: area under native Cerrado (NC) regarded as the equilibrium condition, areas under eucalyptus farming system (EU6 and EU12, respectively, of six and twelve years of cultivation), areas under pasture (PA2 and PA6, respectively, of two and six years of planting), areas under conventional tillage system (CT2 and CT8, respectively, of two and eight years of cultivation), and areas under no-tillage system (NT3 and NT6 respectively, of three and six years of cultivation) (Table I).
Sampling was performed as follows: a point was delimited in the center of the area under cultivation. Subsequently, an area of one hectare was demarcated, using a GPS and a measuring tape. From this area, twenty-five points were defined twenty-five meters apart. A total of thirty-six sample units were collected from the systems in a medium-textured caulinitic soil, with four replicates. After field demarcation, four of the twenty-five points were randomly selected to compose the replication within each area.
Under native Cerrado, a border of 15 meters was established from the margin of the legal reserve, and the same procedure described above was applied (Silva et al. 2017).
Analyzes of organic carbon (OC), humic acid (HA), and fulvic acid (FA) Organic carbon was determined by the dry combustion method using an Organic Carbon Analyzer, Vario ®TOC Cube (Elemental brand). The fractionation of the humic substances was adapted from Benites et al. (2003), in which 1 g of air-dried soil were placed in contact with 10 mL of 0.1 molL -1 NaOH, the mixture was shaken for 30 s and kept still for 24 h. The samples were centrifuged at 18,100 g for 10 min, and the supernatants were stored for 1 h. The samples were centrifuged again and the supernatant was added to the previous extract. The precipitated fraction insoluble in alkaline medium (humin) was lyophilized.
Sulfuric acid (20%) was added to the supernatant obtained from the humin separation (pH ± 2), this acidified extract was used for humic acid fraction decanting process for 18 h. Subsequently, the supernatant was   Average fertilization applied in the planting rows was of 100 kg ha -1 P 2 O 5 and 120 kg ha -1 K 2 O in all soybean crop seasons, adding 130 kg ha -1 N in the corn cultivation year.

NT3
Area converted into farming in 1999/ 2000, with initial application, and every three years of 5,000 and 2,000 kg ha -1 limestone, respectively. In the first crop year, the area was grown with rice and then with soybeans, using millet as second crop or mulching in most years of cultivation by 2008/ 2009. In the 2009/ 2010 season, a no-till system was adopted (NTS) using millet for straw interspersing crops of corn and soybeans until the crop year of 2012/2013, with average fertilization using 100 and 120 kg ha -1 P 2 O 5 and K 2 O, respectively, adding 130 kg ha -1 P 2 O 5 N 2 for corn crops. In 2007/2008, no-till system was adopted using millet for straw. Average fertilization was performed using 100 kg ha -1 P 2 O 5 and 120 kg ha -1 K 2 O for all crop seasons of soybeans and corn, adding 130 kg ha -1 N during corn cultivation.
filtered through a membrane with 0.45-μm diameter pore, in a vacuum fi ltration system, to separate the fractions related to fulvic acid and humic acid.

Sulfuric attack and granulometry
The determination of soil granulometry was performed according to the pipette method (Embrapa 2011), and the granulometric fractions were classifi ed according to the international system (International Society of Soil Science), adopted by the Food and Agriculture Organization of the United Nations (FAO 2006) ( Figure 2). Extraction of silica (SiO 2 ), aluminum (Al 2 O 3 ), iron (Fe 2 O 3 ), and titanium (TiO 2 ) were carried out via sulfuric attack, according to Resende et al. (1987), for soil samples from the native Cerrado (Table II). The molecular ratios Ki (Eq.1) and Kr (Eq.2) were calculated according to the following equations: The Soil was classifi ed as caulinitic (Ki> 0.75, Kr> 0.75) (Table II), with medium texture (clay content varying between 150 and 350 g kg -1 ) ( Figure 2) (Fao, 2006, Embrapa, 2018.

Analysis of soil bulk density, porosity, and water retention curves
For the determination of soil bulk density (Bd), total porosity (TP), macroporosity (Ma), microporosity (Mi), and water retention curves (WRC), samples were collected with a preserved structure in the 0-0.05 m layer, four trenches were opened per area and depth, constituting four replicates in each soil management system. Samples were collected in volumetric rings (approximately 100 cm 3 ) and were saturated by capillarity action from the base.
The points for fi tting the WRC were obtained from the four replicates of each study. Nine matrix potentials were tested: -0.01 kPa as saturated soil condition; -2, -4, -6, and -10 kPa through tension tables; and -33, -100, -500 and -1,500 kPa, employing the Richards extractors (Embrapa, 2011). After reaching water equilibrium, the samples were weighed and dried at ± 105 °C for 24 h to defi ne Bd and the volumetric soil water content (Ɵ) linked to the levels of stress.
Subsequently, the aggregates that were retained in each sieve were dried at 110 °C for 24 h. The geometric mean diameter (GMD) was computed from the model proposed by Mazurak (1950), and the aggregate stability index (AEI) was calculated according to Castro Filho et al. (1998).

Models of water retention in soil
The use of soil water retention curve models under different management systems is presented in Table III: where is the volumetric water content (cm 3 cm -3 ), h is the soil water tension (kPa), and represent the saturated and residual water content respectively, c and d are dimensionless parameters. In the biexponential model (Omuto 2009), shows the difference between the saturated volumetric moisture and the residual volumetric moisture in the texture pore space at a given stress h, is the sum of structural residual moisture in the pore space and the textural pore space in the tension h, pore space and the textural pore space in the shows the inverse of the air-entry potential in the structural pore space, and symbolize the inverse of the air-entry potential in the pore space of textured-soil.

Criteria for the comparing adjustment water retention models
The nonlinear regression models are defi ned by the general formula The nonlinear regression models are defi ned by , where represents the observed value of volumetric water content, represents the observed value of volumetric is a nonlinear function that associates the value of to the soil water tension h (Table III), is the parameter vector  given in Table III, and is the error normal random variable (Gallant 1987).
The assumptions of normality of the error term with zero mean, homogeneity of variance term with zero mean, homogeneity of variance , and independence were verified by the Shapiro-Wilk, Durbin-Watson, and Bartlett tests, considering 5% as the nominal significance level (Fox 2016). The performance of the fi tted models was assessed considering three criteria, presented in the following sequence: The fi rst principle was the residual standard error (RSE) (Eq. 7) to evaluate the deviation between estimated and observed volumetric contents. When comparing two or more models, the one that exhibits the lowest RSE was preferred.
(Eq. 7) where  i is the observed volumetric water content in the i-th sampling unit and is the volumetric water content predicted by the model in the i-th sampling unit.
The second criterion of the adjustment is the coeffi cient of determination ( ) (Eq. 9). The use of the coeffi cient of determination (R 2 ) (Eq. 8) does not contemplate the number of parameters of the model, which can expand it. Thus, a correction was performed on R 2 , adding the mathematical term (Eq. 7) as a penalty factor. Therefore, the equation that showed the largest as a penalty factor. Therefore, the equation that was considered more appropriate (Archontoulis & Miguez 2013) and was obtained as follows: (Eq.8) (Eq. 9) The Bayesian information criterion (BIC) (Eq. 10) (Schwarz 1978, Emiliano et al. 2014, which estimates the quality of the fi tted model, is based on the maximum likelihood function (MLF) that is dependent on the number of observations and parameters. The model that presented the lowest BIC was chosen. BIC is defi ned as: (Eq. 10) In the equations (Eq. 7, Eq. 8, Eq. 9 and Eq. 10): p is the number of parameters of the adjusted model, is the number of parameters of the is the predicted value of Table III. Models of water retention in soil used for different management in a medium-textured kaolinitic soil.

Model of water retention in the representative soil Parameters
van Genuchten-4P (1980) the volumetric water content, is the MLF evaluated on the parameters that return the maximum of the MLF. This function is based on the product of n functions such as those in Table  III, where n is the number of experimental data and ln is the natural logarithm. More details on this process can be seen in Ritz & Streibig (2009).
In the comparison of two or more models for water retention, the one that showed more accuracy in most of the established criteria was considered better adjusted, as the predicted values demonstrated greater agreement with the observed values.

Principal Component Analysis (PCA)
Because the variables present measurements in very different units, prior to sample agglomeration (Eq. 11) units, the data were standardized according to Cao et al. (1999) (Eq. 11). (Eq.11) where is the i th transformed variable, is the value of the ith variable (treatment) before it was transformed, is the mean of the ith variable, and DP i is the standard deviation for i th variable in all scenarios under study.
The principal component analysis (PCA) allowed the correlation between diverse characteristics in each management system and the adjustment of the WRC, characterizing the variables that better discriminated the distinguishing structure in each soil system. Consequently, the initial set of seventeen variables was characterized by two new orthogonal latent variables, which admits the location in two-dimensional fi gures (Hair Junior et al. 2009).

Computational resources
The parameters of the models were estimated by the least squares method for nonlinear models (Table III) (Ritz & Streibig 2009) through the nls function in R (R Core Team 2018) and the initial values, required to estimate the parameter, were obtained using the iterative graphical method in non-linear regression using the manipulate package (Allaire 2014). The manipulate function accepts a plotting expression and a set of controls (e.g. slider, picker, or checkbox) which are used to dynamically change values within the expression, in our case the models in Table III. When a value is changed using its corresponding control the expression is automatically reexecuted and the plot is redrawn. So in this way it is possible to obtain initial values for the parameters of the models that return the curve as close as possible to the experimental data dynamically and these initial values are then used in the estimation process.

RESULTS AND DISCUSSION
From different criteria (Table IV), it was observed that the fitted WRC, in general, presented differences regarding the type of adjustment for different soil management systems. However, the best performance for each system was obtained by the Fredlung-Xing (FX) and van Genuchten-4P (VG4P) models, whereas the other models showed worst performance in all circumstances (Table IV).
The results showed that FX and VG4P present lower BIC and EPR values and higher values were associated with these models. For the adjustment of six of the nine management systems, the FX model was superior, in which the values of BIC, EPR, and systems, the FX model was superior, in which ranged from -65.22 to -54.42, 0.005 to 0.008, and 0.985 to 0.998, respectively (Table IV). The VG4P model showed superiority regarding the adjustment  (1) Soil managements: native Cerrado (NC); eucalyptus plantation of six years (EU6); eucalyptus plantation of twelve years (EU12); pasture of two years (PA2); pasture of six years (PA6); conventional tillage of two years (CT2); conventional tillage of eight years (CT8), no-tillage of three years (NT3), no-tillage of six years (NT6). *P-value for the adequacy tests of the models: (2) Bartlett Test, of three of the nine management systems, in which the values of BIC, EPR, and of three of the nine management systems, in ranged from -56.10 to 41.59, 0.008 to 0.018, and 0.980 to 0.994, respectively (Table IV). It should be noted that the adjustment criteria vary according to the management and the model, which suggests that the WRC should be compared before they are effectively used in each soil management system.
The type of management (Table IV) influenced the adjustment. The model VG4P, proposed by van Genuchten (1980), exhibited a greater adjustment for the CT2, CT8, and EU12 systems of WRC. The FX model (Fredlung & Xing 1994) provided a better fi t for the systems EU6, PA2, PA6, NT3, and NT6, in addition to the reference condition (Native Cerrado -NC) ( Table  IV).
Due to a smaller number of parameters, models such as VG4P showed an enhanced adjustment potential, which translate into a parsimonious model (Omuto 2009). Studies involving WRC have demonstrated that few adjustable parameters are considered a prerequisite to obtain the best-adjusted model (Campbell 1974, Groenevelt & Grant 2004, Too et al. 2014).
The FX model was better fitted for the reference condition (RC) and for the management systems that prioritize soil cover maintenance (EU6, PA2, PA6, NT3, and NT6) ( Table IV). The success in adjusting WRC is related to the fact that this model was developed as a function of the soil pore size distribution (Fredlung & Xing 1994), discerning the pore texture and structure in WRC models (Durner 1994, Fredlund & Xing 1994, Omuto 2009).
Thus, these results showed that the WRC for medium-textured soils are sensitive to the adopted management. Thus, it is important to select models based on scientific principles (Burnham & Anderson 2004), considering soil physics, instead of applying equations in different situations, even without comparisons to other models (Barros et al. 2013, Too et al. 2014. After analysis the main components PC1 and PC2, it we selected the ones with eigenvalues of 2.99 and 1.87 respectively, and a cumulative variance of more than 73% (PC1 = 52.65% and PC2 = 20.70%), factor that confer appropriateness to the analysis (Hair Junior et al. 2009).
In Figure 3 it can be observed that the management systems showed a grouping tendency, in which discriminatory variables led the creation of groups I and II. In this regard, the pore volume (PV), the larger diameter pores (Ma), aggregation indexes (GMD and ISA), larger aggregates (AG1), organic carbon (OC), and the more stable fractions of organic matter (HA) are responsible for the discrimination of group I, where the management systems best fi tted the FX model (NC, EU6, PA2, PA6, NT3 and NT3) (Table IV and Figure 3). The smaller aggregates fractions (AG4, AG5, and AG6) are the most relevant factors in the formation of group II (EU12, CT2, and CT8), which is constituted of improved management systems better suited to the VG4P model (Table  IV and Figure 3).
These results illustrate that the NC, EU6, PA2, PA6, NT3, and NT3 systems promote soil structure building, which may be a consequence of non-tilting and the addition of organic residues on the surface (Castro Filho et al. 1998). On the other hand, the EU12, CT2, and CT8 systems cause a greater fractionating of soil aggregates, which is related to crop rotation (Oliveira et al. 2003) and to lower OC content (Figure 3). In this context, Dufranc et al. (2004) reported that OC contributes to the formation and stabilization of soil aggregates, an energy source for microorganisms, important aggregation agents, which directly impact aggregation indices.
The analysis of the principal components identifies a relationship between the soil structural condition and the adjustment of water retention models. In this respect, management systems with more desirable conditions ( Figure   3) better fit the FX model (Table IV). However, the WRC was efficiently adjusted employing the VG4P model, in which management systems provided greater fractionation (Figure 3 and Table IV).
At noticeable variation in moisture levels in the soil, the FX model is more efficient in modeling WRC, which caused a steep slope of the curve. The VG4P model, however, seems suitable under homogeneous pore distribution of WRC, displaying a smooth slope. These results suggest that the FX model is more efficient for predicting the low suction potential of water, retained in larger pores while the VG4P model is more efficient under high suction potential in poorly structured soil conditions. Nevertheless,  a detailed analysis of the residuals to identify specifi c points of the adjustment is required.
Overall, the analysis of the residues indicated that the FX model showed smaller residues at lower suction values (0.01, 2, 4, and 6 kPa), while the VG4P models were capable of estimating the humidity at higher suction values (100, 500, and 1500 kPa) (Figure 4 and Figure 5). It is important to mention that well-adjusted models hold small residuals, since higher values are indicative of a poorly adjusted data representation, as shown in Figure 4.
Particle size distribution, structure, and mineralogy also interfered in the shape of WRC (Omuto 2009). Relatively large pores in well-structured soils cause abrupt variations when fading away. In bad-structured soils the pore distribution is more uniform, so that with the respective matric suction, it drains a portion of the pores and a certain amount of water remains, which reduces the variation in the curve.
Porosity is especially important to predict soil WRC behavior since it is the most sensitive to operation. Most models ignore structural pores, even though porosity under field conditions is composed of textured and structural pores (Dexter et al. 2008, Too et al. 2014, resulting in an inadequate representation of WRC (Durner 1994).

CONCLUSIONS
Soil management influences the selection of the water retention curve model. Being its greatest influence, observed near the soil saturation zone.
The adjustment criteria for nonlinear models were efficient in choosing the most suitable model for each management system.
The Fredlung-Xing model was more efficient for adjusting the water retention curve of the medium-textured caulinitic soil under native Cerrado and conservationist management systems.
The four-parameter Van Genuchten model was more efficient in modeling the water retention curve of soil under management systems that provide fractionation of soil aggregates.