Glyphosate Transport in Two Louisiana Agricultural Soils: Miscible Displacement Studies and Numerical Modeling

Abstract

Glyphosate (N-(phosphonomethyl) glycine) (GPS) is currently the most commonly used herbicide worldwide, and is generally considered as immobile in soils. However, numerous reports of the environmental occurrence of the herbicide coupled with recent evidence of human toxicity necessitate further investigation as to the behavior of GPS in the soil environment. Batch sorption studies along with miscible displacement experiments were carried out in order to assess the mobility of GPS in two Louisiana agricultural soils; Commerce silt loam and Sharkey clay. Batch results indicated a high affinity of both soils for solvated GPS, with greater affinity observed by the Sharkey soil. GPS sorption in the Commerce soil was most likely facilitated by the presence of amorphous Fe and Al oxides, whereas the high cation exchange capacity of the Sharkey soil likely allows for GPS complexation with surface exchangeable poly-valent cations. Miscible displacement studies indicate that GPS mobility is highly limited in both soils, with 3% and 2% of the applied herbicide mass recovered in the effluent solution from the Commerce and Sharkey soils, respectively. A two-site multi-reaction transport model (MRTM) adequately described GPS breakthrough from both soils and outperformed linear modeling efforts using CXTFIT. Analysis of extracted herbicide residues suggests that the primary metabolite of GPS, aminomethylphosphonic acid (AMPA), is more mobile in both soils, although both compounds are strongly retained.

1. Introduction

Glyphosate (N-(phosphonomethyl) glycine) (GPS) is the active ingredient in several herbicidal formulations and is the most widely used herbicide in the world. Such widespread use can be attributed to the compound’s high weed killing efficiency coupled with its incorporation in agricultural systems utilizing Roundup ready cultivars. Although generally regarded as relatively immobile in soils, several investigations indicate that GPS transport from application sites to the surrounding environment is significant [1,2,3,4]. These results coupled with recent evidence of human toxicity [5,6,7,8] require further investigation into the mobility of GPS in soils.

Several field studies have been conducted in order to assess the potential for GPS leaching in agricultural settings. Kjær et al. [9] conducted an eight-month field trial to assess GPS mobility and detected the herbicide in tile-drained leachates at concentrations that exceeded European Union maximum admissible levels for up to seven days after rainfall events, which the authors attributed to macropore transport. Tiles were located one meter below the ground surface in an area with a shallow water table (1–3 m below the surface), thus it was concluded that well developed soil structure could allow strongly sorbed solutes to contaminate shallow groundwater. GPS contamination of groundwater was also reported by Torstensson et al. [10], the magnitude of which was related to application rate, although the majority of applied herbicide was retained in the upper 30 cm of the soil profile. Additionally, Strange-Hansen et al. [11] reported leachate concentrations of up to 1300 µg L⁻¹ from various types of gravel commonly used in surfacing applications in Denmark. Despite evidence of GPS mobility in field soils, several studies have indicated very limited leaching potentials. Conducting a lysimeter study over the course of 748 days, Bergström et al. [12] recovered only 0.009–0.019% of applied GPS in leachates taken at 1 m below the surface despite high precipitation amounts shortly after application, and Al-Rajab et al. [13] recovered only 0.28% of the applied amount from 25 cm below the surface over an 11 month period. It is therefore implied that leaching potential in a field situation is highly dependent upon soil type and structure, a conclusion that is in agreement with others [14,15].

In order to gain insight upon the physiochemical mechanisms that dictate GPS transport through soils, laboratory studies coupled with numerical modeling have been undertaken. To assess the effect of the pore water velocity, and therefore indirectly the effect of chemical kinetics, Beltran et al. [16] obtained GPS breakthrough curves from a sandy soil at different flow rates. The authors applied a model that incorporated a time dependent Fruendlich-type partitioning coefficient to the observed data and concluded that adsorption was limited by diffusion processes. Candela et al. [17] applied a two-site model that accounted for equilibrium and kinetic sites along with a first order sink and was able to successfully describe GPS breakthrough from two soils. It was subsequently determined that GPS sorption is a kinetic process dependent on pore-water velocities. Magga et al. [18] applied a similar model along with additional terms to estimate microbial biomass dynamics, which yielded a reasonable description of GPS breakthrough from an agricultural soil from Western Greece. Additionally, Zhou et al. [19] applied a two-site model accounting for kinetic and irreversible reactions to GPS breakthrough data, however, limitations on the extent of the observed data make it difficult to assess the performance of the model.

The maximum admissible level of GPS in drinking water is 700 µg L⁻¹ according to water quality standards dictated by the US Environmental Protection Agency (EPA) [20]. In the European Union, concentrations for any single pesticide are not to exceed 0.1 µg L⁻¹, while total pesticide concentrations are not to exceed 0.5 µg L⁻¹ [21]. The different water quality standards are attributed to varying policies. The European Union uniformly applies water quality standards to all pesticides, whereas the EPA sets limits specific to individual compounds based on toxicity study results.

The goals of this study were (1) to quantify GPS breakthrough from two Louisiana agricultural soils; (2) describe observed data using linear as well as non-linear kinetic modeling; and (3) characterize the distribution of GPS and its primary metabolite AMPA in both soils.

