A molecular dynamics simulation study of paraquat intercalated montmorillonite

Abstract

Paraquat, one of the most widely utilized herbicides globally, causes a significant environmental challenge due to its poor degradation rate and tendency to adsorb into clay interlayers. Several remediation methods have been proposed but their effectiveness remains suboptimal. The primary reason for this is the lack of microscopic understanding of paraquat–montmorillonite interactions. In this work molecular dynamics simulations were applied to study the interlayer structures and mobility of paraquat intercalated montmorillonite. Two stable hydration states were identified from the calculated immersion energy curve, which corresponded to a water content of 185 mg_water/g_clay and 278 mg_water/g_clay (the most stable). Paraquats remained in direct contact with the clay surface in both the anhydrous and hydrated states. At the water content of 185 mg_water/g_clay, paraquats formed π-π stacking while at 278 mg_water/g_clay, they were separated by a layer of water. Paraquat showed very small self-diffusion coefficients in the interlayer space of montmorillonite, indicating rather limited motions. The results in this work provide a basis for a better understanding of the interaction of paraquat with clay minerals.

Introduction

Paraquat (PQ) (1,1′-dimethyl-4,4′-dipyridylium dichloride) is one of the most widely used herbicides in the world (Santos et al., 2013; Rashidipour et al., 2019). 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., 2019). During the adsorption process, soil particles, including organic matter and clay minerals, play a dominant role in binding PQ (Knight and Tomlinson 1967; Burns et al., 1973). The strong binding impedes the leaching of paraquat, and inhibits its breakdown or utilization by microorganisms in the soil solution (Roberts et al., 2002; Gondar et al., 2012), 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., 2009; Frimpong et al., 2018). Research indicates that the immobility and persistence of PQ in soil is attributed primarily to its adsorption by clay, rather than organic matter (Khan, 1974; Senesi et al., 1995; Bromilow, 2004). 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 1968). 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., 1965; Knight and Denny, 1970; Tsai et al., 2003; Tsai et al., 2004; Muhamad et al., 2011). 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, 1970; Raupach et al., 1979). 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., 2004a; Greenwell et al., 2005; Liu et al., 2007; Cygan et al., 2009; Liu et al., 2009; Zhu et al., 2012a; Ngouana and Kalinichev, 2014; Teich-McGoldrick et al., 2015; Scholtzová, 2016; Scholtzová, 2020; Cygan et al., 2021). Various force fields have been utilized to describe the interactions between herbicides and clay minerals. For example, Cosoli et al. (2010) 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. (2005). The CHARMM (Brooks et al., 2009) force field and INTERFACE (Heinz et al., 2013) 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., 2019). 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. (2019) discerned that the interlayer adsorption configuration of PQ remains unaffected by pH fluctuations. In the study conducted by Zhu et al. (2012a), ClayFF (Cygan et al., 2004b) and CVFF (Dauber-Osguthorpe et al., 1988) 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., 2007). 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. (2021), i.e. Na_0.84[Si_7.74Al_0.26][Al_2.88Fe³⁺_0.54Mg_0.58]O₂₀(OH)₄ (denoted as 0.84-Mnt). The other was a common Fe³⁺-free montmorillonite with a chemical formula of Na_1.0[Si₈][Al_3.0Mg_1.0]O₂₀(OH)₄ (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 Mg²⁺ for Al³⁺ and 13 Fe³⁺ for Al³⁺ isomorphic substitutions in octahedral sheet and 6 Al³⁺ for Si⁴⁺ 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 Mg²⁺ for Al³⁺ 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 1954). Figure 1c shows the molecular structure of PQ cation (C₁₂H₁₄N₂)²⁺. 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. (2021). 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 (N_water=~0–600). For reference purposes, the water contents were converted to mg_water/g_clay, e.g. the water content for N=600 is 555 mg_water/g_clay, 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.

Figure 1. (a) Snapshot of the anhydrous PQ-0.84-Mnt; (b) snapshot of the paraquat-1.0-Mnt at 278 mg_water/g_clay; (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, 1995). 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, 1991). The flexible single point charge (SPC) water model (Berendsen et al., 1981; Cygan et al., 2004b) 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, 2002) 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:

(1) $$ \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, 1987; Chang et al., 1997):

(2) $$ \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, r_i (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. (2017) in our study. The interlayer water viscosity, as provided by Ho et al. (2020), was selected for calculation purposes. In the PQ-185 systems, the original self-diffusion coefficient of water was 3.388×10^–10 m² s^–1 with a correcting amplitude of –5.3×10^–12 m² s^–1. In the PQ-278 systems, the original self-diffusion coefficient of water was 5.663×10^–10 m² s^–1 with a correcting amplitude of –6.7×10^–12 m² 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., 1995; Smith, 1998), which is defined as:

$$ Q=U(N)-U\left({N}_0\right)-\left(N-{N}_0\right){U}_{\mathrm{bulk}}, $$

Figure 2. 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., 2021).

where U (N) is the potential energy of hydrated montmorillonite, U (N ₀) is the energy of the reference state (370 mg_water/g_clay 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., 1987), 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 mg_water/g_clay (200 water molecules in each interlayer region) and global minimum at 278 mg_water/g_clay (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 mg_water/g_clay state to the 278 mg_water/g_clay state was estimated to be ~3.0 J g_clay^–1, indicating that PQ-1.0-Mnt can swell. (Information on the interlayer structure when the water content was <100 mg_water/g_clay is provided in the Supplementary material; see Fig. S3.)

Figure 3. 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.

Figure 4. 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.

Figure 5. 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., 2020; Hou and Wang, 2021; Macii et al., 2021). 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.

Figure 6. Probability distribution of angles between pyridine rings connected to each other within PQ cations in different systems.

Figure 7. 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 m² s^–1. They are much smaller than the self-diffusion coefficient (magnitude of 10^–9–10^–8 m² 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, 2022; Zhao et al., 2023), and it is also significantly smaller than the self-diffusion coefficient of PQ in bulk water (~10^–10 m² s^–1).

Figure 8. 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, 1978; Leonard et al., 1979; Ouyang et al., 2004; Amondham et al., 2006) 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, 2001) (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., 2015; Zhang et al., 2014).

Figure 9. The self-diffusion coefficient (D_total) 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 mg_water/g_clay and the other of 278 mg_water/g_clay, 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 mg_water/g_clay, PQs formed π-π stacking while at the water content of 278 mg_water/g_clay, the PQs were separated by a layer of water. The self-diffusion coefficient of PQ ranges from 10^–12 to 10^–11 m² 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.