Conjunctive Use of BTC and Batch Methods for Heavy Metal Transport

Marginal sediments can be used to combat punctual pollution by heavy metals in industrial zones. Such practice requires information on metal-concentration in the workshop discharge water. Knowledge about the reaction of the heavy metal with the sediment available in the landscape is of utmost importance. Modeling of batch experiments and breakthrough curves, BTC, supplies relevant information in this regard. We modeled the static batch-data by Freundlich isotherms for testing CdCl 2 aqueous solutions equilibria with sandy loam sediment and the dynamic column-data by two codes, CfitM and CfitIM, under saturated water flow conditions. Three Cd-concentrations (5, 20, and 40 ppm) were employed to investigate the conjunction of using two procedures for obtaining the pertinent parameters for the transport of such a heavy metal and the design of the adequate Cd-trap. The results showed the deviations of the two techniques due to differences in their theoretical concepts, mathematical formulation, and performance. The batch method showed utility in supplying first guesses for the retardation factor, R, to insert into the 4-parameter analytical code, CfitIM, applied for column BTC modeling. The iteratively-obtained-parameters of the Freundlich equation were then employed to generate the distribution coefficient, k d . The generated value was, in turn, used to get more fair guess for the retardation factor, R, to use as a fixed-value in the CfitIM code to get an in-depth insight into the BTC dynamics and to obtain the other pertinent model parameters. The BTC runs indicated that the concentration controls the distinctive adsorption and transport rate and behavior of the heavy metal in the sediment column. The most dilute solution offered the highest Cd impediment, as shown by the most significant values for the distribution coefficient, k d , and retardation factor, R. The malfunction of the sediment as a trap appeared at Cd-concentrations four to eight folds higher than the most dilute solution. However, the loamy sand trap is successful when fed with a dilute aqueous solution. A set of successive traps is to arrange in tandem lines when moderate to high concentration is to discharge from an industrial workshop. The results emphasize the utility of the mutual use of these two lab procedures for the design of adequate traps and landfills and the simulation of more complex situations in the field. The point-pollution control needs to continue running batch and BTC experiments and to carry out their corresponding modeling. non-linear best fit values of the code parameters (2-parameters, P and R the case of the CfitM code, 4-parameters, P, R, β ω the case of the CfitIM code). conjunction of the replica of batch CfitM CfitIM models of BTC postulated


