Estimation of soil hydraulic and solute transport properties from field infiltration experiments

To model water flow and solute transport in soils, hydrodynamic and hydrodispersive parameters are required as input data in the mathematical models. This work aims to estimate the soil hydraulic and solute transport properties using a ponded axisymmetric infiltration experiment using a single-ring infiltrometer along with a conservative tracer (Cl) in the field. Single ring infiltration experiments were accomplished on an Oxisol in Areia in the state of Paraíba, Brazil (6o 58' S, 35o 41' W, and 645 m), in a 50 x 50 m grid (16 points). The unsaturated hydraulic conductivity (K) and the sorptivity (S) were estimated for short or long time analysis of cumulative three-dimensional infiltration. The single tracer technique was used to calculate mobile water fraction (Ф) by measuring the solute concentration underneath the ring infiltrometer at the end of the infiltration. Two solute transfer numerical models based on the mobile-immobile water concept were used. The mobile water fraction (Ф), the dispersion coefficient (D), and the mass transfer coefficient () between mobile and immobile were estimated from both the measured infiltration depth and the Cl concentration profile underneath the infiltrometer. The classical convection-dispersion (CD) and the mobile-immobile (MIM) models showed a good agreement between calculated and experimental values. However, the lowest standard errors to the fitted parameters were obtained by the CD model.


Introduction
In many applications dealing with environmental sciences, the knowledge of soil hydraulic properties is essential: (1) to diagnose the hydrodynamic functioning of soils in relation to the natural and/or anthropogenic constraints which are applied, and (2) to simulate the physical processes to establish a prognosis on the order of magnitude of the hydraulic fluxes able, for instance, to make water and nutrients available for the plant rooting system, or to advect chemicals leading to point or diffuse pollution of the groundwater table. Physically-based modeling of coupled water and solute transport requires the knowledge of soil properties that have to be estimated with adequate accuracy. In this sense, the development of in situ methods to determine both the saturated and unsaturated hydraulic properties of the soil surface has received increasing attention in recent years, in particular, to assess the role of preferential flow pathways in soil management practices (Köhne, Schlüter, & Vogel, 2011;Legout, Legout, Nys, & Dambrine, 2009;van der Linden, Tordesillas, & Narsilio, 2019;Zhang, Cao, Hou, & Cheng, 2021).
Laboratory and field studies have been conducted for the last three decades to estimate the vadose zone's hydraulic and chemical transport properties. But the field transport studies are more complex; the field soils are usually structured by containing large continuous macropores, such as inter-aggregate pores, earthworm burrows, drying cracks, or decayed root channels (Al-Jabri et al. 2006;Arora et al., 2019;Köhne et al., 2011;Tabarzad, Sepaskhah, & Farnoud, 2011). These macropores are generally characterized by distinctly different hydraulic properties from the soil matrix and may result in significantly more rapid solute movement through the unsaturated zone as indicated by average flow estimations, due to preferential flow (Ersahin et al., 2002; Gerke & Maximilian Köhne, 2004;Kamra & Lennartz, 2005;Li, Yao, Yan, & Cheng, 2021).
Complex theories govern the dynamics of solutes and the processes of transformations to which they are subject. For this reason, more and more solutions are proposed for problems concerning the distribution of nutrients in soil solute as to its salt leaching and the pollution of the groundwater by the use of toxic chemical products, pesticide residues, etc. The hydrodispersive characteristic of the soil is necessary to understand the mechanisms of the unequal movement of chemical substances through a soil profile under natural conditions on field study. The destination and movement of these dissolved substances in different soils and aquifers have caused considerable interest in the quality of their unsaturated zone. However, it is challenging to obtain values from the parameters of transport having direct physical signification, such as the water flow in the pores, the retardation factor, the coefficient of dispersion, and its degradation or production. In this sense, the spatial and temporal scales have an essential role in determining these values (Selim, Persson, & Olsson, 2017;Wallis & Manson, 2019).
The reliable accomplishment in predicting the destination and transportation of the chemical substances in the unsaturated zone of the soil resides in the capacity of accurately determining the parameters of their transport.
Tension disk infiltration methods for characterizing hydraulic and chemical transport properties are often timeconsuming, and large areas cannot be easily sampled. Nevertheless, it was successfully used to estimate simultaneously unsaturated hydraulic and solute transfer properties in soils (Angulo-Jaramillo et al., 1997;Angulo-Jaramillo, Gaudet, Thony, & Vauclin, 1996;Clothier, Kirkham, & McLean, 1992;Netto et al., 2013;Roulier et al., 2002). The method herein presented is a simple, easy, and cheap one to implement, known as Beerkan Lassabatère et al., 2006;Mubarak et al., 2010;Souza, Antonino, Angulo-Jaramillo, & Netto, 2008;Yilmaz, 2010), pioneered by Haverkamp et al. (1994), and relies on particle-size analysis, dry bulk density and single ring axisymmetric infiltration tests. The experimental protocol and the method of data analysis, leading to the estimation of parameters describing hydraulic properties, are herein described. This method can also be used to estimate the solute transport properties, known as Beerkan-Solute (Netto et al., 2013), when the water reaches the steady-state infiltration rate, and then three or four equal volumes of solute tracer are applied.
One objective of this study consists in the hydrodynamic and hydrodispersive characterization of an Oxisol, by combining the measures obtained in the field with those of a numerical model (Beerkan-Solute method), from an easy and efficient way of axisymmetric infiltration by using a single ring infiltrometer with a conservative tracer (Cl-).

