Hydrodynamic dispersion of sodium by equilibrium and non-equilibrium methods through undisturbed sandy soil columns from the Adamantina Formation

Column tests using undisturbed samples of residual sandy soil from the Adamantina Formation (K) were performed to determine the sodium hydrodynamic coefficient by equilibrium and non-equilibrium methods. Analyzing the porosity and soil characteristics, as well as the breakthrough curves, it was possible to note a typical non-equilibrium transport behavior, related to the soil dualporosity, particularly at early and late time behavior in the breakthrough curve tails. The experimental data were best fitted to the non-equilibrium method than to the equilibrium one, and the sum of the square errors was up to five-fold lower. The conceptual model analysis of the methods led to observations on the limitations of each method and a careful analysis of the obtained values, which are also related to the soil characteristics.


Introduction
The hydrodynamic dispersion coefficient (Dh) is one of the most important solute transport parameters, and its correct quantification is essential to accurately estimate groundwater contamination and, therefore, analyze, assess and manage the possible risks related to cont aminated sites.Typically, the experimental Dh is obtained inversely, by fitting the analytical solution of the differential equation of a transport model, and the convection-dispersion equation (CDE), for specified boundary conditions (Lapidus and Amundson, 1952) to the experimental breakthrough curves (BTC).In addition, Dh can be determined directly from column tests, using one of the available methods (Bear, 1961;Ogata & Banks, 1961;Brigham, 1974).
When column tests are performed in soil columns with unpreserved intact structures, these methods work very well (Carmo, Antonino, Netto, & Corrêa, 2010;Magga et al., 2008) since they were developed empirically from tests in homogeneous soils and, often, non-reactive solutes.In these conditions, the equilibrium advective-dispersive transport model can represent the regular behavior of the solute.
However, those models are based on simplifications that may not accurately reflect natural conditions.It is well known that undisturbed soil columns are subject to the presence of macroporosity, dual-porosity, and other heterogeneities.These characteristics may strongly affect solute transport parameter values by creating non-uniform flow fields with widely different velocities, phenomena often referred to as nonequilibrium transport, i.e. when an immobile (related with micropores) and a mobile (related with macropores) liquid phase is present (van Genuchten & Wierenga, 1976;Brusseau & Rao, 1990;Ray, Ellsworth, Valocchi, & Boast, 1997;Schwartz, Juo, & McInnes, 2000;Lamy, Lassabatere, Bechet, & Andrieu, 2009).Usually, for non-equilibrium transport, the CDE is valid in the mobile domain, and diffusion is assumed in the immobile one (van Genuchten and Wierenga, 1976).
In this context, this study applied equilibrium and non-equilibrium methods to determine the Dh of the sodium ion, comparing both.To achieve this aim, column tests were performed using undisturbed columns of residual sandy soil from the Adamantina Formation (K) in the city of Cabrália Paulista, in the state of São Paulo, Brazil.Five methods that consider equilibrium advective-dispersive transport model and one that considers a non-equilibrium model were used.Our emphasis is on examining the characteristics of each method and the errors resulting from the use of a less suitable method for the studied soil, demonstrating the importance of choosing a correct model for determining the Dh, and the complexity of this choice.

Material and methods
This study was carried out through the following basic steps: 1) Soil sampling; 2) Soil physical-chemical characterization; 3) Column tests and porosity analysis; 5) Chemical analysis; and 6) Determination of the hydrodynamic dispersion coefficients.These six steps will be detailed after the basic characteristics of the study area section.

Basic characteristics of the study area
The study area is located in Cabrália Paulista city, in the state of São Paulo, Brazil.The main pedological soil types identified therein are Oxisol, Ultisol, and Latosol.This area is constituted by the residual soil of sandstones from the Adamantina Formation, Bauru Group (K).The area was chosen because its geological materials are spread at about 117,000 km 2 in the state of São Paulo (Silva, Kiang, & Caetano-Chang, 2003) and are subdivided into several local aquifer units, which are very important in terms of public water supply.

Soil sampling
Disturbed samples were collected from the study area to characterize soil mineralogy and physical, chemical and physicochemical properties.Undisturbed soil columns were carefully collected from handexcavated trenches, using rigid PVC cylinders, for column testing.

