Adsorption process of polar and nonpolar compounds in a nanopore model of humic substances

Humic substances (HSs) as a major part of soil organic matter (SOM) are complex combinations of natural organic matter, which are ubiquitous in the environment and occur predominantly in soils, residues and natural water. They play a crucial role in the fate and behaviour of contaminants in the environment. In this work, an HS nanopore model was based on a structural motif containing polar and nonpolar domains. The polar domain was represented by the carboxyl groups that predominate in the composition of HSs, whereas the nonpolar domain was constructed by aliphatic chains. The effect of hydrophobic/hydrophilic interactions on sorption in HSs was analysed by the inclusion of the nonpolar polycyclic aromatic hydrocarbon naphthalene and the polar herbicide 2‐methyl‐4‐chlorophenoxyacetic acid (MCPA) in the nanopore. Water cluster was also included to allow the formation of the water molecular bridges (WaMB), which stabilized the whole system. Thermodynamic potential‐of‐mean‐force (PMF) calculations by means of classical molecular dynamics simulations were performed to investigate the effect of polar and nonpolar domains in the adsorption mechanism. The polar MCPA molecule was stabilized at the hydrophilic water–HS interface, whereas nonpolar naphthalene was stabilized inside of the hydrophobic nanovoid of the HS nanopore model. The obtained results demonstrated that the adsorption strength of HSs regarding polar and nonpolar species is similar, but the adsorption mechanisms differ.