Studied site and soil characterization
This study was conducted in a plateau (619 meters above sea level) in Areia, a county of the state of Paraíba, northeast of Brazil (6 o 58'24''S, 35 o 41'59,8''W). The annual average precipitation rate in the region is 1450 mm. The highest precipitation rates in the region are usually observed during June and July.
The soil studied was a sandy clay loam (Oxisol). The soil characterization was based on physical analysis of 80 samples collected in a grid 25x25 m 2 long. Ten samples were used for the chemical analysis. The soil profile samples were taken from a 0 to 15 cm depth. The mean and standard deviation of soil dry bulk density from the 80 points were 1.126 and 0.086 g cm -3 , respectively. The mean and standard deviation for the sand, silt, and clay percentages were 60.8 ± 4.2, 14.3 ± 2.3, and 24.9 ± 4.6, respectively. This soil is designed as a well-structured and very well-drained one. The aggregate-mean size of the 0-15 cm layer is 1.76 mm in diameter when dried and 1.36 mm when wetted. The ratio between the wet and the dry aggregate size is 0.77, indicating high aggregates stability.
As shown in Table 1, this soil is deficient in nutrients (Ca, Mg, Na, and P) and is very acidic, despite the high organic matter and clay contents. The CEC value reflects that this soil, when subjected to natural acidity conditions, presents a low capacity to retain cations and a low activity of the clays.