Soil physical-chemical characterization
In the laboratory, the disturbed samples were air-dried and sieved through a #10 mesh sieve (2 mm openings).The mineralogical analysis was performed by X-ray diffraction by the powder method in the clay fraction (Azaroff and Buerger, 1953), while the particle size analysis was performed according to the ASTM D422 standard and the density of the solids was determined according to the ASTM D854-10.Subsequently, the following parameters were also determined: pH in H2O and in KCl, Eh and electrical conductivity (EC) (Donagema and Campos, 2011), delta pH (pHKCl -pHH2O) (Mekaru and Uehara, 1972), point of zero charge (PZC) (2pHKCl -pHH2O) (Keng and Uehara, 1974), organic matter content according to the ASTM D2974 -00 and cation exchange capacity (CEC).A porosity analysis by mercury intrusion was also performed (Washburn, 1921).

Column tests
The PVC cylinders used for collecting the undisturbed soil samples (150 mm long and with 97.2 mm internal diameter) were used as rigid-wall permeameters, in which column tests were conducted.The soil columns were first slowly saturated from the bottom with distilled water in order to remove entrapped air.This saturation was performed by positioning the bottom of the soil columns at the almost same height as the distilled water level.To ensure complete saturation, the distilled water was allowed to overflow 1 cm on the soil surface during two uninterrupted days.After column saturation, the flow was reversed, and the test was performed under a constant hydraulic gradient of 1.5.This hydraulic gradient can be considered high for real conditions but was chosen to reduce the test time.
The steady state condition was achieved when no flux variation was observed with time.Subsequently, the following water flow parameters were obtained from each column: saturated hydraulic conductivity, Ksat; Darcy velocity, q ; flux, Q and average linear velocity, v (Freeze and Cherry, 1979).Total porosity, n, that represents the total of voids (connected and non-connected pores), was measured by considering the ratio of the pore volume to the total volume of each soil column.On the other hand, connected pores are represented by the effective porosity, ne, which is used in the calculation of v.The percentage of pores larger than 10 µm was adopted based on the testing by mercury intrusion porosimetry (Figure 1) (Ahuja et al., 1984).
The pH and electrical conductivity, both of the input solution and of the percolated effluent solution, were monitored.Percolation occurred in a downward flow.Three columns were percolated with NaCl and the effluents were collected at each pore volume and stored in refrigerated plastic containers.The concentration of the initial solution inleted to columns 1, 2 and 3 were, respectively, 0.050, 0.090 and 0.177 kgm -3 .At the end of the tests, the columns were percolated using three pore volumes of distilled water, with a pH ranging from 6.8 to 7.2.

Chemical analysis of the solutions
Na + concentrations were determined using a MICRONAL B26 flame photometer after appropriate dilution.Based on the relative concentration data (C/Co) versus the respective pore volume percolated with the effluent, the BTC for the Na + were obtained.

Determination of hydrodynamic dispersion coefficient
The Dh was determined using equilibrium methods proposed by: 1) Ogata and Banks (1961); 2) Brigham (1974); 3) Singh (1998); 4) Bear (1961); and 5) Delgado (2007).The non-equilibrium method applied was that proposed by van Genuchten (1981) , while the Dh was determined using the CFITIM code for the third type boundary condition (constant flux).This code is part of the Windows-based computer software package Studio of Analytical Models (STANMOD).Table 1 displays the methods and their main characteristics.