2. Methods

2.1. Soils

In this study, experiments were carried out on two soils of varying physiochemical properties. The first soil, Commerce silt loam, is a fine silty, mixed superactive, thermic Fluvaquentic Endoaquept formed from alluvial deposits on the flood plains of the Mississippi River in the Southeastern United States. This soil is characterized as deep, somewhat poorly drained with moderately slow permeability. Soil specifically used in this study was sampled from the surface horizon (0–5 cm) from the LSU AgCenter Sugarcane Research Station located in St. Gabriel, Louisiana (12 miles south of Baton Rouge). The second soil used in the experiments was Sharkey clay, which is classified as a very fine, smectitic, thermic Chromic Epiaquert. This soil formed from clayey alluvial deposits on low terraces of the Mississippi River and was sampled from the surface horizon (0–5 cm) near Iberville Parish, Louisiana. Soils were air dried and made to pass through a 2 mm sieve prior to experiments. Selected physiochemical properties of both soils are given in

Table 1. Select physiochemical properties of the soils used in this study.

Soil	Commerce Silt Loam	Sharkey Clay
Taxonomic Classification	Fine Silty, Mixed Superactive, Nonacidic, Thermic Fluvaquentic Endoaquept	Very-Fine, Smectitic, Thermic Chromic Epiaquert
pH	6.98	6.52
TOC a (%)	1.31	1.41
CEC b (cmol kg−1)	16.5	29.6
Sand c (%)	30	3
Silt c (%)	54	36
Clay c (%)	16	61
Ammonium Oxalate Extractable
Fe (g kg−1)	4.4 *	0.83
Al (g kg−1)	0.7 *	0.23
Citrate Bicarbonate Dithionate Extractable
Fe (g kg−1)	6.08 *	7.77
Al (g kg−1)	-	2.42

* Values obtained from Harrell and Wang [22]. ^a TOC = Total Organic Carbon; ^b CEC = Cation Exchange Capacity; ^c Grain size distribution: Sand (2.00–0.05 mm), Silt (0.05–0.002 mm), Clay (<0.002 mm).

.

2.2. Sorption

The quantification of the sorption of GPS in the two different soils was carried out following the batch technique described by Selim [23]. Three grams of soil were weighed into a 40 mL Teflon centrifuge tube and 30 mL of GPS solution at an initial concentration of either 1, 2, 5, 10, or 20 mg L⁻¹ were added to each tube. Input solutions were spiked with ¹⁴C labeled GPS (phosphonomethyl group) such that the initial radioactivity of each solution was approximately 1.67 × 10² Bq mL⁻¹. Radiolabeled GPS was purchased from American Radiolabeled Chemicals, St. Louis, MO. All solutions were prepared in 5 mM CaCl₂ in order to maintain constant ionic strength. Each initial concentration treatment was repeated in triplicate. Soil/solution mixtures were then individually vortexed for ten seconds and immediately transferred to a platform shaker and were continuously shaken for 24 h. Following this, tubes were removed from the shaker and centrifuged at 5000 rpm for ten minutes. A 1 mL aliquot of supernatant solution was then transferred to a 7 mL scintillation vial whereupon 4 mL of scintillation cocktail (Packard Ultima Gold) were added, and mixtures were subsequently vortexed until thoroughly homogenized. Samples were analyzed on a liquid scintillation counter (LSC) (Perkin Elmer TriCarb 4810 TR, Waltham, MA, USA) using five minute count times. No quench correction was made and radioactivity was recorded as counts per minute. GPS concentrations were determined by the relative radioactivity of the sample to that of the input solution. The amount of sorbed GPS was calculated as the difference between supernatant and initial concentrations.

2.3. Miscible Displacement Studies

In order to investigate the transport behavior of GPS in soils, miscible displacement experiments as described by Selim [23] were carried out. Acrylic columns (10 cm length, 6.4 cm inner diameter) were uniformly packed with air dried homogenized soil. Columns were then slowly saturated with 5 mM CaCl₂ by upward flow supplied by a variable speed piston displacement pump (Fluid Metering Inc., Syosset, NY, USA). Once a steady state flow regime was established, an approximately 10 pore volume pulse of 50 mg L⁻¹ GPS was applied to the column and was then subsequently leached with 5 mM CaCl₂ solution. Effluent was collected with an ISCO Retriever II fraction collector (Teledyne Isco Inc., Lincoln, NE, USA) and GPS concentration in solution was determined by LSC as described above. The pH of the effluent solution was monitored throughout the duration of the experiments. Immediately following the termination of experiments, columns were sectioned into 1 cm increments and residual GPS was extracted from subsamples in duplicate following a modified version of the method described by Miles and Moye [24]. Three grams of dried soil were weighed into a 40 mL Teflon centrifuge tube and 30 mL of 0.2 M KOH were then added. Tubes were shaken for 24 h and supernatants were sampled and analyzed for radioactivity by LSC as described previously. In an effort to quantify the degradation of GPS in the sorbed phase, concentrations of GPS and AMPA in the extracting solution were also determined by UPLC-MS/MS.