Infiltration Experiments
Infiltration experiments using water and solute were performed in a 4 ha cropped area, employing a 50x50 m 2 grid, with 10 sampling points ( Figure 1). Tillage and sowing (beans (Vigna Unguinculata (L.)) were carried out in March 2006, and the measurements were conducted in the early development stage of the crop. Infiltration experiments were performed using a ring infiltrometer with a diameter of 15 cm, water, and a Chloride (KCl) solution as a soil-water tracer. The experiments were performed by installing the device just below the soil surface to minimize the disturbance in the soil structure and ensure an axisymmetric water flow on the soil surface. Source: Authors.
The initial volumetric water content, ini (cm 3 cm -3 ), initial chloride concentration, Cini (mg L -1 ), and dry bulk density, d (g cm -3 ), were measured from undisturbed soil samples taken a few centimeters from the infiltration location. Cores were taken by using a cylinder 50 mm high and 795.2 cm 3 in volume for three soil layers corresponding to 0-5, 5-10, and 10-15 cm in depth.
The infiltration experiments were conducted following the Beerkan method (Lassabatère et al., 2006;Lassabatère et al., 2019;Mubarak et al., 2010). Fixed water volumes were poured into the ring, and the time required for the infiltration of each volume was measured until steady-state. Volumes of water ranging from 70 to 250 mL were infiltrated. Different water volumes were chosen to avoid a high positive pressure head at the soil surface and thus prevent an induced water flow. When the steadystate infiltration rate was reached, three equal volumes (100 mL) of KCl (0.1M) solution were applied. From the geometric configuration, it is then assumed a three-dimensional infiltration at null pressure head, i.e., natural saturated volumetric water content on soil surface. In contrast to tension disc infiltrometer experiments ( Clothier et al., 1992;Roulier et al., 2002), no surface preparation is required, so the soil surface is not disturbed. In the other grid points (Figure 1), the infiltration experiments were performed with water only.
At the end of the solute infiltration, the soil was sampled at a depth of 15 cm by using a sampler device that allows obtaining fifteen undisturbed soil samples one centimeter high and 15.9 cm 3 in volume along the vertical axis of the ring infiltrometer. These samples were used to determine the final profiles of both volumetric water content f (cm 3 cm -3 ) and total resident chloride concentration Cr (mg L -1 ). Table 2 shows the soil textural classes of all experimental points, while Table 3 presents the soil dry bulk density and porosity for the experimental points at three different depths. It is observed that the soil dry bulk density shows a tendency of increase with depth.

Estimation of the shape parameters by analysis of F(D)
The shape parameters n and m from van Genuchten´s equation (1980) can be obtained from the particle size distribution F(D), assuming that the porous radius is inversely related with the soil water pressure and a similarity between the cumulated porosity size distribution S(R) and F(D) exists, and consequently between F(D) and h(). Haverkamp and Parlange (1986)  (1) where F(D) is the particle size distribution, D is the effective diameter of one particle of the soil [L], Dg is the scale parameter of the particle size [L], and M and N are the shape parameters for the particle size distribution. The relationship between M and The parameters M, N, and Dg were obtained by fitting the experimental particle size distribution data.
In order to obtain the shape parameters m and n of the soil-water retention curve h() given by van Genuchten (1980), the following relationship was proposed (Fuentes, Vauclin, Parlange, & Haverkamp, 1998): ( 2) and (3) ( 4) where  is a coefficient defined by Fuentes et al. (1998), and s is the relative fractal dimension. Knowing that s = Df/E, where Df is the soil fractal dimension, and E = 3 is the Euclidian dimension, the dependence of s regarding the global soil porosity () is implicitly defined as follows: For the hydraulic conductivity equation of Brooks and Corey (1964), the shape parameter  can be expressed as a function of the soil-water retention parameters (for instance,  and m) and the tortuosity factor (p), by: with p =1 (Burdine, 1953).

Estimation of hydraulic properties by Beerkan method
The unsaturated hydraulic parameters of the soil-water retention curve, h() and soil hydraulic conductivity, K() were estimated using the BEST_1.3.0 program (Excel version) (Di Prima, 2013) by the following explicit transient two-term (Eq. (7a) and (7b)) and steady-state expansions (Eq. [8a] and [8b]) ): a short time and long time analysis of cumulative three-dimensional infiltration equations : where constants A B, and C can be defined by : , , and where S is the sorptivity (L T -1/2 );  is a constant parameter of 0.75 ;  is a constant shape equal to 0.6; r is the radius of the ring (L); Ks and K 0 are the hydraulic conductivity values corresponding to s (the volumetric water content at the end of the infiltration experiment), and  0 (the volumetric water content at the beginning of the experiment).
Hydraulic parameters Ks and S are estimated considering the initial conditions by the BEST algorithm (Beerkan Estimation of Soil Transfer parameters; Lassabatère et al., 2006). The associated standard errors Ks (mm s -1 ) and S (mm s -1/2 ) were calculated considering residues between the fitted infiltration equations, Eqs. (7a, 7b) and Eqs. (8a, 8b), for the experimental cumulative infiltration curves.