Introduction
The dynamic study of pollution in the field is a complicated and costly task. It needs a lot of technical resources that are not available to many research institutions. The computer codes based on the numerical solution of solute transport problems in the porous media offer an exciting arena to consider for running useful simulations that avoid costly field-study of pollution by heavy metals. The use of such codes requires, however, a priori knowledge about the values of the model parameters corresponding to the involved sediments and solutes. The use of the breakthrough curve, BTC, technique, and batch runs in the lab may offer a smart answer to the situation. The two procedures are outstanding methods in studying the reactive and non-reactive solute transport in the porous media at low costs. The interest in applying the two techniques has recently increased as the soil, and groundwater pollution problems are escalating worldwide. However, the BTC experiments are less employed than the batch runs due to the convoluted mathematics of the column experiments. The conjunction of the two techniques may, however, be more useful than using any of them apart. This work deals with the conjunctive use of the two procedures.
The pollutant fate differs according to its type, concentration, and the material of the reception sediment. Practically, each model will have its own set of parameter-values that belong to the specified problem. Fortunately, such parameters can be estimated in a preliminary pace, from the BTC runs and batch experiments. Based on the batch data, an iteration protocol can be run to obtain a best-fit adsorption curve (e.g., the Freundlich isotherm). Special-purpose software is also often used to fit the non-linear BTC labobservations using the inverse solution of an analytical code (e.g., CFITM and CFITIM codes, igwmc, 2001) to get the required values of model parameters.
Heavy metals can move with the wastewater flowing from the outlet of an industrial workshop. The discharged wastewater may infiltrate into the soil and then percolate down to the local groundwater when the value of the distribution coefficient, k d, for the surface sediment is small. The k d parameter expresses the mobility of the solute in the porous media. In contrast, the heavy metal may be stopped and trapped in the top centimeters of the soil surface, by adsorption and cation exchange, when the value of such coefficient is high enough.
Consequently, the distribution coefficient is a crucial parameter needed for the design of adequate traps to control the propagation of heavy metal pollution via the wastewater discharge out of workshops in many industrial zones. The k d value is regulated by solute concentration in the aqueous solution, mineralogy of the solid-phase, and chemistry of the concerned solute. Metal adsorption on the solid-phase particles can take place concurrently as well during the BTC, and batch runs to different degrees as the compositional factors, and other aspects (e.g., temperature, pH.), govern the chemical reactions.
The k d parameter, as an indicator of the adsorbed solute to the remaining dissolved solute at equilibrium, offers a new linkage between the batch and BTC methods. In the mathematical formulation of the BTC, the distribution coefficient, k d, is another way of expressing the retardation factor, R, which expresses the ratio of pore-water velocity to solute velocity. In batch run data processing, by Freundlich isotherm, the capacity coefficient, K f, and the sorption-intensity exponent, N, may be related to the k d and the R parameters. Here, we call this linkage the "conjunctive use of the two techniques.". The value of the k d parameter equals the K f value only when N = 1 (as for a linear isotherm), as shown by the Hashimoto derivative.
We are aware that solute-transport in sediment and its breakthrough curve, BTC, correspond to an open-system. Solute transport in the soil column is a time-dependent, nonequilibrium, transient, and dynamic process that includes diffusion. In contrast, a batch run goes on in a diffusion-free closed-system. A batch run is a time-independent procedure in which we consider reaching a kinetic equilibrium state between the solute and the solidphase. The processing of batch run is usually carried out by the Freundlich empirical equilibrium-isotherm formula. Despite the vast conceptual gap and the sharp experiential differences between the batch isotherm and the column-BTC construct, they are considered, here, in a first approximation, as a set of two "complementary methods." Our purpose is to examine to what extent such assumed complementarity may improve the empirical and calculation load needed for the design of a successful trap. The present work deals with preliminary modeling that concentrates on the promotion of the use of marginal sediments (e.g., sandy loam materials) as a depository for Cadmium retention in a desert area where the increasing mining and industrial activities menace the local groundwater systems that are precious fossil water. Table 1. A. The relationship of Freundlich batch Isotherm and BTC, via k d and K f parameters of the batch, and the R-value to insert as a fixed value in the CfitIM input screen. B. Derivation of a relationship between k d and K f. via linearization by integration. The K f and N were independently obtained by iterative methods on batch data using the input concentrations (C_0), not the equilibrium concentrations (C _ equilibrium) after the batch experiment run A. Table 2. Example of guess value for the R parameter (Ret Fac in the table) inserted in the (CfitIM) code. Any of the four parameters (P, R, B, and ω) may be fixed or fitted. Marking the button beneath the inserted parameter value will induce the code to calculate a fitted value for the parameter. Otherwise, the inserted values will be used as a fixed value for that parameter during code execution. We have deliberately used fixed values for the R parameter as obtained from batch runs and linearization by integration Table 3. Results of batch BTC runs with CfitIM fitness to BTC. The R values were obtained from the batch runs, linearization via integration (INTL), and average slope linearization (AVSL). Equilibrium k d (= K f N C N-1 ) values were calculated for each concentration, using N=0.42 and K f = 62.35 (calculated by iteration on Freundlich formula), and then the R-value was obtained (=1+ρ k d /θ). Fixed values of R, lower than those obtained by INTL, were then inserted into CfitIM, and the code was executed several times to finally get BTC best fit and reasonable estimates for P, B, and ω. The k d values are in units of 10 -3 m 3 solution/kg soil (same as the dimensionless value), whereas the P, R, β and ω values are dimensionless Theoretical background Solute transport in porous media is an overwhelming hydrodynamic phenomenon. Regular reports on heavy metals nearby the industrial zones may give some "snap-shots" of their concentration at given depths and locations in polluted soils by the study time. In the environment, however, solute transport and fate is a continuous dynamic transient process in space and time. The complex field configurations (Fig 1) make a full study of the phenomenon a complex task that not all research workers may afford. Despite the inherent difficulties, studies using proper dynamic and kinetic approaches are available even to the modest laboratories and may supply valuable information on soil pollution problems at low technical requirements. Among these approaches, the batch and BTC methods stand distinguished to address the pollution problem in question. The conjunctive use of the two procedures is especially helpful where the mobility of the heavy metal to local groundwater via soil is a severe hazard (Wang, 2008). Many factors control the destiny of solutes in soils, including mineralogy, texture and structure, type of solute, ion concentration in water, and the natural hydrogeochemical processes taking place in the porous system. These actions include dispersion, convection in macro-pores, diffusion in micro-pores, adsorption, cation exchange, dissolution, precipitation, redox potential, and biological processes van Genuchten. 2008 andŠimůnek et al., 2009).
While it is almost impossible to carry out a full simulation of all aspects of the complicated solute-transport phenomenon in the environmental systems, the analytical and numerical models provide excellent tools for understanding solute behavior and fate in the porous systems in the lab through the use of soil columns and sand tanks (Siyal, 2010 andMacintyre et al., 1991).
Retention of heavy metal with low mobility at the soil surface (or in a sandbox trap) would prevent the pollution of the local aquifer, while heavy metals with high mobility represent a severe risk to groundwater (Fig. 2). The distribution coefficient, k d, operationally expresses solute mobility. The k d value is to obtain via the conjunctive use of batch runs and breakthrough experiments in the laboratory, as shown in the present work.
Kinetic studies carried out via batch equilibrium have long ago suggested that ion exchange is an instantaneous process (Hissink, 1924, Borland and Reitemeier, 1950, cited in Nkedi-Kizza, 1979. Batch data are usually processed using the empirical Freundlich isotherm in its raw form or after modification (Coles and Yong, 2006) to estimate the retardation factor, R, which is of primary interest to follow up heavy metal transport in the porous media. Nkedi-Kizza, 1979 in discussing the ion exchange as an instantaneous process, has written the following statement that we cite after being modified for clarity purposes, "This is expected since ion exchange is a stochiometric process with a change of enthalpy of about 2 kcal mol -1 , resembling dipole-dipole interactions (Helfferich, 1962). Diffusion limits the instantaneous equilibrium, and thus it determines the exchange rate. However, diffusion can be eliminated as a rate-limiting step in the batch studies because of vigorous shaking, but in BTC runs, particularly in aggregated soils, diffusion cannot be eliminated since the low ion diffusion rate can damp ion exchange rate to the exchange sites." Several research workers (e.g., Smith, 1968;Nkedi-Kizza et al., 1982), showed that even if ion exchange was an instantaneous process, the global kinetics could be influenced as well by fluid velocity and aggregate-size (that determine the diffusion path length) in the BTC runs.
The BTC results are to model using a code that deals with numerical or analytical equations that describe solute transport in the porous media. The computer code will fit a curve to the experimental column observations obtained in the laboratory. The generated curve belongs to the mathematical relationship between a series of dimensionless concentrations, C/C0, and the corresponding pore-volumes, PV, of the effluent draining out of a soil column. The observed BTC, which is very often a highly non-linear curve, is to optimize by the code, via an inverse solution protocol, to estimate the unknown values of the parameters involved in a conceptual solute-transport model enhanced by a non-linear fitness algorithm (Šimůnek et al., 2008, andCFITM andCFITIM codes, igwmc, 2001). The primary interest of such fitness is its operational relationship to later field applications (Porro and Wierenga, 1993;and Wang et al., 2008). The BTC technique adopted here particularly tries to extend the operating procedure to test the use of marginal sediment as trap under saturated water-flow conditions to restrain punctual pollution at the industrial aqueous outlets.
In a straightforward model of solute transport, there are usually two parameters: Péclet number, P, and retardation factor, R, in the used Convection-Dispersion Equation, CDE (van Genuchten, 1981). The retardation factor, R, is a ratio of pore-water velocity to solute velocity (at one pore volume for both). The R-parameter has no scale-dependent value. Larger than unity R-value is obtained when solute velocity is smaller than water velocity (due to cation reaction with the solid phase). When working with anion tracers, R could be small than unity. Such a small value will show anion exclusion; otherwise, it may show the presence of a stagnant moisture zone. The R formula (Eq 1) includes the distribution coefficient, k d, that accounts for solute interaction with the solid phases.
In the elaborate models, such as the Mobile-Immobile, MIM, models (van Genuchten and Wierenga, 1986), there are two more parameters, β and ω ( Fig. 3 and 4). The β parameter stands for the instantaneous retardation (governed by many factors such as the mass fraction of the solid-phase surface available to adsorption, the total and the free moisture content, retardation factor, and bulk density). The ω parameter is the rate constant (mass-transfer coefficient) of solute exchange between the mobile and immobile moisture (i.e., a product of the rate coefficient α, hr -1 , and water residence time, L/q, hr). We do not use the parameter "Pulse" in the two models used in this study since the solute was applied continuously rather than in a brief contaminant beat. The question of whether cation retention in soil column is due merely to adsorption, or depends only on cation exchange, or attributed simultaneously to the two processes, seems perplexed. If adsorption solely governs cation retention, what would be the instantaneous electric balance in the aqueous solution during a BTC run?
The so-called "salt adsorption" could be an answer to that legitimate question in some cases. For example, Qafoku and Sumner (2002) presented evidence that simultaneous adsorption (in equivalent amounts of cations and anions of an electrolyte without a net release of other ions into the soil solution, i.e., "salt adsorption" of an "indifferent electrolyte,") occurs in some acid tropical and subtropical soils in BTC runs conducted using dilute inflow solutions.
However, the same authors had previously (Qafoku and Sumner (2001) shown that the "salt adsorption" had caused a simultaneous depletion of Ca 2+ and NO3from the input solution, with no other cations or anions released in the leachate in the first 2.5 PV, and no inner-sphere complexes of cations were formed when "salt adsorption" occurred.
Later, Qafoku and Sumner (2002) assumed that overlapping of the diffusion double layers around the positively-charged Fe and Al oxides, and the negatively-charged silicate minerals, is related to "salt adsorption" of indifferent electrolytes as caused by simultaneous adsorption of the corresponding ions, in the oppositely-charged diffuse layers. Such layers get compressed in response to ionic strength increase in the next soil solution. This compression suggests that "salt adsorption" occurs when two oppositely charged solid phases are present. These authors concluded that Cs and Cl were most probably "simultaneously adsorbed" in the outer spheres (diffuse layers) on oppositely charged soil particles, while the adsorption of such ions in the inner Stern layer was not promoted under their experimental conditions. The "indifferent-ions" are those who are only absorbed through Coulomb forces by surfaces of opposite charge, not absorbed on an uncharged surface, and unable to reverse the sign of the surface charge of the solid-phase particle, but only to compress the electrical double layer. In contrast, the "specifically-absorbed" ions possess a "chemical affinity" for the surface of the solid phase particles, in addition to the Coulomb interaction. Such affinity is related to three forces, a) van der Waals, or hydrophobic bonding, b) pi-electron exchange, and c) complex formation (http://what-when-how.com/nanoscience-and-nanotechnology/mineralnanoparticles-electrokinetics-part-2-nanotechnology/). Source http://www.epa.gov/radiation/docs/kdreport/vol1/402-r-99-004a.pdf Pollution of vulnerable soils and precious groundwater resources may go with the mining activities and industrial installations. Contamination of the local aquifers can take place in the long run by heavy metals through fluid leakage. The inflow of aqueous solutions to special design traps packed with marginal sediments for temporary deposition at the pointsources can offer one of the answers to such an acute problem. The filling material used in the present work was collected as single-grain samples of surface sandy loam sediments. The sample collection site is in the middle of a 2000 km 2 oasis, nearby Al-Kharj city, KSA. The location coordinates are 24.148°N, 47.305°E (i.e., 24°8′54″N, 47°18′18″E) in the oasis found at some 50 to 77 km to the southeast of Riyadh city, Fig 1. This oasis is nearly lying at the origin of the central sandy plateau, to the middle and east of the country, with a height of about 500 m AMSL. The local, but extended, fossil groundwater system is one of the significant water resources of the country that is extensively used in the region's historically cultivated areas. Heavy pumping for domestic uses, and large-scale agricultural expansion under pivot irrigation, has dramatically diminished the hydraulic head since 1980, and severe drawdown has taken place during the last three decades beneath the quartz eolian sands and the quartz and carbonate quaternary wadis of the Al Kharj region. The substantial water withdrawal practice has rationally given rise to the disappearance of the local upward flowing springs, and their associated lakes that once existed. The Al Kharj region is part of the significant sedimentary basin of the country. The deposits of that lose sedimentary region have Calthiorthids-Torripsamments deep sandy soils, with 0-3% slopes that become nearly level on smooth plains and in the old alluvium of the closed wadi bottoms characterized by very low soil salinity and moderate to high permeability. However, small convoluted and banded surface pattern areas have Calciorthides and Gypsiorthids soils with Calcic and Gypsic horizons (Saeed, 1995). Unpleasantly, the Riyadh and Al-Kharj agglomerations became significantly menaced by industrialization-associated risks via the industrial complexes recently-installed in the outskirts of these two major cities, with the heavy-metals being the major pollutants (Al-Hamad et al., 2012;and Al-Shayeb, 2002). The most menaced site is lying nearby the coordinates 24°03′48″N 47°34′50″E in Al-Kharj area, which is home to about 10% of the local population.
Four sandy loam samples were collected just outside the southeastern wadi bottom, next to Al Kharj city (Fig. 1). The sediments have bulk densities around 1300 kg/m 3 and are made up of quartz particles. However, the solid phase carbonate contribution is in the range from 9 to 35%. The four samples have low salinity (from 0.20 to 1.40 dS/m, in 10:1 extracts), shallow organic matter content (from 0.40 to 1.20%), and moderate saturated hydraulic conductivity (from to 0.30 to 0.50 cm/min). The present work is concerned only with the sample collected from location #1 that has a sandy loam texture with 8.70% Calcite, 0.40% organic matter, and a bulk density of 1298.7 kg/m 3 . This particular soil sample was selected for the present work since it has the lowest solid-phase carbonate content, which may help to avoid potential effects of the high carbonate content on the chemical reaction of the Cadmium aqueous solution with the solid phase during the batch and BTC runs. The pF-curves of undisturbed samples from the four locations were obtained using the conventional pressurecooker technique. Then best representative pF smooth-curves were fitted to the obtained θ(h) observations using RETC code (a special-purpose software for plotting soil moisture retention curves). This moisture retention code is equipped with an advanced curve-fitting optimization algorithm. The corresponding pF curves are not shown here since this issue will be dealt with in another work on water flow in the sampled sediments with the other three samples dealt with using other heavy metals. The saturated volumetric moisture content of the sediment used for the present work is 0.5089. The pH in the 1:1 water to soil mass ratio mixtures for the collected samples was ranging between 7.4 and 8.4.
The experimental work started with the closed-system batch adsorption experiments (Amacher et al., 1988) using Cadmium Chloride aqueous solutions on three sub-samples from location #1 in three runs, each in duplicate, using three concentrations. The three Cdsolutions used for the batch and the BTC experiments had nominal concentrations of 5, 20, and 40 ppm Cd, but the real concentrations were slightly different, as shown in Table 1. The soil material and CdCl2 solution were mixed for the batch run in a 10:1 ratio, as often used (Qualls, R, andHaines B., 1992 andVenkata et, 2007), and exposed to vigorous electric shaking for 12 hours. Then the solid phase particles were forced to settle-down by centrifugation (for 20 minutes, at 6000 rpm) to get a clear supernatant for chemical analysis of the residual soluble Cadmium, assuming that equilibrium has taken place (Amacher, 1991 andDaniel et al., 2007). Atomic absorption apparatus was used for the measurement of Cd concentrations (Unicam 969 AA Spectrometer, Page 1982). Batch data-points were used in Excel ® for non-linear iterative fitness of Freundlich, Langmuir, and Freundlich-Langmuir formulas using the built-in SOLVER functions in the workspace as given in the link (http://www.ars.usda.gov/services/software/download.htm?software=true&modecode=64-45-00-00). Only the Freundlich equation was considered in the rest of this work.
Breakthrough curve, BTC, runs under open-system saturated flow conditions, were posterior to the batch runs, and conducted using a locally made inflow-outflow setup (Fig. 5) where the PVC column was 50 cm long, and 5 cm in diameter. The sediment was carefully packed to an average fraction of bulk-porosity, ρ, of about 0.50, and its one-pore-volume, PV, was about 500 ml. Effluent fractions were successively collected for 12 hours for each BTC run and stored at 4°C for later chemical analysis. The observed BTC data-points of Cdconcentration in the effluent corresponding to the observed pore-volumes were used in the STANMOD software package (CFITM and CFITIM codes, igwmc, 2001) to conduct the non-linear best fit of the observations to obtain the unknown values of the code parameters (2-parameters, namely P and R in the case of the CfitM code, but 4-parameters, namely P, R, β and ω in the case of the CfitIM code). The conjunction of the Freundlich replica of batch results and the CfitM and CfitIM models of the BTC runs were postulated for this present study with the design of adequate sand traps for industrial pollution control in workshops of the Al Kharj city in mind depending on the values of the obtained parameters.

Results and Discussion
This exploratory work is concerned with Cadmium adsorption on local coarse-grained sediment (sandy loam) for the estimation of the transport parameters of such a heavy metal as a preliminary stage for the design of adequate traps for pollution control using regional marginal sediments and subsequent simulations relative to a comprehensive fieldwork campaign. This study is the first step in a broader work on the adsorption of five heavy metals (Cd, Ni, Cu, Zn, and Pb) on four local sediment textures (ranging from sandy to sandy loam) with different contents of calcium carbonate (in the range for 9 to 35%).
The BTC runs were conducted on sandy loam sediment collected from the Al-Kharj area, Fig 1, KSA, using three arbitrary spiked Cd concentrations of 5, 20, and 40 ppm in aqueous solutions, (Fig. 8 and Table 3) after conducting the batch runs on the same materials. Problems with the constant saturated flow in the columns were met with the used downward flow regime , Fig 4 A, and repeated work (reported elsewhere) was carried out until the most stable water discharge rate was obtained. Given this situation, it will be much better to use upward flow , Fig 4 B, in future work to avoid air bubbling. The analytical data of the BTC runs is not enough in itself to get the required solute transport parameters. The solute transport parameters are to regularly derive from the non-linear best-fit of BTC data using adequate computer code (Table 2). Batch runs were to perform to use the obtained Freundlich constant K f to get the first guess for the R parameter via the k d coefficient (Table 1 A and B, Figs. 6 A and B, and Fig. 7 B). The last figure, Fig 7 B, is a superposition in Fig. 7 A (copied from Selim et al., 1990) for comparison purposes. The conjunction of the batch and BTC results was necessary to render the work on the CfitIM computer-code much less timeconsuming. The seeding value, on an excellent R input, as an affixed value, in the screens of the CfitIM code, was obtained from the k d value via the K f value derived from the iteration on batch data in Excel. However, the fixed R-value was deliberately selected to be slightly lower than that obtained via the processing of the batch data for each Cadmium solution. Without the empirical conjunction of the batch and BTC experimental procedures, and without using fixed R-value, the robust CfitIM code fails to reach any best-fit solution. According to our experience, the correct conjunction of the analytical results of the batch experiments and the BTC data fitness in the CfitIM code (STANMODE software package, Šimůnek et al., 2008) has significantly reduced the time needed to carry out the modeling efforts shown in this study.
Despite the familiarity with the CfitM code for the fitness of the BTC data, the CfitM code did not fit the Cd-transport results obtained using the device shown in Fig. 5. Consequently, we shifted to use the CfitIM code and obtained the three curves shown in (Fig.  8 A). For the CfitIM code, initial guesses for four parameters (the retardation factor, R, column Péclet number, P, β and ω) are needed, while guesses for only two parameters (R and P) must execute the CfitM code. The relevant concepts of the two codes are given in Table 1 A and B.
We used the capacity coefficient, K f, of Freundlich isotherm to get the value of the distribution coefficient, k d. This coefficient was then used to generate the value of the retardation factor, R. We insert such a value, in turn, as a fixed guess in the CfitIM input screen. The other three parameters (P, B, and ω) were introduced each as a free guess. The conjunction of the batch and BTC methods is explained below.
The problem with the CfitIM code is that the user must introduce four guess-parameters to be perfected by the software, with one or more of these parameters to insert as a fixed value. We have reduced the task to get the first guesses by obtaining the value of a given parameter from the data of another experimental work. We obtained the column Péclet number, P, from another experimental work since such a P-argument is only a soil material-dependent, not a chemically dependent parameter when diffusion is negligible. We processed the batch data with Freundlich isotherm iterative fit in an Excel workspace to get the value of k d, and then the value of R, via K f and N of the Freundlich formula. Any linear isotherm for the batch data will give unrealistic k d values. Besides, in Freundlich isotherm processing, we have used the C0_BTC (the initial concentration of the aqueous inflow solution), not the Cequil_batch (the final concentration of the solution after the batch run) since the Cequil_batch values proved to be not adequate for the linearization needed for the CfitIM code as the Cequil_batch concentration has produced very high and unrealistic k d and R values.
Moreover, for each Cd-solution used in the BTC runs, we tested a band of R values (that are all slightly less than the R-value obtained via the linearization by the integration of the batch results). Among these last R-values, we picked only one R-value to insert, as a fixed R-value, in the CfitIM screen, and then the BTC best-fit modeling by the CfitIM code was finally executed. This conjunction procedure, nonetheless, required that we execute the CfitIM code as many times as needed to get the best-fit of every set of BTC data.
A separate table was updated to follow the try and error execution of the CfitIM code by showing the successive guesses for the R parameter. Such an extensive table was not to show up here. Instead, we mention, in the next sentence, an example of its role in such a follow-up. With the free guess of R = 40, the R-value obtained from the non-equilibrium CfitIM code was 44.14, when the first-type boundary condition was in use. Another guess for R-values was provided from the Freundlich isotherm as R = 27.57, under the third-type boundary condition, and the R-value obtained from the CfitIM code was 38.1. This brief example shows that the CfitIM code is highly sensitive for the inserted R-value guess. This code is also sensitive for the inserted guesses for the parameters β and ω but there is no procedure known to get initial guesses for those two specific parameters.
In contrast, the conjunction of the batch and BTC runs has supplied an exciting way to get the best R-value to be introduced in the CfitIM screen as a fixed value (that is not changed during the code execution). However, for the column Péclet number, P, an independent estimation may also be provided via an independent BTC experiment, but this is possible only when diffusion in the sediment is negligible. Negligible diffusion makes the Pargument depend only on the type of soil material, not the used chemical element or solution concentration.
The conjunction of the two experimental procedures implies that we consider the Rvalue, provided by the batch run, as an upper ceiling (a slightly exaggerated R-value) for the BTC run, since k d value is always exaggerated in the batch mangling. Accordingly, the Rvalue to be fixed in CfitIM should be lower than the R-value provided from the linearization of Freundlich isotherm by integration. Such conjunction conforms with the extreme nonconcordance known between the batch method and field conditions. However, the consonance will gain reproach to the field conditions in the case of the BTC technique (to what extent? nobody can exactly say for the time being). An upper limit for the R-value is to use in this context as a simple guideline for selecting a band of R-values to be inserted in the CfitIM screen. With each new R-value in this band introduced to the code, the software execution must be anew run, and the graphical and textual tabulation of the code is to be thoroughly reviewed to see how the best fit of the BTC is approximated. The statistical parameters of the code are to use to judge the goodness of the BTC curve fit precisely. This procedure is better than using free R-value guesses, which make using the CfitIM code a time-consuming task that may lead to no BTC curve best fit. However, an exaggeration of the k d values, and then the R-values, in the batch run, is inherent in the batch method due to the basic concept and the practice of its mangling. Such exaggeration may, however, decrease, or increase, via the linearization by integration, INTL, of the Freundlich isotherm (Fig. 6 A).
During the BTC's lab work, we also observed minor irregularities. This situation later imposed removing four high-end data-points, for the 40.48 ppm Cd-solution, before fitting the curves, Fig. 8, as these data-points were outshooting. Some irregularities showed up in the early runs. Anomaly in the C/C0-Pore volume relationship is due to some events that may take place during the BTC run. Such inconsistency is related to the change, over-time, in some of the following chemical and physical processes: 1) pH change, 2) ionic strength shift by later desorption, 3) change in the redox potential by micro-organisms, 4) formation of some complexes and thin iron-oxide coatings, called cutans, on the quartz and solid-phase carbonates, and a tiny precipitation of Cd-carbonate or hydroxycarbonate, 5) change in the saturated regime that resulted in non-steady-state flow at some column depth, 6) change in the temperature (that resulted in the difference in the flow rate and giving rise to some desorption), and 7) tardily column close-packing deterioration (Huang et al., 1993) or migration of fine solid-phase particles with water flow downward in the column. Despite this long list of theoretical reasons for such a deviation, we may refer to our own experience that says that the most crucial reason for the imperfection is primarily the gradual slight shift from the saturated steady-state water flow conditions. Accordingly, carrying out repeated runs on the same soil and chemical setup is an absolute necessity. The recurrent measurements are to follow in the future using an upward flow regime (Fig. 5 B). Moreover, we must run such a lab work using as many Cdconcentrations as we can on the same sediment. Also, the improvement in the lab work may include better packing to ensure the dominance of the saturated steady-state flow conditions before the application of the heavy metal solution. The present BTC runs were carried out using a downward water-flow regime (Fig 5 A). However, the upward water-flow regime (Fig 5 B) would be better to follow, in the coming BTC experiments, to prevent the entrapment of air bubbles, to remove most of the ambiguity on the C/C0-PV data-points, and to shrink the number of the needed replicates.
The most visible results of our conjugated experiments are that it showed the immense impact of the pollutant concentration on the position and shape of the BTC's curve to fit in the C/C0-PV plot, and the relationship of the values obtained for the involved parameters with the concentration of the used heavy metal (Fig. 8 A and B).
The different positions and shapes of the three BTC curves for the three different Cdconcentrations show how the change of the metal adsorption behavior takes place on the soil material used with the change of the applied concentration. The most dilute Cd-solution showed the highest impediment of the heavy metal, as shown by the highest kd_INTL value (10.42*10 -3 m 3 solution/kg soil, i.e., 10.42 as a dimensionless value) for the 4.96 ppm solution, Table 3, at the intersection of raw 10, and column 2. The kd_INTL value has decreased to 7.26 and then to 6.22 with the increase of Cadmium concentration in the input solution to 21.68 and 40.48 ppm, respectively, as seen at the intersection of rows 11 and 12 with column 2 in Table 3. This result is not astonishing, as it might look at first glance. The value of N=0.75, for the Freundlich isotherm, is a fair value for the sorption intensity, and the K f = 13.55 is a fair value for the capacity coefficient. Besides, many authors have shown that Rvalue increases with the decrease of the solute concentration in the aqueous phase (e.g., Davidson et al., 1980, page 26). The hydrodynamic parameters involved in BTC modeling do not exist in Freundlich isotherm conceptualization of the batch run. The observed increase of sorption in the most dilute Cd-solution agrees with other authors (e.g., Nkedi-Kizza 1979 who observed, for aggregated Oxisol, that the adsorption of 45 Ca increased with the decrease in Ca concentration at any given pH.) We perceived that the modeled CfitIM output of the BTC runs show lower values for the R-parameter as compared to the batch results (Table 3, in the rows 10-12 intersection with the column 3 and 4), i.e., the batch runs generally overestimate the adsorption S value. It is accepted that the BTC procedure gives a better appreciation of cation adsorption than the batch-technique. Dynamic flow in an open-system is prevailing in the BTC run while the batch experiment is based on a closed-system equilibrium. Besides, the formulas and concepts of both methods are different. However, in the present work, we were able to show how these two different techniques may be conjugated to make the use of the CfitIM code much more efficient.
The non-equilibrium solute-convey in the BTC is due to several reasons. The first disequilibrium is chemical. Such a chemical shift from equilibrium is due to kinetic adsorption. The chemical disequilibrium process is dealt with in the CfitIM code by the instantaneous adsorption and first-order kinetics. The second disequilibrium is physical due to heterogeneity in the saturated flow regime. Such heterogeneity is dealt with in the CfitIM code through dual-porosity (macro-and micropores) and solute mass-transfer rate constant ω between the mobile and immobile moisture regions and governed by a first-order process.
The value of β parameter sharply decreases at the highest Cd-concentration (40.48 ppm), Table 3, not due to any decrease in the mobile moisture fraction, θm, but solely by the small fraction, f, of the solid phase surface available to adsorption. Also, the value of the rate constant, ω, sharply decreases at the highest Cd-concentration due to a low solute transfer rate coefficient, α, between the two moisture regions, but not via any change in the water residence time, L/q, in the loose soil moisture region (L is column length, and q is Darcy velocity in the free moisture region.) When applying the most dilute Cd-solution (4.96 ppm) to the soil column, the water in the free moisture zone, θm, Eq 21, in the Supplementary File, will not increase. However, there will be a significant fraction, f, of the solid-phase surface (a fraction that is close to unity) available for adsorption (accompanied by a low sum for the adsorbed solute, S, high distribution coefficient value k d, and high retardation value, R), Table 3.
If these values for such parameters were pertinent to the surface-layer in a soil profile in the field, they would show less solute mobility to the underneath unconfined groundwater reservoir, and only a tiny amount of the heavy metal pollutant would move down and reach the water table, Fig. 2.
Accordingly, a dilute Cd-solution leads to a high value for the β -parameter, Table 3.
As well, there is a high rate of solute transfer, α, Eq 17, in the Supplementary File, between the two moisture regions, but not high water residence time, L/q, in the free moisture zone, and this situation leads to an immense value for ω-parameter in the case of applying a dilute Cd-solution, Eq 16, in the Supplementary File.
The behavior of ω and β parameters (Eq 16 and 21, in the Supplementary File) will sharply get reversed (too low ω and β values) when applying a Cd-solution with a high concentration.
In brief, the values of the non-equilibrium parameters β and ω depend on the values of the fraction, f, of the surface available for adsorption, and the kinetic rate coefficient α, for solute transfer between the two moisture zones. The relevant formulas for β and ω parameters appear in the Supplementary File.
No sensitivity analysis for the present-day Cd-adsorption study is given here as we have only used three spiked concentrations and one sediment for our current purposes. However, we introduce an account of the propagated (combined) uncertainty, σ, that goes with the retardation factor, R, as a function of the distribution coefficient, R(k d), with ρ/θ ratio considered constant, and the column Péclet number as a function of longitudinal dispersivity, P(λ), using the root of the sum of the squares method of the GUM approach (Bell, 2001 andStant et al.., 2016). Uncertainty in the three parameters, R, k d and P In the retardation factor formula, Eq 1, in the Supplementary File, we assume a constant value for the bulk density, ρ, and volumetric moisture content, θ, under the applied saturated flow conditions, while the distribution coefficient, k d , will change with the Cd concentration, The uncertainty in R is i.e., R = 18.60 ± 0.26, and the relative certainty in R = (0.2552/18.60)*100 = 1.37%. As the retardation factor is an indirect expression of the distribution coefficient, the percent relative uncertainty in k d would have a value close to that of the percent relative uncertainty, 1.37%, in R, for the Cd concentration of 4.86 ppm. Accordingly, the value of the σ in the k d = 0.0950 since the value of the k d is 6.90. The approximation for the k d, obtained by Excel's What If iterations, is 0.0947, as shown in Table 4 A. The arbitrarily assigned value of 0.1000 for σ in k d is shown in Table 4 B, instead.  Table 3). B. We assigned the value 0.1000 to σk d, and σ% in k d becomes 1.45%. Worth mentioning, σ% in R visibly increases downward (corresponding to the Cd-concentration increase in the first column of Table 3 at  The function P (λ) depends only on the value of the longitudinal dispersivity, λ., since the column length, L, is a constant value. The derivative of P concerning λ is dP dλ = -L λ -2 . The uncertainty in P is γιϖεν βψ σ If P = 4.42, L = 50±0.10 cm, λ = L/P = 50/4.42, λ = 11.31 ± 0.01 cm, the uncertainty in P is
Using short columns will slightly decrease the uncertainty in P; however, there is a trade-off as using a short column will decrease the proper P-value. In contrast, using a column with adequate length (30 to 50 cm) will favor the minimization of the deviations of the values of the effective longitudinal hydrodynamic-dispersion coefficient, DL, and the distribution coefficient, k d (van Genuchten, 1981, Report No 118.) The low uncertainties going with the retardation factor, R, distribution coefficient, k d, and column Péclet number, P, show the excellent quality for the obtained analytical results, albeit the potential presence of minor secondary sources of uncertainty.   (Selim et al., 1990, in their figure #2-1, page 11) reference diagram for different soils. B. Superimposition of our log-log diagram (Fig 6 B) in the reference diagram (Fig 7 A). Our data-points (represented, in Fig. 7 B, by the oblique bold blue line) have a slightly lower slope than most lines in the reference diagram (Fig 7 A) and occupy a low position. This observation is due to the coarse texture of our sediment. The three vertical bold grey bars, to the right-hand side, in Fig 7 B, show the three concentrations for the Cd solutions used in the present work.

Conclusions
We have tested the conjunctive use of the batch and BTC methods to obtain the retardation factor, R, to use in breakthrough curve (BTC) modeling, via the distribution coefficient, k d, derived from the Freundlich constant K f of the batch before using the CfitIM code in modeling the transport of Cadmium in sandy loam soil material. The approach is superior to the free speculation on the values of the four-parameters to use in the CfitIM code (namely, R, P, β and ω). In contrast, repeatedly introducing incorrect guesses is tedious and may finally make the use of a robust code a very time-consuming task.
We have seen that the CfitIM code is more involved than the CfitM code for simulating the transport of a chemically-reactive solute in the sediment since the first code accounts for non-equilibrium conditions and solute partitioning into the mobile and immobile moisture zones, and addresses convection transport as being dominant in the mobile region. Moreover, the CfitIM code considers diffusion to the available adsorption sites on the solid phase as the dominant process in the immobile moisture zone. We propose the conjunction of the BTC technique (modeled by CfitIM) and the batch method (modeled by the Freundlich isotherm and the linearization by integration) to save time in studying the contamination by reactive solutes. However, the batch's R-value should be used as an upper limit for the R-value to fix in the CfitIM. The user must be sure of the dominance of the saturated steady-state flow conditions during the BTC run to remove ambiguities on the data and to improve the values to assign to the involved solute transport parameters.
The R-parameter for Cd-transport in the used sandy loam soil column showed lower values than those obtained from the Freundlich isotherm for the same sediment and solute concentration in the batch run. The significant decrease of R-parameter in the BTC may show that it is the batch method that overestimates such a parameter due to the conceptual basis, formulation, and high water to soil ratio in the batch. The values of the parameters k d and R of the BTC have increased with the decrease in Cd-solution concentration, showing more impediment of Cd from the more dilute aqueous solution.
The dilute Cd-solution has resulted in a high fraction, f, of the surface of solid-phase mass accessible to adsorption (i.e., the BTC showed a high value for the instantaneous retardation parameter, β). As well, the dilute solution resulted in a high first-order kinetic rate α, for solute transfer between the two moisture zones, mobile and immobile, (i.e., the BTC showed a high value for the mass-transfer coefficient parameter, ω, of solute exchange between the two moisture zones). Consequently, for the practical application of sand traps in the control of Cd at an industrial outlet, it is recommended to discharge the dilute wastewater into a single repository sediment trap. However, in the case of discharging aqueous liquid with high metal concentrations, it is recommended to use a series of traps in a line. Details of trap design and regeneration must be developed in function of the values of the involved transport parameters and the size of the workshop that discharges heavy metals to the environment.
The observed low uncertainty associated with the two parameters, k d, and P, show the excellence of the reported analytical data.