Soil characterization and flow parameters
The mineralogy consists of quartz, kaolinite, and gibbsite.Texturally, the soil consists of 79.4 sand, 10.5 clay and 10.1% silt particles (mean values).The specific weight of the solids was consistent with sandy soil values (26.4 kNm -3 ).Table 2 displays the physical characteristics of each of the undisturbed columns.The columns display physical characteristics and values consistent with the type of soil under investigation.The CEC value (2.903 cmolkg -1 ) suggests a soil with a low capacity to adsorb cations by electrostatic adsorption.According to the soil salinity classification of the Food and Agriculture Organization of the United Nations (FAO), the electrical conductivity values indicate small amounts of dissolved salts (55.70 mSm -1 ) and a nonsaline soil (Brigham, Reed, & Dew, 1961).The soil has values of 5.23 and 4.20 for the pH in H2O and in KCl, respectively.The soil presented values of 5.23 and 4.20 for pH in H2O and in KCl, respectively.The negative delta pH (-1.03) and PZC lower than the pHH2O (4.18) indicate a predominance of negative charges, which can promote cation adsorption.In addition, the soil contains a small amount of organic matter (1.18%).
The porosimetry results were discussed based on Figure 1 and on the porosimetry analysis report.According to the classification scheme proposed by Koorevaar, Menelik, and Dirksen (1983), in which the diameters of micropores, mesopores, and macropores are, respectively, <30 µm, 30-100 µm and >100 µm, the porosimetry test indicated that the pores consisted of approximately 5 macropores, 50 mesopores, and 45% micropores.Besides that, Figure 1 shows the presence of two primary pore sizes, i.e., diameters of 0.01 to 0.03 µm and 10 to 100 µm, indicating the presence of dual-porosity in the form of two frequency peaks.Model used for equilibrium of a one-dimensional transport situation with good approximation.Considers isotropic and homogeneous geological material.Can present problems in short columns, that is, when the mixed zone is about the same length as the porous medium and, so, the dimensionless dispersion, g = (v L Dh -1 ), is small.
; Bear (1961) Simple and fast method.Can present problems due to dispersivity coefficient values adopted and procedure to obtain it.
Considers equilibrium advective-dispersive transport.The Dh is commonly estimated from the BTC and the dispersivity is then obtained by the Bear (1961) equation.Brigham (1974) Based on the Ogata and Banks (1961) method, but uses a different part of the BTC.
Adequate for short columns.Considers sorption aspects.Very easy to use.

Delgado (2007)
Empirical method with a specific equation for different flow conditions.Considers a molecular diffusion coefficient and Péclet, Reynolds, and Schmidt numbers.
Valid for homogeneous geological materials in terms of porosity and particle size.
Considers the equilibrium transport model.Depends on the linear average velocity to obtain the Péclet number.Depends on the tortuosity coefficient and, consequently, on the values of the effective molecular diffusion coefficient.A visual analysis indicated that the soil also contains macropores larger than 100 µm that were not identified in the porosimetry analysis, mainly due to the small size of the soil sample analyzed in that analysis, i.e., cylindrical samples of 1 cm in diameter and 2 cm long.These pores range up to 4 mm in diameter and were very important in determining solute transport mechanisms.

Analysis of breakthrough curves
Figure 2 displays the breakthrough curves for Na + .The time of duration of the tests was set as the dimensionless number of pore volume (T = vt L -1 ), allowing comparisons between the columns.The BTC tails show that at the beginning solute speedily left the soil column, followed by a slow transport and tailing.Thus, Na + concentrations did not reach inflow concentrations in any curve, probably an indication of physical non-equilibrium (Vanderborght, Timmerman, & Feyen, 2000;Dousset et al., 2007;Jarvis, 2007).The observed behaviors can be explained by the presence of dual-porosity in the soil, in which some portions of the solutes can leave the column due to advection and others can be delayed, due to percolation through micropores and non-connected pores, causing the ions to remain in the column throughout the test.This has also been mentioned in other papers (van Genuchten & Wierenga, 1976;Jarvis, 2007;Silva et al., 2016).In addition, according to Gerritse (1996), a non-sigmoidal shape could be due to the presence of heterogeneous adsorption sites in the soil, caused by the small-scale heterogeneity.Finally, the results of the BTC supporting laboratory tests clearly show that the soil samples contain macropores that strongly influences flow and solute transport.