In order to obtain independent estimates of the hydrodynamic dispersion coefficient D, a physical parameter unique to an individual soil column, a one pore volume pulse of tritiated water (³H₂O) was applied to each column prior to the application of GPS pulses and subsequently leached with the background solution. Radioactivity of the effluent solution was determined by LSC in a manner identical to that of GPS. Breakthrough curves were modeled using the CXTFIT model [25] in the inverse mode, within the STANMOD software package to obtain estimates for D. Physical parameters for transport experiments are given in

Table 2. Physical parameters for miscible displacement experiments.

Soil	Bulk Density (g cm−3)	Water Content (cm3 cm−3)	Pore Volume (cm3)	Pore Water Velocity (cm h−1)
Commerce	1.225	0.538	167.57	0.541
Sharkey	1.330	0.496	154.74	0.512

.

2.4. UPLC-MS/MS Analysis

As both GPS and AMPA molecules have the ¹⁴C labeled phosphonomethyl group, LSC analysis alone cannot discriminate between the two compounds. To do this, samples of extracting solution were analyzed using UPLC-MS/MS with a pre-column derivatization and solid phase extraction (SPE) clean-up procedure.

Fluorescent derivatives of GPS and AMPA were obtained through the labeling of the amine group with 9-fluorenylmethyl chloroformate (FMOC). This was done by combining 1 mL aliquots of supernatant or extract solution with 0.5 mL 5% borate solution and 0.5 mL of 3% FMOC in acetonitrile. This mixture was vortexed and allowed to incubate at room temperature for 30 min. To stop the derivatization reaction, 5 mL of 1% KH₂PO₄ adjusted to a pH of 3 using 6 N HCl (pH 3 water) was added to each mixture. These were subsequently vortexed and then centrifuged for 5 min at 3500 rpm to separate any potential precipitates.

Derivatitized samples then underwent an SPE clean-up protocol. SPE columns (Oasis HLB 3cc 60 mg cartridge) were placed into a vacuum manifold set-up and conditioned with 3 mL of methanol, followed by 3 mL of deionized water, followed by 3 mL of pH 3 water. Derivatized, centrifuged samples were then loaded onto the column. After the full volume of the sample had been eluted through the cartridge, 5 mL of pH 3 water was applied, and the column was allowed to dry under full vacuum for 30 min. Derivatized GPS and AMPA were then eluted with 7 mL of methanol into 15 mL polypropylene centrifuge tubes. Tubes were then transferred to a 50 °C water bath and were brought to dryness under a stream of N₂. The resulting precipitates were then reconstituted with 1 mL of deionized water and filtered through a 0.2 µm PVDF syringe driven filter into auto-sampler vials for UPLC-MS/MS analysis.

Concentrations of GPS and AMPA were quantified using a Waters UPLC-MS/MS (TQD) instrument with Waters MassLynx software. Chromatographic separation was achieved with a Waters Acquity UPLC BEH C18 2.1 × 50 mm 1.7 µm particle column operated at 30 °C with a flow rate of 0.30 mL min⁻¹. Mobile phases of 5 mM ammonium acetate and UPLC grade acetonitrile were applied to the column using a gradient flow protocol. Sample injection volumes were 10 µL and expected retention times for GPS and AMPA were 1.06 and 1.23 min, respectively.

The mass spectrometer was operated with electrospray in the negative-ion mode. The electrospray capillary was set to 2.9 kV and the extractor cone was operated at 2.00 V. The source temperature was set to 110 °C and the desolvation temperature was set at 350 °C. The desolvation gas was nitrogen (500 L h⁻¹) and the collision gas was argon (0.20 mL min⁻¹). The multiple reaction monitoring mode of the degradation patterns m/z 390.9 → 167.9 (GPS) and 332.2 → 109.8 (AMPA) were used for quantification. For verification, the degradation patterns m/z 390.9 → 149.8 (GPS) and 332.1 → 135.9 (AMPA) were used.

2.5. Modeling

2.5.1. Sorption Model

The quantification of the behavior of solutes within reactive matrices commonly employs the use of equilibrium isotherm models in order to partition between sorbed and solution phases. One extensively used model is the Freundlich isotherm, expressed in Equation (1) [23]:

$$S=K_f{C}^n$$

(1)

where S is the sorbed concentration (mg kg⁻¹), K_f is the Freundlich distribution coefficient (Lⁿ mg¹⁻ⁿ kg⁻¹), C is the solution concentration (mg L⁻¹), and n is a dimensionless coefficient commonly less than 1.

2.5.2. Linear Modeling

Quantifying the transport and ultimate fate of contaminants in soils and groundwater provides a basis for determining their potential environmental hazard. The results of laboratory miscible displacement studies offer insight upon the physiochemical mechanisms that dictate the behavior of solutes in the subsurface environment and, thus, supply investigators with an initial idea as to the movement of chemicals in a natural or field situation. These studies are critical, as they guide management decisions to maintain the viability of soil and water resources.