Mobile water fraction estimation
The mobile water fraction () is the ratio between the mobile volumetric water content, m, and the total volumetric water content () and was calculated from the total solute concentration at the end of solute infiltration experiments, C * , following the procedure of Clothier et al. (1992) and Angulo-Jaramillo et al. (1996). Assuming that there is no solute exchange between the mobile and immobile water regions at the sampling time and depth (i.e., Cim(z0,tsol) ≈ 0),  can be written as: (12) where Cini, was the initial solute concentration in the soil, and C0 was the concentration of the input solution.

Estimation of solute transport parameters
Concerning the solute transport parameterization, the experimental data were analyzed by the classical convectiondispersion (CD) model and the mobile-immobile (MIM) model, using the CXTFIT 2.0 numerical code (Toride, Leij, & van Genuchten, 1995) for resident concentration. One-dimensional steady-state flow is assumed to occur at the center of the infiltration ring during the solute application time.
The one-dimensional convection-dispersion (CD) equation for a non-sorbing conservative solute in a homogeneous porous media with no source and/or sink, and under steady-state flow conditions, is written as: (13) where C is the relative concentration, T is the pore volume number, X is the relative distance, P=Lv/D is the Peclet number, and R is the retardation factor.
The parameters D and R were estimated by using the CD model, following the condition of mass conservation.
Assuming that for small soil volumes , the following relation can be observed between the retardation factor, R, and the mobile water fraction, Ф, determined according to the method of Clothier et al. (1992), equation (12): The initial value for the pore water velocity v was estimated by: Research, Society and Development, v. 10, n. 14, e195101421764, 2021 (CC BY 4.0) | ISSN 2525-3409 | DOI: http://dx.doi.org/10.33448/rsd-v10i14.21764 where q1Dtsol is the flow density at the end of the infiltration experiment and is the average volumetric water content from the measurements along with the soil profile at the end of the infiltration experiment.
The dimensionless two-region solute transport equation (MIM) is given by: where the subscripts 1 and 2 refer to the equilibrium and nonequilibrium sites, respectively; ' is a partitioning coefficient, and  is a dimensionless mass transfer coefficient (i.e., the Damkhöler number).
The dimensionless parameters may be defined as shown in Table 4, where:  = q/ is the effective pore-water velocity (L T -1 ); q is the Darcy's velocity (L T -1 ); L is the length of the flow system (L); z is the distance (L); t is the time (
The MIM model was used to fit the parameters D, R, and ω. The initial value for v was estimated according to Eq. (15).
The parameter β' was taken as equal to the measured water mobile fraction, Ф (Eq. 12), and was taken as constant because solute is non-reactive. The mass conservation condition was used in the fitting process for both CD and MIM models.
The one-dimensional flux densities are evaluated by: for short-time (Eq. 17) and long-time (Eq. 18) infiltration equations, respectively.

Results and Discussion
Cumulative infiltration data and the fitted curves using the long-time cumulative infiltration equation for four experimental points are depicted in Figure 2. The experimental point C8, with the lowest porosity, presented the highest infiltrability.
Research, Society and Development, v. 10, n. 14, e195101421764, 2021 (CC BY 4.0) | ISSN 2525-3409 | DOI: http://dx.doi.org/10.33448/rsd-v10i14.21764  Table 5 presents the values of the final volumetric water content as well as the variation of the volumetric water content between the end and the beginning of each experiment.
The relatively high values of Ks may be associated with the influence of high clay (Table 2) and organic matter content (Table 1), which may contribute to the formation of aggregates in the soil leading to the occurrence of macropores. Seyfried and Rao (1987) and Allaire et al. (2002bAllaire et al. ( , 2002a highlighted that high clay content associated with high hydraulic conductivity values would frequently indicate flow through macropores. fin -final volumetric water content; Δ − variation of the volumetric water content; Ssorptivity; Kssaturated hydraulic conductivity; hg -is related to the air entry suction; n and mshape parameters of the soil-water retention curve;  − shape parameter of the hydraulic conductivity curve; cp -a parameter which depends only on the shape parameters n, m and  tmax -maximum time.

Source: Authors.
Four different behaviors may be distinguished from the results of the hydrodynamic parameters (Table 5): 1) the experimental points A6, C6, and C8, with high sorptivity and high hydraulic conductivity; 2) the experimental points A4 and C2, with high sorptivity values and relatively low hydraulic conductivity; 3) the experimental points A2, A8, C4, and E2, exhibiting low sorptivity and low hydraulic conductivity; and, 4) the experimental point E4, with the lowest sorptivity value in ten experimental points, also presenting lower hydraulic conductivity. The four different characteristic behaviors may be attributed to differences in the soil structure characteristics at the topsoil layer (0-5 cm) and in the compaction of the lower layers (5-10 cm and 10-15 cm). As indicated by the sorptivity values, the initial conditions at the top layer refer to the influence of climate (rainfall) and machinery and man movement on the soil surface, which may lead to crust formation controlling soil infiltrability. The same factors may cause compaction of the lower layers, thus influencing the hydraulic conductivity values.
The different estimated hydrodynamic characteristics reflect the effect of the soil heterogeneity at the study site, as influenced by natural and soil use conditions. Based on Table 5, the variation of volumetric water content (Δ) was 0.117 to 0.214 cm 3 cm -