Hydrodynamic dispersion coefficient
Analyzing the porosity and soil characterization, as well as the BTC (Figure 2), a typical non-equilibrium transport behavior is noted, particularly at the early (up to 4 pore volumes -Figure 2) and late (from 8 pore volumes -Figure 2) behavior in the BTC tails.This behavior may be explained by the fact that that the mobile region contains large pores that fill with a solute that left the soil column quickly and appears at an earlier BTC time.Subsequently, a significant tailing occurs at the late time, due to solute diffusion into immobile regions.This is also demonstrated with the comparison of errors resulting from the fit of the Dh values for sodium obtained with the models proposed by Ogata and Banks (1961) (equilibrium) and van Genuchten (1981) (non-equilibrium).The experimental data were best fitted to the non-equilibrium method and the sum of the squares errors was up to five times lower, demonstrating the suitability of the non-equilibrium method for the studied soil.Subsequently, the Dh was calculated using all method described in the materials and methods section.Table 3 displays the parameters obtained for each method.The average results of equilibrium methods were compared with those obtained with the non-equilibrium method.The latter method was the most adequate for representing the data, due to the aforementioned analysis.The differences were analyzed based on the study of conceptual models of each method.These results are displayed in Figure 3.
Higher errors were obtained through Ogata and Banks's method (1961), probably due to the fact that the fit of the experimental data using this method is truncated at the zero-initial time, thus assuming that the relative concentration at this time is also zero, or very low and that it increases very slowly.Moreover, the BTC adopted in this method has a sigmoid shape and the inflection point occurs at a relative concentration equal to 0.5, a behavior that was not observed.In addition, this method was developed exclusively for homogeneous soils, which is not the case herein.
Acta Scientiarum.Technology, v. 41, e35419, 2019   Before using the Brigham (1974) method, the dimensionless dispersion was determined for all soil columns.Since the values were 8.9, 6.7 and 9.8 for columns 1, 2 and 3 respectively, this method was deemed suitable.However, the probability plot of the function U did not result in a straight line, demonstrating that the boundary conditions used in this method would not be suitable for representing sodium Dh in the studied soil.The probability plots of the U-functions resulted in concave upward lines at high concentrations, a behavior more consistent with the presence of non-connected pores (Brigham, 1974).It is important to mention that this method is based on the Ogata and Banks equation and showed low error values in relation to the van Genuchten model (1981).As this method is not suitable for the studied soil, the low error can be obtained because the fit of the experimental data using this method has a greater degree of

Sodium
freedom compared to the Ogata and Banks method, and it allows the curve to be shifted to a start time different than zero.Thus, a larger amount of data points passes the shifted curve, causing smaller errors than those resulting from the Ogata and Banks method.The method proposed by Singh (1998), although based on Ogata and Banks, also resulted in small errors, which can be explained by a different mathematical approach used in this method, allowing the consideration of higher concentrations since the beginning of the test, due to the greater degree of freedom of the fit, as occurred with the method proposed by Brigham (1974).Furthermore, with the displacement of the BTC, this method allows the consideration of early and late arrival times resulting from the percolation of the solution by macro and micropores, respectively, therefore justifying their small error values compared to van Genuchten (1981).
The method proposed by Delgado (2007) presented the lowest average errors.In this method, molecular diffusion is considered.In addition, it exhibits a great detail of the flow parameters when considering the Peclet and the Reynolds numbers.Therefore, this method can fairly represent an advective-diffusive transport without delay.However, the equations proposed by Delgado (2007) are based on results of several laboratory tests using nonreactive solutes in homogeneous porous media in terms of porosity and pore size.In addition, the low Dh values obtained should be more related to the importance given to the molecular diffusion, with values in the order of 10E-09 m 2 s -1 , than to the consideration of two pore regions, unlike the method proposed by van Genuchten (1981).
In the method proposed by Bear (1961) the errors also were low.However, the BTC is not regarded in this method.The determination of the longitudinal hydrodynamic dispersion by Bear's method is performed, essentially, by considering only that the spreading of solute is Fickian and that it is affected only by the velocity differences in the pores, mechanical dispersion, and molecular diffusion.However, the results can also strongly be affected by the longitudinal dispersivity coefficient (αl).Since obtaining αl is often based on empirical relationships, which are not based on the breakthrough curve (Fetter, 1999), and do not take into account the actual behavior of the solute in the soil, this method may lead to inadequate Dh estimates.

Conclusion
The interpretation of the results must consider dual-porosity transport models to obtain valid results.Other methods would quite likely produce erroneous results.The analysis of the conceptual models of the methods used to determine the hydrodynamic dispersion led to observations of the limitations of each method and a careful analysis of the obtained values, which are related to the characteristics of the soil.Due to the good fits of experimental data and the soil characteristics, the best method for determining the hydrodynamic dispersion coefficient was that proposed by van Genuchten (1981).It is worth mentioning that the main limitation of this research is that it is site specific.In addition, only a limited amount of soil column was studied.Future research can also be performed by considering different boundary conditions, since this can influence Dh values.Lastly, it was concluded that there is no single perfect method for determining soil transport parameters.Modeling of solute transport should be performed with caution and the modeler should evaluate not only the final value but also the conceptual models and the limitations of each method as they relate to the soil under investigation.

Figure 3 .
Figure 3. Differences between mean Dh values obtained by the six methods applied herein.

Table 2 .
Physical characteristics of the soil column.

Table 3 .
Parameters obtained for each of the six methods applied herein.