The movement of non-reactive solutes in the subsurface is contingent upon three dominant mechanisms [23]. Convective transport of mass refers to transport of a dissolved chemical along with the pore water as it moves through a soil profile, the magnitude of which is a function of the pore water velocity determined by Darcy’s Law. Another important mechanism is mass flux induced by local concentration gradients, commonly referred to as molecular diffusion and governed by Fick’s Law. This is an active process regardless of whether or not there is a net water flux in the system. The final transport mechanism is dispersion, which is the mixing of the soil solution that is not attributed to molecular diffusion. In soils and geological media, this is brought about by a non-homogenous distribution of flow velocities attributed to the pore size distribution, path length fluctuation (tortuosity), and a flow velocity gradient across the width of a conducting pore (Poiseuille’s Law) [23]. Consideration of all three of these processes coupled with mass balance yields the following equation for the transport of a non-reactive solute in one dimension [23]:

$$\theta \frac{\partial C}{\partial t}=\theta D\frac{{\partial }^{2}C}{\partial x^2}-q_x\frac{\partial C}{\partial x}$$

(2)

where $\theta$ (cm³ cm⁻¹) is the volumetric water content, $C$ (mg L⁻¹) is the solution concentration, $D$ (cm² h⁻¹) is the hydrodynamic dispersion coefficient, and $q_x$ (cm³ cm⁻² h⁻¹) is the water flux density in the x direction.

In addition to the transport mechanisms discussed above, the behavior of a reactive solute in the subsurface environment is complicated by chemical interactions with the porous matrix. If it is assumed that solute mass within the solution and solid phases are at an instantaneous chemical equilibrium via a linear partitioning coefficient, the one-dimensional transport of a reactive solute is given by the linear convection-dispersion equation [26]:

$$R\frac{\partial C}{\partial t}=D\frac{{\partial }^{2}C}{\partial x^2}-v\frac{\partial C}{\partial x}-\mu C$$

(3)

$$R=1+\frac{\rho }{\theta }{K}_d$$

(4)

where R is a dimensionless retardation factor, $v$ (cm h⁻¹) is the pore water velocity ($v=q_x/\theta$), $\rho$ (g cm⁻³) is the matrix bulk density, $K_d$ (L kg⁻¹) is the linear partitioning coefficient, $\mu$ (h⁻¹) is the rate coefficient associated with the term $\mu C$ which represents a first-order sink within the system. This can be interpreted as irreversible or precipitation reactions, or mass loss from the system due to degradation or volatilization. For conservative solutes, this term is taken as zero. Analytical solutions satisfying a number of initial and boundary conditions are given by Toride et al. [27] and incorporated into the CXTFIT model within the STANMOD software package [28].

2.5.3. Multi-Reaction and Transport Model

As physiochemical heterogeneity is inherent in soil systems, it is reasonable to expect that a solvated chemical interacts with the soil matrix via a variety of different mechanisms. As such, several approaches to describe solute retention by soils involving the use of multi-reaction models have been introduced [29,30]. Presented here is the nonlinear multi-reaction model put forth by Selim and Amacher [31]. Features of this model are that mass exchanges between solution and solid phases can be conceptualized as a system of various retention mechanisms, including non-linear equilibrium and kinetic reactions, along with first-order irreversible sinks. A schematic representation of the conceptualized system is displayed in Figure 1.

Figure 1 implies a system consisting of three general types of reactive sites: (1) equilibrium sites; (2) reversible kinetic sites; and (3) irreversible kinetic sites. Equilibrium sites are always at chemical equilibrium with the solution phase via the partitioning coefficient K_E, which is analogous to the Freundlich partitioning coefficient K_f. Reversible kinetic sites are characterized by fractional order kinetics, except for first order kinetics associated with $S_3$, described by forward and backward rate coefficients. The third type of reactive site is irreversible, where sorption is considered a first-order kinetic process. These types of sites act as a sink for mass within the system, and may be interpreted as inner-sphere complexation or precipitation reactions. The multi-reaction model can be expressed mathematically as a system of differential equations where equilibrium and kinetic sites are given as [31]:

$$S_E=\frac{\theta }{\rho }{K}_E{C}^n$$

(5)

$$\frac{\partial S_1}{\partial t}=\frac{\theta }{\rho }{k}_1{C}^m-k_{2}S$$

(6)

$$\frac{\partial S_2}{\partial t}=\frac{\theta }{\rho }{k}_3{C}^b+k_6{S}_3-S_2\left(k_4+k_5\right)$$

(7)

$$\frac{\partial S_3}{\partial t}=k_5{S}_2-k_6{S}_3$$

(8)

$$\frac{\partial S_{irr}}{\partial t}=\frac{\theta }{\rho }{k}_{irr}C$$

(9)

where $S_E$ (mg kg⁻¹) is the sorbed concentration within the equilibrium sites, $S_1$, $S_2$, and $S_3$ (mg kg⁻¹) are the sorbed concentrations within kinetic sites, $S_{irr}$ (mg kg⁻¹) is the sorbed concentration within irreversible sites, $K_E$ (L kg⁻¹) is the equilibrium partitioning coefficient, $k_1$, $k_2$, $k_3$, $k_4$, $k_5$, $k_6$, and $k_{irr}$ (h⁻¹) are rate coefficients, n, m and b are reaction orders, $\theta$ is the volumetric water content (cm³ cm⁻³) and $\rho$ is the bulk density of the reactive matrix (g cm⁻³). Given this grouping of reactive sites, the total mass $S$ (mg kg⁻¹) in the sorbed phase is expressed as [31]:

