Introduction
Paraquat (PQ) (1,1′-dimethyl-4,4′-dipyridylium dichloride) is one of the most widely used herbicides in the world (Santos et al., Reference Santos, Schaule, Alves and Madeira2013; Rashidipour et al., Reference Rashidipour, Maleki, Kordi, Birjandi, Pajouhi, Mohammadi, Heydari, Rezaee, Rasoulian and Davari2019). It is favored for its fast-acting, non-selective contact killing of weeds, as well as its inactivation upon contact with soil due to rapid and strong adsorption (Wang et al., Reference Wang, Orr, He, Dalaijamts, Chiu, Tamamis and Phillips2019). During the adsorption process, soil particles, including organic matter and clay minerals, play a dominant role in binding PQ (Knight and Tomlinson Reference Knight and Tomlinson1967; Burns et al., Reference Burns, Hayes and Stacey1973). The strong binding impedes the leaching of paraquat, and inhibits its breakdown or utilization by microorganisms in the soil solution (Roberts et al., Reference Roberts, Dyson and Lane2002; Gondar et al., Reference Gondar, López, Antelo, Fiol and Arce2012), leading to its persistence and potential entry into the food chain through biomagnification, ultimately posing a health risk to humans and mammals upon long-term exposure (Pateiro-Moure et al., Reference Pateiro-Moure, Nóvoa-Munoz, Arias-Estévez, López-Periago, Martínez-Carballo and Simal-Gándara2009; Frimpong et al., Reference Frimpong, Ofori, Yeboah, Marria, Offeia, Apaataha, Sintima, Ofori-Ayeh and Osae2018). Research indicates that the immobility and persistence of PQ in soil is attributed primarily to its adsorption by clay, rather than organic matter (Khan, Reference Khan1974; Senesi et al., Reference Senesi, D’Orazio and Miano1995; Bromilow, Reference Bromilow2004). It was found that the adsorption capacity and binding strength of PQ on swelling clay (e.g. montmorillonite) far outweighs that on non-swelling clay (e.g. kaolinite) (e.g. Weber and Weed Reference Weber and Weed1968). This is due to attractive electrostatic interactions between the negatively charged clay layers and the positively charged PQ.
Several scholars have conducted experimental studies on the adsorption of PQ on montmorillonite (Mnt) (Weber et al., Reference Weber, Perry and Upchurch1965; Knight and Denny, Reference Knight and Denny1970; Tsai et al., Reference Tsai, Lai and Hsien2003; Tsai et al., Reference Tsai, Lai and Hsien2004; Muhamad et al., Reference Muhamad, Ismail, Sameni and Mat2011). The results obtained from these studies have revealed that exposure time has no significant impact on the adsorption of PQ on montmorillonite. On the other hand, solution pH and alkali metal ion concentration have a substantial effect on the adsorption process. Regarding the arrangement of PQ in the montmorillonite interlayer, previous studies have demonstrated that PQ formed a monolayer in parallel with montmorillonite layers in the anhydrous state (Knight and Denny, Reference Knight and Denny1970; Raupach et al., Reference Raupach, Emerson and Slade1979). However, the structures and dynamics of hydrated PQ-Mnt are yet to be explored.
Over the past 30 years, molecular simulation has shown unique advantages in uncovering the microstructure and dynamics of cations-clay and organic-clay systems, providing an important complement to experiments (Cygan et al., Reference Cygan, Guggenheim and van Groos2004a; Greenwell et al., Reference Greenwell, Jones, Coveney and Stackhouse2005; Liu et al., Reference Liu, Lu, Wang, Zhou and Xu2007; Cygan et al., Reference Cygan, Greathouse, Heinz and Kalinichev2009; Liu et al., Reference Liu, Lu, Wang, Zhou and Xu2009; Zhu et al., Reference Zhu, Shen, Ma, Ma, Zhou, Yuan, Liu and He2012a; Ngouana and Kalinichev, Reference Ngouana and Kalinichev2014; Teich-McGoldrick et al., Reference Teich-McGoldrick, Greathouse, Jové-Colón and Cygan2015; Scholtzová, Reference Scholtzová, Madejovaá, Jankovič and Tunega2016; Scholtzová, Reference Scholtzová, Cavallaro, Fakhrullin and Pasbakhsh2020; Cygan et al., Reference Cygan, Greathouse and Kalinichev2021). Various force fields have been utilized to describe the interactions between herbicides and clay minerals. For example, Cosoli et al. (Reference Cosoli, Fermeglia and Ferrone2010) employed the COMPASS force field to investigate the diffusion coefficient of atrazine herbicides in a saturated sand model. The results matched the experimental results obtained by Helmke et al. (Reference Helmke, Simpkins and Horton2005). The CHARMM (Brooks et al., Reference Brooks, Brooks, MacKerell, Nilsson, Petrella, Roux, Won, Archontis, Bartels, Boresch, Caflisch, Caves, Cui, Dinner, Feig, Fischer, Gao, Hodoscek, Im, Kuczera, Lazaridis, Ma, Ovchinnikov, Paci, Pastor, Post, Pu, Schaefer, Tidor, Venable, Woodcock, Wu, Yang, York and Karplus2009) force field and INTERFACE (Heinz et al., Reference Heinz, Lin, Mishra and Emami2013) force field, which describe organic materials and inorganic materials, respectively, were adopted to probe the adsorption mechanisms of PQ and glyphosate on montmorillonite (Wang et al., Reference Wang, Orr, He, Dalaijamts, Chiu, Tamamis and Phillips2019). By designing the variation in the protonation state of siloxane groups on the clay mineral surface at pH2 and pH7, and analyzing the herbicides’ residence time on the surface under different pH conditions, Wang et al. (Reference Wang, Orr, He, Dalaijamts, Chiu, Tamamis and Phillips2019) discerned that the interlayer adsorption configuration of PQ remains unaffected by pH fluctuations. In the study conducted by Zhu et al. (Reference Zhu, Shen, Ma, Ma, Zhou, Yuan, Liu and He2012a), ClayFF (Cygan et al., Reference Cygan, Liang and Kalinichev2004b) and CVFF (Dauber-Osguthorpe et al., Reference Dauber-Osguthor, Roberts, Osguthorpe, Wolff, Genest and Hagler1988) force fields were employed to investigate the adsorption of tetrachlorodibenzo-p-dioxin on tetramethylammonium and tetrapropylammonium modified montmorillonite. The findings indicated that dioxin exhibited a stronger interaction with tetra-alkylammonium species than with the clay surfaces. Additionally, the calculated adsorption energy of dioxin on the montmorillonite surface was comparable to that obtained from ab initio calculations (Austen et al., Reference Austen, White, Marmier, Parker, Artacho and Dove2007). These results suggest that the combined ClayFF–CVFF force field is suitable for herbicide–montmorillonite systems.
The aim of this study was to characterize in detail the microscopic structures and mobility of PQ–Mnt. Systematic molecular dynamics simulations were performed for the PQ–Mnt systems with varying levels of water content. The swelling curves and swelling thermodynamics were obtained. The interlayer configurations of PQ and water molecules as well as PQ’s mobility at different levels of water contents were characterized. The findings in the present study provide an atomic-level understanding for the environmental fate of PQ in soils.
Materials and methods
Models
Two montmorillonite models with different chemical compositions were built in this study. One had the same composition with Ilari et al. (Reference Ilari, Etcheverry, Waiman and Zanini2021), i.e. Na0.84[Si7.74Al0.26][Al2.88Fe3+0.54Mg0.58]O20(OH)4 (denoted as 0.84-Mnt). The other was a common Fe3+-free montmorillonite with a chemical formula of Na1.0[Si8][Al3.0Mg1.0]O20(OH)4 (denoted as 1.0-Mnt).
The supercell used in this study consisted of two clay platelets with each containing 24 unit cells: six in the x-dimension and four in the y-dimension (Fig. 1a,b). In the 0.84-Mnt system, 14 Mg2+ for Al3+ and 13 Fe3+ for Al3+ isomorphic substitutions in octahedral sheet and 6 Al3+ for Si4+ isomorphic substitutions in tetrahedral sheet were introduced in each layer, resulting in a layer charge density of –0.84 e/unit cell. In the 1.0-Mnt system, 24 isomorphic substitutions of Mg2+ for Al3+ in octahedral sheet were imposed in each layer, which resulted in a layer charge density of -1.0 e/unit cell. The isomorphic substitutions obey Loewenstein’s rule, i.e. two substitution sites cannot be adjacent (Loewenstein Reference Loewenstein1954). Figure 1c shows the molecular structure of PQ cation (C12H14N2)2+. Initially, we placed 8 PQ cations randomly in each interlayer region of the 0.84-Mnt system (denoted as PQ-0.84-Mnt) according to the experimental loadings in Ilari et al. (Reference Ilari, Etcheverry, Waiman and Zanini2021). The samples measured in the experiment are considered to be anhydrous due to the drying process. Na+ were used to compensate the net charges of the system. For the 1.0-Mnt system, we placed 12 PQ cations randomly in each interlayer region to compensate the net layer charge, denoted as PQ-1.0-Mnt. To study the effect of water content on the interlayer structure, a different number of water molecules were inserted into the interlayer region of PQ-1.0-Mnt (Nwater=~0–600). For reference purposes, the water contents were converted to mgwater/gclay, e.g. the water content for N=600 is 555 mgwater/gclay, and the system is thus denoted as PQ-555. The compositions of paraquat-0.84-Mnt and paraquat-1.0-Mnt systems are shown in Table S1.
(a) Snapshot of the anhydrous PQ-0.84-Mnt; (b) snapshot of the paraquat-1.0-Mnt at 278 mgwater/gclay; (c) structure of paraquat cation. Key: C, green; N, blue; H, white; O, red; Si, yellow; Mg, cyan; Al, pink; Fe, violet.
Computational details
All molecular dynamics simulations (MD) were carried out using the LAMMPS package (Plimpton, Reference Plimpton1995). Clayff-CVFF was employed to describe the interactions in the systems. The partial charges of PQ (Table S1; Fig. S1) were obtained using the charge equilibration method (QEq) (Rappe and Goddard, Reference Rappe and Goddard1991). The flexible single point charge (SPC) water model (Berendsen et al., Reference Berendsen, Postma, van Gunsteren, Hermans and Pullman1981; Cygan et al., Reference Cygan, Liang and Kalinichev2004b) was used to describe the water molecule. A 10.0 Å cut-off was used for the short-range interactions. The long-range electrostatic interactions were calculated using the Ewald summation method (Frenkel and Smit, Reference Frenkel and Smit2002) and the number of K-space vectors was determined to reach a precision of 1.0×10–4.
The isothermal-isobaric (constant particle number, pressure, and temperature; NPT) simulations were performed at 298K and 1 atm for 20 ns to obtain the basal spacing (Fig. S2). The equilibrium was confirmed by monitoring the evolution of simulation box and energy. Nosé-Hoover thermostat and barostat for temperature and pressure controls were used in this study, and the coupling constants were set to 0.1 ps and 1.0 ps, respectively. The equilibrium was checked by monitoring the evolution of simulation box and energy. The microscopic structures, immersion energy and dynamics of interlayer species were obtained from a 30 ns production run in the canonical ensembles (constant particle number, volume, and temperature; NVT) at 298K following the previous NPT simulations. In all simulations, the time step was set to 0.5 fs and an interval of 500 fs was used to record trajectories. All atoms were allowed to move freely.
Data analysis
The distribution of interlayer species in space is characterized by the atomic number density distribution along the z direction (normal to the clay basal plane), which is calculated by averaging the NVT trajectories. Consider the plane where the oxygen atoms in the silicon-oxygen ring at the bottom of the montmorillonite layer are as the starting plane (z=0).
The formula for distribution of species number density between layers is as follows:
$$ \unicode{x03C1} (z)=\frac{\left\langle \mathrm{N}\left(z-\frac{\Delta z}{2},z+\frac{\Delta z}{2}\;\right)\right\rangle }{\Delta \mathrm{z}\times \mathrm{S}}. $$
Here
$ \mathrm{N}\left(z-\frac{\Delta z}{2},z+\frac{\Delta z}{2}\;\right) $
is the average number of atoms appearing in the height interval
$ \left(z-\frac{\Delta z}{2},z+\frac{\Delta z}{2}\;\right) $
(∆z = 0.2 Å in this work).
The self-diffusion coefficient D of interlayer species can be obtained using the Einstein relationship (Allen & Tildesley, Reference Allen and Tildesley1987; Chang et al., Reference Chang, Skipper and Sposito1997):
$$ \frac{1}{N}\sum \limits_{i=1}^N\left\langle {\left|{r}_i(t)-{r}_i(0)\right|}^2\right\rangle =2 dDt, $$
where N is the number of atoms, ri (t) is the center-of-mass position of the ith atom at time t, and d is the diffusion dimension. For example, d = 3 represents the total coefficient, and d = 1 represents the coefficient components in the x, y, and z directions. The left-hand side of Eqn (2) is the mean square displacement term (MSD). The components of MSD in each individual direction are denoted as XX, YY, and ZZ. The motion parallel to the x–y plane (denoted as ||) is calculated as (XX+YY)/2. The standard error was taken as the error. The self-diffusion coefficients were calculated from the MSD within 2 ns intervals in the linear portion of the MSD curve, which spanned from 2 ns to 18 ns. We divided this 16 ns period into eight equal segments and utilized these segments to calculate the standard error. Given that simulation size may have an impact on the diffusion coefficient, we utilized the correction method proposed by Simonnin et al. (Reference Simonnin, Noetinger, Nieto-Draghi, Marry and Rotenberg2017) in our study. The interlayer water viscosity, as provided by Ho et al. (Reference Ho, Wang, Jové Colón and Coker2020), was selected for calculation purposes. In the PQ-185 systems, the original self-diffusion coefficient of water was 3.388×10–10 m2 s–1 with a correcting amplitude of –5.3×10–12 m2 s–1. In the PQ-278 systems, the original self-diffusion coefficient of water was 5.663×10–10 m2 s–1 with a correcting amplitude of –6.7×10–12 m2 s–1. The magnitude of the correction in both systems is considered negligible and can be disregarded.
Results and Discussion
Swelling behaviors
In the simulated swelling curve of PQ-1.0-Mnt (Fig. 2), the experimental value for the anhydrous PQ-0.84-Mnt is also shown for comparison. It can be observed that the experimental value is very close to our calculated basal spacing, i.e. 12.5 Å vs 12.7 Å. This agreement indicates that the ClayFF-CVFF force field can describe the swelling behavior of the PQ-Mnt system. Swelling thermodynamics is described using the immersion energy (Boek et al., Reference Boek, Coveney and Skipper1995; Smith, Reference Smith1998), which is defined as:
$$ Q=U(N)-U\left({N}_0\right)-\left(N-{N}_0\right){U}_{\mathrm{bulk}}, $$
Simulated swelling curve of PQ-1.0-Mnt. For comparison, the calculated basal spacing for the anhydrous PQ-0.84-Mnt is also presented alongside the corresponding experimental value (Ilari et al., Reference Ilari, Etcheverry, Waiman and Zanini2021).
where U (N) is the potential energy of hydrated montmorillonite, U (N 0) is the energy of the reference state (370 mgwater/gclay in this study), and U bulk is the mean interaction potential of the bulk flexible SPC water (−41.12 kJ mol–1) (Teleman et al., Reference Teleman, Jönsson and Engström1987), which was calculated from the simulation of a pure water system. The immersion energy Q represents the energy released when the system with N water molecules is transformed into the reference state through water absorption from the bulk water.
In the calculated immersion energy curve of the PQ-1.0-Mnt systems (Fig. 3), two minima can be identified on the curve: local minimum at 185 mgwater/gclay (200 water molecules in each interlayer region) and global minimum at 278 mgwater/gclay (300 water molecules in each interlayer region). These minima suggest that PQ-1.0-Mnt with the corresponding water contents are thermodynamically stable and could exist in natural or experimental conditions. Therefore, these two hydrated systems and the anhydrous system were selected to investigate the interlayer structures and dynamics in the following sections. The energy barrier for the transformation from the 185 mgwater/gclay state to the 278 mgwater/gclay state was estimated to be ~3.0 J gclay–1, indicating that PQ-1.0-Mnt can swell. (Information on the interlayer structure when the water content was <100 mgwater/gclay is provided in the Supplementary material; see Fig. S3.)
Simulated immersion energy curve of paraquat-1.0-Mnt. The statistical error was estimated as the standard error of average values obtained from 10 evenly divided fragments of the trajectory (each spanning 3 ns intervals).
Interlayer structure
To investigate the arrangement of interlayer species, we calculated the atomic number density along the z-axis direction, which was obtained by averaging the density in the two interlayer regions. The PQ cations were represented by N and C atoms, and water molecules were represented by O atoms. In the anhydrous system, PQ presented a sharp peak in the density profile (Fig. 4a), suggesting that PQ formed a monolayer, with the bipyridine ring plane parallel to the montmorillonite surface (Fig. 5a). In PQ-185 and PQ-278 systems, three peaks can be identified for water molecules (Fig. 4b,c), suggesting a triple-layer-like structure. PQ cations show two sharp peaks in the two hydrated states. It was found from the trajectory that PQs are in direct contact with the clay mineral surface, and there are no water molecules between them.
Z-density distribution of the paraquat and water molecules in the PQ-1.0-Mnt systems with different water contents (the gray line represents the location of the inner surface of the montmorillonite layer). (a) PQ-0 system; (b) PQ-185 system; (c) PQ-278 system.
Snapshots of the PQ-1.0-Mnt systems at different water contents. (a) PQ-0 system; (b) PQ-185 system; (c) PQ-278 system. Key: C, green; N, blue; H, white; O, red; Si, yellow; Mg, cyan; Al, pink. Hydrogen atoms of paraquat molecules are not shown for clarity.
As presented above, PQ cations adopted a parallel orientation to the clay surface in both anhydrous and hydrated states. The PQ molecules were nearly parallel to the clay surface. In all systems, the angles between the two PQ pyridine rings peaked at approximately 14–30° (Fig. 6). Due to the influence of clay surfaces and water molecules on PQ cations, the angle between the pyridine rings deviated from the equilibrium of 30° (Mendes et al., Reference Mendes, de Freitas, Brown and de Souza2020; Hou and Wang, Reference Hou and Wang2021; Macii et al., Reference Macii, Detti, Bloise, Giannarelli and Biver2021). The hydration shell of PQ in hydrated PQ-Mnts were depicted using the spatial distribution function of water around PQ (Fig. 7). In both hydration states, there is no water distribution between PQ molecules and montmorillonite basal surface, indicating that PQ is in direct contact with montmorillonite. However, no correlation between clay substitutions and PQ location was found from the trajectories. In PQ-185, no water distributed above PQ molecules along the z direction (the hydration shell exists only in the xy plane), which suggests that PQ is also in direct contact with other PQ molecule. Together with the equilibrium structure, this observation suggests that PQ molecules formed stable π-π stacking in PQ-185. In PQ-278, the hydration shell is closed above the PQ molecule, which is consistent with the density distribution curve in Fig. 4c where the PQ molecules were separated by a water layer.
Probability distribution of angles between pyridine rings connected to each other within PQ cations in different systems.
The spatial distribution maps of water molecules (hydration shell) around PQ on the lower surface of the montmorillonite. (a) Front view for PQ-185 system; (b) bottom view for PQ-185 system; (c) front view for PQ-278 system; (d) bottom view for PQ-278 system. Key: C, green; N, blue; H, white; Si, yellow; Mg, cyan; Al, pink; O, red. Clay surfaces are represented by thin lines.
Mobility of interlayer species
To investigate the mobility of PQ, we calculated the self-diffusion coefficient of PQ at different water contents (Fig. 8) (see Supplementary material for more details). The diffusion coefficients of PQ at these water contents are all of the magnitude of 10–13 to 10–11 m2 s–1. They are much smaller than the self-diffusion coefficient (magnitude of 10–9–10–8 m2 s–1) of various polycyclic aromatic hydrocarbons (e.g. PAHs, including naphthalene, anthracene, chrysene, and benzopyrene) intercalated in clay minerals in a wide range of water contents (Chen and Hu, Reference Chen and Hu2022; Zhao et al., Reference Zhao, Tan, Zhang, Zhen, Song, Ju and Ling2023), and it is also significantly smaller than the self-diffusion coefficient of PQ in bulk water (~10–10 m2 s–1).
The self-diffusion coefficient (D||) of PQ in PQ-1.0-Mnt at different water contents.
As PQ cations almost only move in the xy plane direction, the self-diffusion coefficients diagram (Fig. 8) revealed that the self-diffusion coefficient (D ||) were similar for the anhydrous state and PQ-185, and both were significantly smaller than that at PQ-278 that lacked π-π stacking interactions. These results indicate that PQ has very limited movement in the interlayer region and can be effectively trapped by montmorillonite. It is consistent with previous conclusion (e.g. Smith and Mayfield, Reference Smith and Mayfield1978; Leonard et al., Reference Leonard, Langdale and Fleming1979; Ouyang et al., Reference Ouyang, Mansell and Nkedi-Kizza2004; Amondham et al., Reference Amondham, Parkpian, Polprasert, DeLaune and Jugsujinda2006) that the transport of paraquat is extremely limited in clay-rich soils. In the PQ-175 and PQ-278 systems, the self-diffusion coefficients of water molecules increased with water content, and was much smaller than the mobility of bulk water (Mark and Nilsson, Reference Mark and Nilsson2001) (Fig. 9). Moreover, these coefficients were found to be within the same order of magnitude as the self-diffusion coefficient of interlayer water in Na and Ca-montmorillonite with two or three layers of water (Greathouse et al., Reference Greathouse, Hart, Bowers, Kirkpatrick and Cygan2015; Zhang et al., Reference Zhang, Lu, Liu, Zhou and Zhou2014).
The self-diffusion coefficient (Dtotal) of water in PQ-185, PQ-278, and bulk water.
Conclusion
Paraquat is used widely in agriculture production, and poses a health risk to humans due to its long-term stability in soils. In this study, molecular dynamics simulations were conducted to explore the structure and dynamics of PQ in the interlayer space of montmorillonite. Two stable hydration states were identified, one with a water content of 185 mgwater/gclay and the other of 278 mgwater/gclay, with the latter being the most stable. The interlayer structures and mobility of PQ in the stable hydration and anhydrous states were obtained. It was found that PQs were in direct contact with clay surfaces in both anhydrous and hydrated states. At the water content of 185 mgwater/gclay, PQs formed π-π stacking while at the water content of 278 mgwater/gclay, the PQs were separated by a layer of water. The self-diffusion coefficient of PQ ranges from 10–12 to 10–11 m2 s–1, with a slightly greater value observed at the high water content. Notably, the diffusivity of PQ is significantly smaller than that of other aromatic organic compounds. The findings presented in this study provide an atomic-level insight into the stability of PQ in soils and forms a microscopic basis for development of efficient remediation strategy.
Supplementary material
The supplementary material for this article can be found at http://doi.org/10.1017/cmn.2024.13.
Data availability statement
The datasets generated and/or analyzed during the current study are available from the authors upon reasonable request.
Acknowledgements
We are grateful to the High-Performance Computing Center (HPCC) of Nanjing University for doing the numerical calculations in this paper on its blade cluster system.
Author contribution
Writing-original draft, Visualization, Methodology, Investigation: Haotian Su; Writing-review and editing, Investigation: Yingchun Zhang;Writing-review and editing, Investigation, Supervision: Xiandong Liu; Writing-review and editing: Qingfeng Hou, Xiancai Lu.
Financial support
This study was supported by the National Natural Science Foundation of China (grant nos. 42125202 and 42172050).
Competing interests
The authors declare that they have no competing interests.