• The work elucidates the adsorption mechanism of polar and nonpolar species in SOM at the molecular scale. • The derived HS model was suitable to describe and explain the adsorption properties of SOM by means of constrained molecular dynamics simulations. • Adsorption mechanisms for polar and nonpolar species in SOM differ. Polar molecules prefer the water-SOM interface, whereas nonpolar molecules are trapped in hydrophobic nanovoids of SOM. • Molecular dynamics simulation is an effective method providing a detailed description of the adsorption processes in SOM. Humic substances (HSs) as a dominant part of soil organic matter (SOM) are biogenic, heterogeneous, large, relatively stable organic species mainly found in soils and water residues (Aiken, Mcknight, Wershaw, & MacCarthy, 1985). Humic substances originate from the decomposition of dead biological material, mostly from plants, as well as from microbially-derived material and animals in the milieu matter (Filella, Parthasarathy, & Buffle, 2005;Rodriguez & Nunez, 2011). They play a fundamental role in many biogeochemical processes in natural systems and constitute the most abundant pool of carbon (C) on earth (Stevenson, 1994). It is argued that HSs control the fate of environmental contaminants through acting as sorbents for toxic metal ions, radionuclides and organic pollutants. Humic substances have a composition of a variety of organic compounds (dominantly humic and fulvic acids) bound by various intermolecular forces such as electrostatic, van der Waals and hydrophobic interactions, hydrogen bonds and cation bridges (Piccolo, 2001;Sutton & Sposito, 2005). The high variability of HSs, as a result of their composition depending on the specific ecosystem where they originate (e.g., the vegetation, climate and topography), is a complicating factor in the determination of their complete molecular structure. Therefore, the chemical composition and structure of HSs (or more generally SOM) remains an object of debate (Gerke, 2018;Lehmann & Kleber, 2015) and several models have been developed in the literature to describe principles of their structure and evolution. In early times, a polymer structure was proposed in which HSs were described as an assembly of several organic macromolecules with high molecular weights of up to even hundreds of thousands of grams per mole (Hayes & Clapp, 2001;Stevenson, 1994). However, the present hypothesis based on recent experimental findings defines HSs as heterogeneous supramolecular assortments of rather small molecules, with a molecular weight of a few thousand grams per mole stabilized mainly through noncovalent interactions, electrostatic and hydrogen bonds, and hydrophobic interactions. Starting from this hypothesis, HSs can also form self-assembled micellar constructions, with an inner hydrophobic region protected from outside water through interfacial hydrophilic regions (Nebbioso & Piccolo, 2012;Piccolo, 2001;Piccolo, 2002;Piccolo, Conte, & Cozzolino, 2001). Additionally, protection of HSs by mineral associations and HSs' surface adsorption on reactive soil surfaces is an important process in SOM dynamics and stabilization (Lehmann & Kleber, 2015).
Molecular simulations, which include classical force-field (FF)-based methods and/or quantum chemical (QC) methods, have become an essential tool in the theoretical investigation of the molecular structure and properties of HSs. They allow the determination of chemical, physical and electronic properties at the molecular level, as well as the interactions of HSs with other systems in the environment (e.g., organic contaminants) (Bruccoleri, Sorenson, & Langford, 2001;Haberhauer, Aquino, Tunega, Gerzabek, & Lischka, 2001;Kubicki & Apitz, 1999;Orsi, 2014;Schulten, 1999;Tunega, Gerzabek, Haberhauer, Lischka, & Aquino, 2016). Pioneering work in this area regarding the development of the macromolecular structure of HSs through simulations was based on classical FF simulations by Schulten and Schnitzer (1997). Computational tools were also used to develop "building-block" structures, such as the widely known Steelink model (Steelink, 1985) and the TNB (Temple-Northeastern-Birmingham) model (Davies et al., 1997). Diallo et al. (2003) used integrated experimental and computational data to build up threedimensional (3-D) structures of humic acids. Further, molecular mechanics based on experimental results were used to create a fulvic acid model taking the TNB monomer as a reference (Alvarez-Puebla, Valenzuela-Calahorro, & Garrido, 2006). Recently, a general HS building tool was developed: the Vienna Soil-Organic-Matter Modeler (VSOMM) (Petrov, Tunega, Gerzabek, & Oostenbrink, 2017;Sündermann et al., 2015). The VSOMM uses a database of molecular fragments generated on the base of knowledge on the elemental composition and structural/functional groups of various SOM species obtained from experiments. The VSOMM can create a supramolecular aggregate of the SOM, comprising thousands of atoms with a desired chemical and functional composition, pH, water and cation content. A comprehensive and detailed review of the application of molecular modeling simulations of HS models is available in the work by Orsi (2014).
The occurrence of HSs in soils has an important effect on the sorption of chemically very different contaminants, including pesticides, heavy metals, radioactive substances and other organic compounds (Stevenson, 1985;Sposito, 1986). The flexible structure, composition and organization of a form of supramolecular and organo-mineral aggregates are responsible for a high sorption capacity of soils. The unique properties of HSs are mainly due to the presence of different types of organic functional groups, from polar (carboxyl, hydroxyl, amino, methoxy and phenolic) to nonpolar (aliphatic and aromatic), which are responsible for the occurrence of hydrophilic and/or hydrophobic domains in the HSs' structures. It is known that carboxyl groups (in neutral or anionic form) belong to the most abundant and active functionalities of HSs (Leenheer, Wershaw, & Reddy, 1995a;Leenheer, Wershaw, & Reddy, 1995b), and are accountable, for example, for ion exchange properties of HSs in relation to metal ions. Despite the importance of HSs for the environment, the chemical and structural complexity of HSs as a supramolecular association of smaller molecules makes the modelling of molecular interaction with contaminants a difficult task. In addition, factors such as pH, water content and natural inorganic cations in SOM have a strong impact on sorption of various species, making simulations even more complicated. Several models of HSs characterized by a relatively low number of atoms have been used previously to analyze the interactions, formation energies and thermodynamic properties of adsorption processes (Ahmed et al., 2015;Gros, Ahmed, Kühn, & Leinweber, 2017;Kubicki & Apitz, 1999;Schulten, Thomsen, & Carlsen, 2001;Tunega et al., 2016). QC methods are rarely used in the simulations of HSs as they are time consuming, which limits the number of atoms in the models. Indeed, in our previous works on modelling HSs' interactions with organic pollutants, we have developed structurally simple models (functional and oligomeric) containing polar functional groups typical for HSs (Aquino et al., 2008;Aquino, Tunega, Haberhauer, Gerzabek, & Lischka, 2007).
In this work, we focus solely on HSs. Generally, it is practically impossible to create a reference model of HSs due to their extremely high compositional and structural variability. Additionally, there is a lack of complete structural data information. In this work, we built up a representative model in which polar (hydrophilic) and nonpolar (hydrophobic) domains, typical for HSs, were explicitly included. We developed an HS nanopore model containing polar and nonpolar domains. The goal of this work was to investigate the effect of hydrophobic/hydrophilic interactions on sorption of nonpolar polycyclic aromatic hydrocarbons (PAHs), particularly the naphthalene molecule and a polar herbicide molecule (2-methyl-4-chlorophenoxyacetic acid (MCPA)), by performing thermodynamic potential-of-mean-force (PMF) calculations using classical molecular dynamics simulations.