$$S=S_E+S_1+S_2+S_3+S_{irr}$$

(10)

Incorporation of Equations (5)–(9) into Equation (2) yields the multi-reaction transport model (MRTM) describing reactive solute transport under steady-state water flow conditions ($\theta$ and $q_x$ are spatial and temporal constants within the system) [31]:

$$\frac{\partial C}{\partial t}+\frac{\rho }{\theta }\frac{\partial S}{\partial t}=D\frac{{\partial }^{2}C}{\partial x^2}-v\frac{\partial C}{\partial x}$$

(11)

where $S$ (mg kg⁻¹) is the solid phase concentration given by Equation (10). In order to simulate measured GPS break through curves, best-fit parameter values are obtained via a least-squares optimization procedure. Criteria for evaluating model performance are the coefficient of determination (r²) and the root mean squared error (RMSE) given by Equation (12) [31]:

$$RMSE=\sqrt{\frac{RSS}{N-P}}$$

(12)

where RSS is the residual sum of squares, N is the number of observations and P is the number of fitted parameters.

3. Results and Discussion

3.1. Sorption

The sorption isotherms of GPS in both soils are displayed in Figure 2 with optimized values of Freundlich coefficients given in

Table 3. Optimized values of Freundlich coefficients for GPS sorption onto Commerce and Sharkey soils.

Soil	Kf (Ln mg1−n kg−1) ± 95% CI	n ± 95% CI	SSE (mg2 kg−2)	r2
Commerce	158.2 ± 13.9	0.66 ± 0.12	104.955	0.995
Sharkey	395.7 ± 34.0	0.93 ± 0.08	20.516	0.999

. In general, measured data are well described by the Freundlich model (coefficient of determination values >0.99 in both cases). Optimized values for the K_f coefficient are 158.2 and 395.7 Lⁿ mg¹⁻ⁿ kg⁻¹ for the Commerce and Sharkey soil, respectively. These values are indicative of a high affinity of GPS for both matrices and are consistent with those obtained from similar studies performed on a variety of different soils [32,33,34].

The Freundlich isotherm for the Commerce soil is characterized as highly non-linear (n = 0.66), suggesting a more heterogeneous distribution of site affinities across the range of input concentrations. Because interactions between the solute and high affinity sites are, by definition, more thermodynamically favorable than interactions involving low affinity sites, high affinity sites are occupied first. As these sites are filled, interactions among less energetically favorable sites predominate. This effectively reduces the overall affinity of the reactive matrix with increased GPS loading, resulting in reduced “rates” of partitioning at higher input concentrations. The strong retentive behavior of this soil can be attributed to its relatively high content of amorphous Fe and Al oxides [35,36,37], as determined by extraction with ammonium oxalate. Additionally, the content of crystalline Fe and Al oxides may also lend to the high retentive capacity of this soil. It has been reported widely in the literature that these materials display a high affinity for solvated GPS, where sorption likely occurs via mono or bidentate inner-sphere complexation [38,39].

GPS sorption onto the Sharkey soil is nearly linear, implying a more homogeneous distribution of site affinities where the “rate” of partitioning of the solute into the solution and sorbed phases is more or less independent of the solution concentration [40]. A higher value of the Freundlich partitioning coefficient in this soil is indicative of a greater affinity of the Sharkey soil compared to the Commerce. This increased affinity relative to the Commerce soil may be attributed to several physiochemical properties of this matrix. The high cation exchange capacity of this soil may partially explain this result. While performing a regression analysis on 101 different soils, Dollinger et al. [41] found that both linear and Freundlich partitioning coefficients were highly dependent upon, and positively correlated with CEC, more so than the content of Fe and Al oxides. De Jonge and de Jonge [42] reported increased sorption with increasing ionic strength of the solution, which could potentially be attributed to GPS complexation with surface exchangeable poly-valent cations. Coupling the results of these two studies, enhanced di-valent cation (Ca²⁺) sorption brought about by the higher CEC of this soil could provide additional surface sites for the sorption of GPS, resulting in a greater affinity. Additionally, Paradelo et al. [43] performed a similar regression analysis and determined that GPS adsorption was positively correlated with the clay content of soils, which is consistent with results reported here. Since electrostatic repulsion between the negatively-charged clay surfaces and the herbicide molecule will likely prevent the formation of direct surface complexes, GPS sorption onto clay minerals within this soil likely takes place via ligand exchange with surface coordinated hydroxyl groups on the edge sites of layer silicates [44]. In accordance with the above discussion involving sorption by the Commerce soil, the relatively high amounts of crystalline Fe and Al oxides in this soil may also largely contribute to the high retentive behavior of the Sharkey soil [39]. The results of several recent studies have emphasized the importance of organic compounds in the sorption of solvated GPS, which may be of importance here due to the higher amount of organic matter in this soil (as indicated by TOC). It has been demonstrated that GPS interacts strongly with humic and fulvic acids [45], as well as with various peptides [46], which are abundant in soils. Furthermore, the lower pH of this soil compared to the Commerce likely increases the propensity for the sorption of solution phase GPS [35,47].

3.2. Transport Studies

3.2.1. Tritium Breakthrough