3
. For all experimental points, it was observed that the geometry of the source was not dominant on the initial character of the capillary in the infiltration process.
The estimates of Ks and S allowed the calculations of the one-dimensional flux density, q1Dtsol, and the pore water velocity (Eq. 15), hereafter called calculated velocity (vcal). The values of the mobile water fraction, Ф, calculated through the solute concentration at the end of the experiment (Eq. 12) are presented in Table 6, along with the calculated velocity and q1Dtsol. Ciniinitial solute concentration; q1Dtsol -flow density at the end of the infiltration experiment;  = m/ − mobile water fraction, vcal -effective pore-water velocity; Isolcumulative solute infiltration; Zfront = Isol/ -wetting front. Source: Authors.
The values of the mobile water fraction varied from 0.42 to 0.85 within the study site, with the highest value related to point C8, which showed the highest infiltrability. The values of  obtained (Table 6) are in accordance with the values obtained in the existing literature. In the field experiments, Clothier et al. (1992) found values of  up to 0.49 ( 0,12), by applying bromide (Br -) in fine sand. Gvirtzman and Magaritz (1986) found  values ranging from 0.45 to 0.60 in an experiment with a natural tracer. Rice, Bowman, and Jaynes (1986) found average  values of 0.80 in soil under irrigation. Jaynes, Logsdon, and Horton (1995), also using bromide, found values of  ranging from 0.25 to 0.98. Angulo-Jaramillo et al. (1996), using Cland Roulier et al. (2002, using 18 O as a tracer, found  values ranging from 0.5 to 0.93.
For the point E4, with the lowest infiltrability, the lowest sorptivity, and low hydraulic conductivity, which may be associated with the presence of a surface crust and compaction of the underlying soil, but having the mobile water fraction relatively low, indicating that only part of the entire pore volume was functioning during transfer. Point C8, with high hydraulic conductivity and sorptivity values, presented the highest mobile water fraction among the ten experimental points, indicating that even with smaller pore size, the pores were more connected in a soil where surface crust and compaction might not occur.
For point C2, with high sorptivity, a surface crust may not occur. However, with relatively low hydraulic conductivity, the relatively low value of the mobile water fraction was found. Table 7 shows the results of the hydrodispersive parameters estimated by using the CXTFIT 2.0 program with 95% confidence interval, considering the CD model. The Peclet number shows that the predominating process in the transportation Research, Society andDevelopment, v. 10, n. 14, e195101421764, 2021 (CC BY 4.0) | ISSN 2525-3409 | DOI: http://dx.doi.org/10.33448/rsd-v10i14.21764 13 of chloride in five points (A2, A8, C4, C6, and E4) is the diffusion (P < 1). About the effects of the aggregates in the transportation of the solute, while undergoing experiments with colorful tracers, Seyfried and Rao (1987) showed that only a fraction of the soil matrix was in contact with the solute, concluding that the diffusion within the soil matrix would be the means of the mass transfer to these systems. The coefficient of hydrodynamic dispersion (D) ranged from 64.3 to 1950.0 cm 2 h -1 , thus pointing out the high spatial variability of this parameter. These values can be considered high when compared to other studies. For example, Roulier (1999) found values up to 7.8 cm 2 h -1 in sandy soil, using the tension disc infiltrometer. Regarding soil column studies or similar (glass front Dtsol Z ω q = α  1 spheres, for example), Rao et al. (1980a) found values ranging from 1.57 to 12.46 cm 2 h -1 . By using soil with aggregates, Seyfried and Rao (1987) found values ranging from 1.65 to 60.9 cm 2 h -1 ; Bejat et al. (2000) found values of 6.64 cm 2 h -1 for a silty-clay soil and of 20.8 cm 2 h -1 for a sandy-clay soil; Comegna, Coppola, and Sommella (2001) found values of about 9.54 cm 2 h -1 in a sandy-clay soil. According to Renard et al. (1977), the hydrodynamic dispersion coefficient is one of the most sensitive parameters to measurement errors.

CD Model
Regarding dispersivity ( = D/v), the condition imposed by Snow (1999), Isol > , for the validation of the hypothesis Cm(zo,tsol) = Co, was considered for points C2 (Isol = 1.53 cm;  = 1.5 cm) and C8 (Isol = 1.73 cm;  = 0.4 cm). However, Clothier et al. (1995) showed that a water layer of up to 15 mm could be used to validate the hypothesis Cm(zo,tsol) = Co. Moreover, during the soil-water infiltration experiment, there was a concern to avoid a possible hydraulic head when applying a large water layer.
This would invalidate the infiltration equations since they refer to an infiltration without a hydraulic head. By applying the method of Clothier et al. (1992), for the determination of mobile water fraction (), Roulier (1999) also found values of a water layer of infiltrated solution (Isol) smaller than dispersivity ().
The values of Isol presented in Table 6 were calculated using the short-time infiltration equation and the fitted values of S and Ks. The values in Table 6 show that the conditions posed by Snow (1999) to meet the hypothesis of solute-free immobile water content can be applied for the experimental points.
The root mean square errors (RMSE) of the fitted parameters using both models and the Peclet number are presented in   The fitted values of the dispersion coefficient, D, were very high in some cases compared to values reported in the literature (i.e., Rao et al., 1980a;Seyfried and Rao, 1987;Roulier, 1999;Beja et al., 2000;Comegna, Coppola, and Sommella, 2001). According to the estimated Peclet number, the dispersion was the dominant effect in the solute transport.

Conclusion
This study applied a combined analysis of water infiltration, and solute movement in an Oxisol cultivated with beans (Vigna Unguinculata (L.)) to determine the hydrodynamic and hydrodispersive parameters.
A simple ring infiltrometer showed to be a simple and efficient method to obtain the hydrodynamic and hydrodispersive parameters.
The hydrodynamic characteristics (Ks, S, hg and characteristic time and length) were obtained from the threedimensional accumulated infiltration equations for a long time. The determined parameters permitted a good estimation of the one-dimensional fluxes and, also, the assessment of the average velocities applied in estimating the hydrodispersive parameters.
A good agreement between calculated and measured values was achieved for the CD and MIM models. However, the lowest standard errors to the fitted parameters were obtained by the CD model. The estimate of the dispersion coefficient (D) values from MIM model were underestimated when compared to those obtained from the CD model.
In any case, due to the general poor definition of parameters involved in the MIM model, it is very arduous and hazardous to attach them to any physical meaning. Apparently, the increased number of parameters invoked to represent additional processes produced unreliable estimates and unproved goodness of fit when nothing clearly corroborated the existence of those processes. When several parameters are optimized simultaneously, the problem of non-uniqueness arises because of the correlation between parameters. However, one should recall that when the number of parameters decreases, the optimization technique allows us to identify unique parameter values, and, in all those cases, the estimated parameter set is strong enough for the considered prediction. It must be recognized that the application of the dual-porosity model to many field systems will be an approximation of reality, although it is a useful one.