| METHODS
The model consists of three layers of connected aliphatic chains terminated by carboxyl groups. Figure 1 shows the front and side view of this model. Each layer consists of one aliphatic chain (15C length) forming a backbone of the model (bottom part in Figure 1). Two additional aliphatic chains (11C length) are connected to this chain at both ends in a nearly perpendicular configuration. The distance between these perpendicular chains is about 13 Å. Further, these chains have carboxyl groups connected to the last but one C in the chain in a face-to-face configuration ( Figure 1). The whole HS nanopore model was optimized prior to any molecular dynamics (MD) simulation to obtain optimal distances between the layers. The optimization of the distances between the layers of the HS model structures (d in Figure 1) was performed by using a density functional tight binding method (DFTB, Elstner et al., 1998) implemented in the DFTB+ code (Aradi, Niklasson, & Frauenheim, 2015). The distance between the layers was changed step by step with an interval of about 0.2 Å, keeping the terminal C atoms of the bottom aliphatic chain (15C) of all three layers fixed and reaching the optimized value of 4.7 Å. Our HS nanopore model explicitly includes three characteristic features typical for HSs: (a) hydrophilic and hydrophobic domains represented by the carboxyl groups (shown as oval in Figure 1) and by aliphatic chains (rectangle in  Figure 1), respectively; (b) such a geometric arrangement forms a nanopore space; and (c) the model also incorporates wetting conditions, represented by the inclusion of a certain amount of water molecules (microhydration), which are important for overall stabilization of the HSs' structures. The water molecules are usually concentrated in the hydrophilic domains of HSs, establishing among the polar functional groups (Aquino et al., 2009;Aquino et al., 2011a;Aquino et al., 2011b;Schaumann & Bertmer, 2008). In Figure 1, such a water molecular bridge is shown as a chain of water molecules connecting carboxyl groups through hydrogen bonds. A water molecular bridge can lock the nanopore space, forming a nanocavity in which host molecules can be trapped.
To simulate a microhydration effect on adsorption in the HS model, as mentioned above, a cluster of 50 water molecules was added near the polar carboxyl groups. The initial HS geometry is shown in Figure 2. The model was completed by the inclusion of two molecular species to demonstrate the capability of HSs to adsorb both polar and nonpolar species. Naphthalene was selected as a nonpolar moiety representing polycyclic aromatic hydrocarbons (PAHs), which are organic contaminants. MCPA was chosen as a representative polar organic pollutant, which can enter the environment in large amounts. MCPA is a widely used herbicide that can remain in relatively high concentrations in soils because its chemical decomposition is a relatively slow process. Therefore, it is important to know how MCPA behaves in soils and how it can interact with various soil components.
The simulations aimed to predict the change of free energy in a process of trapping MCPA/naphthalene molecules in the nanopore of the HS model ( Figure 2) by means of constrained classical molecular dynamics calculations. They were performed as follows. The adsorbed molecule was placed near the backbone of the aliphatic chains of the HS model. The C atoms located at the corner of the aliphatic backbone (yellow balls in Figure 3) were fixed to keep the model free from unexpected changes during the simulation steps. The free energy change arising from the interaction of the different molecular species was evaluated by means of potential of mean force (PMF) calculations (Roux, 1995). Details of the method used can be found in the Supporting Information. Briefly, in the PMF method, a free energy profile (FEP) is calculated with respect to a defined reaction coordinate (rc) of the system under investigation, which represents a constraint in the MD simulation. To evaluate the PMF from this statistical manifold, the Blue Moon method (Giovanni, Kapral, & Vanden-Eijnden, 2005;Sprik & Ciccotti, 1998) was applied in using a constraint imposed on the distance between the adsorbed molecule and the HS model. Particularly, the reaction coordinate was defined as a distance between the central C atom of the aliphatic backbone and the centre of mass of the adsorbed molecule (brown balls in Figure 3). Then, the molecule was moved in a perpendicular direction away from the cavity in steps (Figure 3). The whole rc length was about 40 Å. The length of the step varied, depending on a particular position of the molecule with respect to the HS cavity. In critical regions the step was 0.5 Å, but in other regions it was longer. All structural models were placed in a cubic supercell box with an edge length of 100 Å to avoid artificial interactions between periodic images when performing calculations with periodic boundary conditions. For each selected distance, the constrained MD simulation was performed using two phases. In the first phase of each MD run, equilibration was accomplished with a duration of~1 ns, followed by a production phase of an additional 1 ns. The equilibration phase was controlled by a temporal evolution of the potential energy. The production phase was used for the calculations of the PMF. All MD calculation steps were carried out using the canonical NVT ensemble (Andersen thermostat (Andersen, 1980) at T = 300 K). Newton's equations of motion were integrated using the velocity Verlet algorithm with a time-step of 1 fs. Interatomic interactions were described using the OPLS-AA force field (Jorgensen, Maxwell, & Tirado-Rives, 1996), whereas the TIP3P FF model (Jorgensen, Chandrasekhar, Madura, Impey, & Klein, 1983) was selected for water molecules. All FF-MD calculations were performed with the program TINKER (Ponder & Richards, 1987).