Tritium breakthrough curves (BTC) from both columns are displayed in Figure 3, with fitted CXTFIT parameters and fraction of recovered radioactivity given in

Table 4. Optimized values for the hydrodynamic dispersion coefficient and fraction of applied radioactivity recovered for both miscible displacement columns.

Soil	D (cm2 h−1) ± 95% CI	RMSE	r2	Fraction of Radioactivity Recovered
Commerce	1.705 ± 0.241	0.041	0.989	0.907
Sharkey	4.108 ± 0.772	0.056	0.950	0.965

. In general, the classical convective-dispersive equation (CDE) (Equation (2)) provided a good description of the measured data, with coefficient of determination values of 0.989 and 0.950 for the Commerce and Sharkey soil columns, respectively. The BTC for the Sharkey column is characterized as mostly symmetrical, indicating the presence of a single flow domain. Conversely, the Commerce BTC displays some degree of tailing during leaching, suggesting that physical non-equilibrium conditions may exist. In this case, flow is not restricted to just inter-particle porosity and, therefore, intra-particle diffusion may play a more prominent role. This result is further emphasized by the relatively low recovery of applied radioactivity (~91%) for this column, which can be attributed to limited diffusion of tritium from the intra-particle domain back to the bulk soil porosity. Recovery of applied radioactivity from the Sharkey column (~96%) is indicative of limited mass loss to the intra-particle domain and is, therefore, consistent with the overall shape of the BTC. Although evidence of physical non-equilibrium within the Commerce column exists, the magnitude is such that any major effect on GPS transport is not to be expected. Shapes of tritium BTCs are contrary to expectations, as Sharkey is well known to have highly stable aggregates. Therefore, evidence of physical non-equilibrium was anticipated in the Sharkey column and not in the Commerce column. Additional investigations are needed to further understand the dominant mechanisms of tracer transport in these soils.

3.2.2. GPS Breakthrough

GPS BTCs from both soil columns are displayed in Figure 4. Results indicate that GPS mobility in both soils is highly limited due to their high adsorptive capacities, with 3% and 2% of the applied mass recovered in the effluent solution from the Commerce and Sharkey columns, respectively. Although observed breakthrough was very low, GPS mobility in the Commerce soil was greater than that in the Sharkey soil, with peak concentrations ~4% that of the influent solution. In the Sharkey soil, the maximum effluent concentrations were approximately 2% of the influent concentrations. Qualitative differences of GPS mobility in both soils are consistent with trends determined by batch sorption studies, with more limited mobility occurring in a soil with greater GPS affinity. Very low GPS mobility in soils has been well documented in the literature. Monitoring of leachate from 1 m below a field soil over a 2 year period, Kjær et al. [14] reported complete retention of GPS within the soil column, and Napoli et al. [48] recovered only an average of about 0.82% of applied GPS from the same depth over the course of a year. Additionally, Al-Rajab et al. [13] reported less than 0.28% of applied GPS was recovered from a 25 cm soil profile, and Bergström et al. [12] recovered only 0.009–0.019% of applied GPS in leachate sampled at 90 cm depth over the course of 748 days. As field conditions are inherently variable, results reported from laboratory studies may provide a more direct comparison to results presented here. Using flow rates similar to those used in this study, Candela et al. [17] observed maximum concentrations of GPS in the effluent from a Spanish surface soil at approximately 1% that of the influent concentration, even though columns were only 2 cm long and 150 pore volumes of GPS solution were applied. Significant breakthrough was observed in this study when flow rates were increased by two orders of magnitude; however, it is unclear whether such rates are realistic in a field setting.

It is noteworthy that GPS breakthrough from both soil columns is rapid with no apparent lag time, suggesting high connectivity in a small proportion of the pores, which allows for a limited fraction of the bulk pore water (and therefore solvated GPS) to move relatively quickly. Therefore, the observed rapid breakthrough is attributed to physical properties of the system rather than chemical properties. These results are in contrast with other GPS laboratory miscible displacement studies. Using a flow rate greater than those employed here and applying a continuous pulse, Magga et al. [18] reported a lag time of over 35 days before any GPS was detected in the effluent solution despite input concentrations twice those used in the current study. However, columns were 80 cm long and, therefore, the proportion of relatively non-tortuous path lengths is expected to be less. Additionally, analytical methods allowed for a lower limit of detection of 0.75 mg L⁻¹, so very low breakthrough concentrations may have gone undetected. Similarly, Beltran et al. [16] observed a lag time of approximately 300 pore volumes before GPS breakthrough was detected, although input concentrations were two orders of magnitude lower than those used in this study. Consistent with our results however, Dousset et al. [49] observed rapid breakthrough of GPS from vineyard soils with high concentrations of Cu when applying a finite pulse.

3.2.3. Multi-Reaction Transport and Linear Modeling

Measured GPS breakthrough curves along with MRTM and CXTFIT modeling simulations are displayed in Figure 5, with optimized parameter values and evaluative statistics given in

Table 5. Optimized parameter values for CXTFIT and MRTM models.

	Soil	R ± 95% CI	µ (h−1) ± 95% CI		RMSE	r2
CXTFIT	Commerce	6.645 ± 0.491	0.399 ± 0.007		0.005	0.827
	Sharkey	8.690 ± 0.828	0.926 ± 0.015		0.003	0.601
	Soil	Kirr (h−1) ± 95% CI	K1 (h−1) ± 95% CI	K2 (h−1) ± 95% CI	RMSE	r2
MRTM	Commerce	0.349 ± 0.003	0.434 ± 0.050	0.045 ± 0.006	0.002	0.965
	Sharkey	0.673 ± 0.015	0.255 ± 0.033	0.028 ± 0.008	0.002	0.901

. Upon statistical evaluation of several versions of MRTM, a two site multi-reaction model incorporating reversible and irreversible kinetic sites (Equations (6) and (9)) provided the best description of the data. In general, MRTM was able to describe observed data quite well from both soils with r² values of 0.97 and 0.90 for the Commerce and Sharkey soils, respectively. Additionally, an ocular assessment of Figure 5 indicates that the model was able to predict the overall shape of each BTC to a reasonable extent as well. Optimized parameter values indicate an order of magnitude higher rate of mass transfer to the solid phase relative to rates of release, with higher rates of sorption onto irreversible sites for the Sharkey soil and similar rates for the Commerce soil. Since GPS retention in both soils was so high and no extensive tailing occurred, this result is expected.

In general, linear modeling was able to predict the timing and magnitude of the effluent peak to some degree, although peak effluent concentrations were under predicted for both cases. Because the linear model does not account for non-equilibrium release from the solid phase, predicted concentrations in the effluent at the advanced stages of leaching were under predicted. Rates of mass transfer to the conceptual ‘sink’ (values for µ in Equation (3)) were very similar to rate coefficients determined by MRTM, which is expected as mass retention is a dominant mechanism in GPS transport within both soils. Overall, MRTM performed better than linear modeling due to its capability to account for a variety of retention mechanisms.

The MRTM model used in this study is a model that is simpler than those used by Candela et al. [17] and Magga et al. [18]. These scientists both employed models consisting of equilibrium and kinetic sites along with first order sinks and were able to describe GPS breakthrough reasonably well. Zhou et al. [19] also used a two-site non-equilibrium model to describe GPS breakthrough, however the dataset was very small (six points) and only adsorption was considered.

3.2.4. Distribution in the Soil Column

Measured GPS distribution obtained from KOH extractions of column sections along with MRTM and CXTFIT predicted distributions are displayed in Figure 6. Consistent with what would be expected for a strongly sorbing solute, the majority of the extracted mass is concentrated near the input port. In fact, 76% and 59% of total herbicide extracted from the columns was recovered from the first 2 cm for the Sharkey and Commerce soils, respectively. Larger quantities of retained GPS observed at lower depths in the Commerce column are reflective of greater herbicide mobility in this soil relative to the Sharkey soil. This finding is in agreement with our batch and BTC results. Recovered mass decreases rapidly with depth, with only minimal GPS recovery from the latter half of the columns for both soils. This distribution profile is consistent with those reported by a number of other studies. While conducting a mobility study through undisturbed soil columns, Okada et al. [47] recovered 68% of applied GPS in the top third of a 15 cm column, where Yang et al. [50] recovered the majority of GPS and AMPA residues from the upper 2 cm of soil in a field plot study. Landry et al. [15] reported that although no residual GPS was extracted after a yearlong field study involving agricultural soils in France, residual AMPA was concentrated in the top half of 20 cm profiles. Similar GPS distribution profiles were also described by Al-Rajab et al. [13], where the vast majority of GPS was extracted from the top 5 cm of a 25 cm profile of all three soils at all seven sampling times used in the study.

The CXTFIT model fails to predict the distribution of retained GPS in either column, with very low estimates of residual mass. This can be attributed to how the parameter µ is interpreted within the context of Equation (3). As discussed by van Genuchten et al. [51], µ is taken as a first-order degradation coefficient when utilizing CXTFIT within the STANMOD software package, and not as a rate coefficient for irreversible retention reactions. Low recovery of applied GPS in the effluent solution is therefore accounted for by relatively large optimized values for µ. Since this is a degradation rate coefficient, residual GPS within the column is not conserved, but is rather taken as mass lost from the system resulting in very low estimates of residual concentrations. Conversely, MRTM predicts the general shape of the distribution curve well, with greater amounts of residual herbicide located in the portion of the column closest to the inlet. The increased performance of this model relative to that of CXTFIT is due to complete conservation of mass throughout the duration of the numerical simulation. Here, there is no mechanism present to account for degradation and any type of mass ‘sink’ within the system is attributed to irreversible reactions. Although the shape of the distribution is approximated to a high degree, MRTM results in a somewhat over-prediction of measured residues, the reason for which is twofold. Extraction efficiencies from the soil are expected to be less than 100%, so measured residual GPS will automatically be less than what actually exists. Additionally, degradation of solid phase GPS is expected, the consequence of which being that mass loss due to biological activity is unaccounted for by MRTM, which will bias predictions higher than the actual amount of residual herbicide. Further effort to modify this model such that degradation is accounted for by the incorporation of various biological functions would improve estimates of residual mass. It is important to note that extraction data is not included in the model optimization procedure, and therefore, the ability of MRTM to estimate the general shape of GPS distribution within the soil profile further lends to the mechanistic validity of the model.