| RESULTS
In the calculations of free energy changes of chemical processes and also for comparative purposes it is necessary to define a reference state. In our simulations, we had complex systems consisting of several regions; that is, hydrophobic and hydrophilic domains represented by the aliphatic chains and carboxyl groups, and water molecular clusters representing microhydrated regions. Moreover, the studied molecules differ in their nature (i.e., MCPA is a polar organic acid molecule, whereas naphthalene is nonpolar moiety); thus, they behave differently in hydrophobic/hydrophilic surroundings. Therefore, the reference state selected corresponded to a configuration in which the molecules were placed far away from the HS-water model at rc distance of 40 Å from the aliphatic backbone of the HS model. The free energy value for the reference point was defined as zero for both molecular species; therefore, both energy profiles could be compared.
Another important aspect of the PMF calculations is a reliable prediction of the free energy changes. This method is frequently used in the simulation of very large and complex (bio)chemical systems and experience shows that reliable prediction of the free energy changes needs a very good configurational space sampling, which is not an easy task. It means that, at least, long MD simulations need to be performed or it is necessary to use some accelerated MD techniques (e.g., replica exchange). In any case, as the duration of the MD production phase was limited to 1 ns, it was important to validate the quality of the PMF simulations. To verify the consistency of the prediction of the FEP of the molecule trapped in the nanopore of the HS model, several time lengths of the production phase of the MD simulations were analyzed. Particularly, the PMF profiles were calculated for the last 1.0, 0.9 and 0.8 ns lengths of the production phase. The results obtained for the naphthalene molecule are shown in Figure S1. Comparison of the three PMF profiles demonstrated their similarity within the whole range of the reaction coordinates, with relatively small differences. All three curves have two distinct minima (details discussed later), differing by about 0.5 kcal/mol. Thus, it can be concluded that the 1 ns production phase was sufficiently long to provide a reliable description of the FEPs for the molecular trapping in the nanopore of the HS model. Similar trends were also observed for the MCPA-HS system (not shown in the SI). Figure 4 shows the FEP for the naphthalene molecule trapped in the HS nanopore model. The curve plot is accompanied by four structural snapshots (A-D) taken from the constrained FF-MD simulations. These four snapshots illustrate four distinct regions of the FEP curve. Two minima were clearly observed separated by a relatively small barrier (about 2.5 kcal/mol from the side of the global minimum, A and C snapshots in Figure 4). The first deepest minimum was about −5.0 kcal/mol at rc length of about 5 Å. The naphthalene molecule was trapped in the nanopore cavity, having a mostly parallel configuration with respect to the aliphatic backbone of the HS model. The aliphatic chains terminated by the carboxyl groups were flexible and formed multiple hydrogen bonds with water molecules from the wetting microregion of the HS model. This observation reflected the formation of the water molecular bridges (WaMB) between polar functionalities, as suggested by Schaumann and Bertmer (2008). Water molecules also formed a multiple hydrogen bond network themselves. The formation of the WaMB locked the naphthalene molecule in the HS nanocavity and contributed to the overall stability of the naphthalene adsorption. When the naphthalene molecule was pushed out of the cavity, the WaMB link was broken. This step required work to break hydrogen bonds, reflected by the appearance of an energy barrier on the FEP curve at a distance of about 10 Å. This barrier corresponds to the second snapshot (b) in Figure 4, which, moreover, also clearly demonstrated the hydrophobic nature of the naphthalene molecule. The molecule penetrated into the water cluster and water molecules formed a hydrophobic solvent cavity around the naphthalene, while also keeping hydrogen bonds with the carboxyl groups of the HS pore. The second shallow minimum (about −3.0 kcal/mol, (c) in Figure 4) was observed at a distance of~14-15 Å. In this position (water-HS interface), the naphthalene molecule was slightly stabilized by a weak solvent effect. Moreover, the terminal -CH 3 groups of the -CH chains also contribute to the naphthalene stabilization. In this configuration, the nanopore started to lock again by forming attractive hydrogen bond interactions between the carboxyl groups, which additionally contributed to the stabilization of the whole system. In the fourth domain, the naphthalene molecule travelled through the wetting microenvironment represented by the water cluster. The last snapshot in Figure 4 (d) shows the molecule in a region of the water-vacuum interface. The FEP curve approached zero, the limit of the reference state.