Discrepancies between modeled and measured results brought about by extraction inefficiencies and degradation can be handled by normalizing both data sets based upon calculated center of mass (COM) of GPS in the column, the results of which are given in

Table 6. MRTM predicted and measured center of mass for GPS distribution within Commerce and Sharkey soil columns.

Soil	MRTM COM (cm from Inflow)	Measured COM (cm from Inflow)
Commerce	2.46	2.64
Sharkey	2.29	1.48

. In order to do this, it must be assumed that extraction efficiency and rate of degradation from each column section is the same. The validity of these assumptions is uncertain, as it was determined in a separate study that degradation rates are dependent upon sorbed phase concentrations in the Commerce soil and that no clear trend was evident in the Sharkey soil (Unpublished results). In addition, extraction efficiencies from soil with differential sorbed phase concentrations are expected to be different, as extraction from soils with low sorbed phase concentrations will most likely be lower due to GPS association with higher affinity sites. Yet, this assumption may be valid for the Sharkey soil, as there was a more homogenous distribution of reactive site affinities (Freundlich n closer to 1), whereas the effect of differential extraction efficiencies will be greater in the Commerce soil due to a more heterogeneous distribution of reactive site affinities. However, if it is taken that these assumptions are valid, MRTM over-predicts GPS mobility in the Sharkey soil (predicted COM of 2.29 cm vs. a measured COM of 1.48 cm), and under-predicts the mobility of GPS in the Commerce soil (predicted COM of 2.46 cm vs. a measured COM of 2.64 cm). Additionally, a measured COM deeper in the soil profile for the Commerce soil relative to the Sharkey soil is consistent with observed BTC data and batch sorption results.

As biological degradation of GPS in soils is expected to occur over the time period in which these studies were conducted [44], efforts were made to quantify the amount of the primary metabolite (AMPA) in the residual extracts. A lack of prolonged measured radioactivity in the effluent solution suggests that the ¹⁴C remains associated with the highly reactive phosphonomethyl functional group, indicating that the dominant mechanism of microbial degradation is through the AMPA pathway, consistent with the findings of others [44,52]. Assuming that degradation beyond AMPA will result in a metabolite lacking the phosphonomethyl group and, therefore, a compound that will be readily mineralized to ¹⁴CO₂, all measured radioactivity in the extracting solution is taken to be either GPS or AMPA. As such, fractions of both GPS and AMPA determined via UPLC-MS/MS were applied to concentrations determined by LSC to produce GPS and AMPA distribution profiles displayed in Figure 7. Again, assuming that rates of degradation are identical throughout the soil profile despite differential solid phase concentrations of GPS, determining the COM of both compounds provides a basis to assess the mobility of AMPA relative to that of GPS in each soil. Calculated COM based on UPLC-MS/MS analyses are given in

Table 7. Calculated COM for GPS and AMPA within the Commerce and Sharkey soil columns.

Soil	GPS COM (cm from Inflow)	AMPA COM (cm from Inflow)
Commerce	2.55	3.80
Sharkey	1.40	1.95

. These results coupled with the above assumption indicate that AMPA is more mobile than GPS in both soils, although to a lesser extent in the Sharkey soil. This is consistent with the findings of Báez et al. [53], where reported K_f values for AMPA were lower than those for GPS in six out of eight soils studied. However, these conclusions must be considered carefully, as they are contingent upon the validity of a number of assumptions.

4. Conclusions

Batch studies of GPS sorption onto Commerce and Sharkey soils indicate a high affinity of both soils for the herbicide, although Sharkey exhibited a greater affinity. The high GPS retentive behavior of the Commerce soil is most likely due to the greater presence of Fe and Al amorphous oxides, whereas the high CEC of the Sharkey soil is likely responsible for GPS affinity in this matrix. Miscible displacement studies indicate that GPS mobility is highly limited in both Commerce and Sharkey soils, although the herbicide is slightly more mobile in the Commerce soil. A two-site multi-reaction model consisting of reversible and irreversible kinetic sites provided an adequate description of GPS breakthrough data from both soils. Predicted GPS breakthrough from both soils employing linear modeling adequately predicted the timing of peak effluent concentrations, but underestimated the magnitude of the peak in both cases. In general, MRTM outperformed linear approaches to describe GPS BTCs from both soils. Additionally, MRTM was able to accurately predict the general shape of the distribution of residual herbicide in the soil column, with residual mass concentrated towards the inflow point. Characterizing predicted and extracted residual herbicide based on COM calculations indicate that MRTM predictions overestimate the mobility of GPS in the Sharkey soil, whereas it was underestimated in the Commerce soil. Lack of prolonged release of radiolabeled compounds strongly suggests that microbial degradation of GPS is via the AMPA pathway, consistent with results reported by other authors. Assuming that degradation rates of solid phase GPS are constant, a deeper COM of residual AMPA in both soils is an indication of greater mobility of the primary GPS metabolite in the subsurface environment. Subsequent studies involving the development of models that consider additional processes will likely prove useful. For example, the incorporation of additional parameters that account for the degradation of solution and sorbed phase herbicide may improve the description of GPS breakthrough from soil columns, as well as the distribution of residual GPS and AMPA.