| MCPA-HS adsorption
In this simulation, we used the MCPA molecule in a neutral form. MCPA is a weak organic acid with pK a of~3.1. The neutral form of MCPA was used for two reasons: (a) to have a similar molecular state to naphthalene, which is a neutral molecule at normal conditions, and (b) several experiments revealed that the adsorption of MCPA increases with decreasing pH, showing highest adsorption at pH close to pK a (Ayar, Bilgin, & Atun, 2008;Kersten, Tunega, Georgieva, Vlasova, & Branscheid, 2014;Thorstensen, Olav, Ole, & Agnethe, 2001). Figure 5 displays the calculated FEP curve for the MCPA adsorption in the HS nanopore model, which shows several differences from the shape of the FEP curve of naphthalene (Figure 4). For MCPA only a distinct minimum was observed, and three domains could be resolved (A-C in Figure 5). In contrast to naphthalene, the position of the MCPA molecule inside the nanopore (snapshot A in Figure 5) did not correspond to a stable configuration. The relative energy difference with respect to the reference configuration for the MCPA trapped F I G U R E 4 Calculated FEP of naphthalene inserted in the HS nanopore model with water molecules included in its hydrophilic region. FEP: free energy profile; HS: humic substance [Color figure can be viewed at wileyonlinelibrary.com] in the nanopore was about 2.5 kcal/mol. The MCPA mostly interacted through nonbonding dispersion interactions with the hydrophobic -CH chains. When the MCPA molecule was brought out of the nanopore hole, the whole system became more stable. The molecule broke the molecular water bridge between the carboxyl groups, which were involved in the formation of the stable hydrogen bonds with the carboxyl and -C-O-C-functionalities of the MCPA. All snapshots in Figure 5 document the flexibility of the aliphatic chains terminated by the carboxyl groups, which contributed to the effective formation of the hydrogen bonds. The hydrogen bond lengths were in a range of 1.7-2.2 Å, corresponding to hydrogen bonds of medium strength. The middle snapshot in Figure 5 (B) refers to the stable MCPA configuration, with a distance of about 10-12 Å between the molecule and the reference point of the aliphatic HS backbone. This snapshot also shows formation of a hydrophobic cavity of water molecules around the aromatic part of the MCPA molecule. The estimated minimum of the stable configuration of the MCPA adsorption was about −4.5 kcal/mol, being a similar value to that for the naphthalene trapping inside of the nanopore hole (−5.0 kcal/mol). This value is not far from the DFT-calculated free energies of the exchange reactions of the MCPA molecule interacting with several polyacrylic acid fragments in a polar solvent (−4.5 to −5.9 kcal/mol) (Aquino et al., 2007). Further motion of the MCPA molecule through the water microcluster stabilized it by about −2.5 kcal/mol with respect to the reference state (snapshot C in Figure 5). This stabilization was mainly due to the solvation effect of water molecules, even though any distinct minimum is visible at the FEP curve. This effect includes hydrogen bond interactions of water molecules with the carboxyl and C-O-C groups of the MCPA and weak dispersion interactions with the aromatic ring. In addition, in this configuration, the hydrogen bonds between the carboxyl groups of the HS model and water molecules were reconstructed as well.

| DISCUSSION
The calculated PMF curves of the naphthalene and MCPA molecules have similarities and/or differences, as discussed in the following. Firstly, the PMF depths were comparable (−5.0 vs. −4.5 kcal/mol), documenting that adsorption of polar and nonpolar species in SOM can be of similar strength. In our case, the interaction of naphthalene with the HS pore model was slightly stronger than the interaction of F I G U R E 5 Calculated free energy profile (FEP) of MCPA inserted in the HS nanopore model with water molecules included in its hydrophilic region. FEP: free energy profile; HS: humic substance; MCPA: 2-methyl-4-chlorophenoxyacetic acid [Color figure can be viewed at wileyonlinelibrary.com] MCPA. A similar trend was observed in the work by Ahmed et al. (2015). The authors performed FF-MD simulations to analyze the interactions of nonpolar hexachlorobenzene (HCB) and polar sulphanilamide (SAA) molecules with several models of HSs containing hydrophobic and/or hydrophilic domains. They found that the interactions of HCB with all HS models were always stronger than the interactions with SAA.
The main difference in the interaction of naphthalene and MCPA with the HS nanopore model was the form of binding and the domain in which the molecule was stabilized. In the case of MCPA, the main binding mechanism was represented by hydrogen bonding of the MCPA carboxyl group with the carboxyl groups of the HS model and with water molecules forming wetting microregions in the HS moieties. As expected, polar molecules preferred binding to the hydrophilic domains. On the other hand, the naphthalene molecule was encapsulated inside of the nanopore of the HS model, which is mainly formed by nonpolar aliphatic chains. Thus, the dominant mechanism has its origin in nonbonding dispersion interactions. Generally, nonpolar species will prefer binding in hydrophobic domains and micro(nano) cavities of HSs. Similar observations were recorded by Ahmed et al. (2015).
Our observations also indicated that kinetics of adsorption can be different between polar and nonpolar species. Generally, SOM is structurally very complex organic matter, having a form of supramolecular aggregates with varying density. Denser sectors contain cross-linking domains with concentrated polar functional groups, water and counter ions. These sectors can alternate with less dense regions, mainly represented by nonpolar hydrophilic functionalities, nanovoids and nanopores. Therefore, it is difficult to determine a distinct surface of SOM. It also means that adsorption sites in real SOM are not only at the surface but also inside the structures, differing by chemical affinity. Thus, small molecules cannot only be adsorbed at the surface but also penetrate (diffuse) inside the SOM structure. Therefore, it is challenging to develop a general model for the adsorption of organic species in SOM sorbents.
One of the proposed concepts for adsorption of organic species in SOM is the dual-mode model combining the dissolution mechanism in SOM with pore/hole-filling adsorption (Chefetz & Xing, 2009;Pignatello & Xing, 1996). The simulations carried out in this work showed that polar organic species could be (relatively rapidly) adsorbed directly at outer water-SOM interfaces, which can be formed by polar functional groups representing hydrophilic regions. Therefore, further penetration of polar species deeper inside the SOM structure can be significantly hindered or at least reduced. On the other hand, nonpolar species would prefer sorption inside the SOM structures, in hydrophobic domains present inside the structure in a form of nanopore void. However, some (faster) adsorption can also occur at the water-SOM interface, as indicated by a less deep minimum at the PMF curve of naphthalene in Figure 4. The process of penetration into the SOM can be relatively slow as a barrier has to be also overcome. This process can be described as a pore space opening during which hydrogen bonds from the WaMB have to be broken (Aquino et al., 2009;Schaumann & Bertmer, 2008). After transporting molecules inside nanovoids the WaMBs can be reconstructed again. Consequently, nonpolar species can be trapped inside the SOM structures and their eventual release (desorption) can be more difficult than for polar species at the water-SOM interfaces.
Although numerous studies exist on sorption of polar and nonpolar species in soils, SOM and/or pure humic/fulvic acids, it is hard to compare our results directly with other experiments. Most of the outputs from adsorption experiments are not adsorption energies (enthalpies) but distribution (adsorption) coefficients, such as K D or K OC; see, for example, several comprehensive works collecting data for various organic species, including Schwarzenbach (2006a), Niederer, Goss, andSchwarzenbach (2006b), Kah and Brown (2006), Franco & Trapp (2008), Shao et al. (2014) and Wang et al. (2015). The problem is deriving adsorption energies, which require determination of equilibrium constants at several temperatures and fitting the results using the van't Hoff equation. Another difficulty hampering direct comparisons is the fact that adsorption experiments are mostly conducted in a solution phase, whereas in this study the reference state is defined as an isolated molecule in a gas phase. In addition, experimentally determined adsorption parameters can differ considerably because of the chemical and structural variation of SOM sorbents (Grathwohl, 1990;Niederer, Schwarzenbach, & Goss, 2007;Xing, 1997). For example, the data collected by Niederer et al. (2007) for logK OC /logK DOC values of naphthalene vary in the range 2.69-3.74 L/kg. Nevertheless, experimental values of the K OC constants for naphthalene are generally higher than those for MCPA (Franco and Trapp (2008) provided logK OC values of 1.82 and 1.87 L/kg for MCPA). This documents that the naphthalene adsorption in SOM is higher than for MCPA, which is in good agreement with our observations.
There are only a few papers with experimentally determined energies (enthalpies) of MCPA/naphthalene adsorption on SOM in the literature. Adsorption free energy of MCPA on bituminous shale (model adsorbent containing solid organic matter in a mineral matrix) was determined in a range of −2.4 to −3.2 kcal/mol from fitting to four different temperature measurements (Ayar et al., 2008). These values agree quite well with our estimated value of about −2.0 kcal/mol, related to the difference between the minimum and the PMF value at about 20 Å distance, corresponding to MCPA solvated by water molecules (snapshot C in Figure 5). For naphthalene, Niederer et al. (2006a) and Niederer et al. (2006b) provided enthalpy of adsorption to Leonardite humic acid from the gas phase of −13.7 kcal/mol. The authors did not specify the entropy or free energy values. Using the TΔS value of~5.5 kcal/mol from the adsorption experiment of MCPA on bituminous shale (Ayar et al., 2008) reveals an estimated adsorption free energy for naphthalene of about −8.8 kcal/mol. This value is larger than our calculated minimum from the PMF calculations (−5.0 kcal/mol) but the difference is acceptable if all factors influencing adsorption on SOM are taken into consideration.

| CONCLUSIONS
In this work, classical molecular dynamic simulations were performed to obtain FEPs of hydrophilic/hydrophobic interactions between polar and nonpolar molecules encapsulated in the nanopore of the HS model, which includes both polar and nonpolar domains. The calculated PMF curves of the naphthalene and MCPA have comparable depths of −5.0 and − 4.5 kcal/mol, respectively, supporting that the adsorption of polar and nonpolar species in SOM can be of similar strength. These similar stabilities are explained by the fact that the polar MCPA interacts with the carboxyl groups of the HS model by hydrogen bonding, whereas the naphthalene is stabilized through nonbonding dispersion interactions with aliphatic chains representing hydrophobic domains of SOM.
Regarding the kinetics in the adsorption process, we can infer from our results that polar organic species could rather be adsorbed directly at the outer water-SOM interfaces, which contain polar functional groups. This can prevent or reduce additional diffusion of polar species deeper inside the SOM structure. For nonpolar species like PAHs, quite a rapid adsorption at the water-SOM interface can also be expected, as deduced by the second (less deep) minimum of the PMF curve of the naphthalene molecule. However, diffusion into the SOM is expected to be relatively slow as this process would involve a pore space opening in which the protecting hydrogen bonds of the WaMB must be broken and, thus, an energy barrier be overcome. Therefore, also desorption of nonpolar species from inside regions of the SOM can be more difficult than for the polar species at the water-SOM interfaces. The presented work documents the capability of molecular dynamic simulations to provide detailed molecular-scale information about the mechanism and nature of adsorption processes of organic pollutants, such as herbicides in SOM, giving us better understanding of their behaviour and fate in soil environments. We believe that our work can be helpful for further experimental studies focusing on revealing adsorption mechanisms of various organic contaminants in SOM.

DATA AVAILABILITY STATEMENT
The data that support the findings of this study are available from the corresponding author upon reasonable request.