Application of Density Functional Theory in Soil Science

Abstract

Soil is the basis for life and soil science is regarded as the final frontier; however, as compared to chemistry, physics, biology, and other disciplines, soil science undergoes an obviously slower development and remains almost stagnant in the past few decades, mainly due to two reasons: (1) wrong and outdated perceptions for a large portion of soil researchers; (2) complexity of soil systems that are difficult to characterize by current experimental techniques. Computer simulations have unique advantages to handle complex systems while currently, its role during soil researches is far from being recognized. In this chapter, several examples are given with respect to application of density functional theory (DFT) calculations to soil science, focusing on the adsorption of uranyl ion and SO2 onto mineral surfaces and reaction mechanisms to form acid rain. In this way, insightful clues at the atomic level are provided for the adsorption, interaction, and reactions regarding soil systems. We believe that computer simulations including DFT are the right key to unravel the complicated processes occurring in soils. More efforts of computer simulations are anticipated for soil science with aim to decipher the experimental results and probe the uncharted principles that may result in a revolutionary in the near future.

1. Introduction

According to Natural Resources Conservation Service (NRCS), soil is defined as a natural body comprised of solids (mainly minerals and organic matters), liquids, and gases that occurs at the intermediate surface of the Earth, occupies space and is characterized by one or both of the following properties: horizons and layers, which are distinguishable from the initial materials as a result of addition, loss, transfer, and transformation of energy and matter or the ability to support rooted plants in natural circumstances. Soil constitutes the basis for life and bridges the biosphere, atmosphere, hydrosphere, and geosphere. Despite these facts, apparently less attention has been given to soil science than to other disciplines such as physics, chemistry, and biology. As said by Gardner (the past president of Soil Science Society of America) [1], “not a few people mistakenly perceive that everything worth knowing about soils has already been understood, and all we need to do is merely to apply that knowledge properly. Even knowledgeable scientists assume that principles and theories developed from other disciplines can be applied to the researches of soil science in a straightforward way, requiring little imagination or creativity. In their opinions, soil science is just one of expressions for the applied physics, chemistry, or biology.”

The situation of soil science research is alarming. The core concepts of current soil science textbooks remain almost unchanged as compared to those of half a century ago. Obsolete or even incorrect standpoints are a commonplace [2]. Fortunately, a few researchers have recognized such a crisis. On the other hand, because of the complexity of systems and co-function of multiple factors, it seems challenging for experimental techniques to in-situ characterize “real” soil properties and processes; in addition, the experimental results from one lab may not be reproducible by others, since soil samples of different areas or even different batches may vary significantly. Computer simulations have unique advantages within this context: (1) probing the various influencing factors one by one; e.g., six factors (identity of heteroatoms, crystallographically distinct T sites, structural alterations, quantity of negative charges, distance from charge centers to metal ions and source of negative charges) were identified to affect the adsorption of metal ions at clay surfaces, and their respective contributions were estimated by density functional theory (DFT) calculations. The quantity of negative charges is the foremost factor that controls the adsorption processes, while other factors in certain circumstances can also play a critical role [3]. The adsorption strengths and numbers of all metal ions increase in a direct proportion to the intensities of electric fields [4]; (2) providing useful and detailed information at the femtosecond scale such as how ions from aqueous solutions diffuse toward to clay surfaces [5]; (3) understanding the adsorption, interaction and reaction processes at the atomic level such as how metal ions interact with surface-O atoms and respond to the increase of electric fields. Based on Hirshfeld, Mulliken, and NBO charge analyses, we found that polarization rather than electrostatic interactions are more likely to result in the pronounced cation-specific effects at clay surfaces [6]; (4) unraveling the exact reaction mechanism by comparing the structural and (especially) activation barriers of competing paths. This can be considered as an extension of (3). There are a plethora of competing reactions occurring in soils; e.g., with DFT calculations, it was clarified that Mn⁴⁺ rather than Mn³⁺ sites are more reactive for the oxidation of As³⁺ and the oxidation processes are significantly blocked by As⁵⁺ complexes [7]. A more convictive example is the mechanistic study of Brönsted acid-catalyzed conversion of biomass sugars [8]. More than 120 reaction paths were explored, and the low reactivity and selectivity of glucose conversions were clearly addressed: unlike fructose that prefers to dehydrate at the anomeric O₂H group and initiates a sequence of facile reaction steps toward 5-hydroxymethyl-2-furfural (HMF), the less reactive sites in glucose (O₂H and O₃H) produces levulic acid not involving fructose and HMF intermediates, while the most reactive O₁H site leads to humin precursors or reversion products [8].

The relationship between system size and computational accuracy for representative theoretical levels is shown in Figure 1. For huge systems (100,000 or even more atoms), classical methods (Monte Carlo and Molecular Dynamics) seems to be the good choice, although the computational accuracy is relatively low; In contrast, ab initio quantum mechanical methods are restricted to relatively small systems (up to several hundred atoms or even fewer atoms) while achieve high accuracy. The computational accuracy of quantum mechanics/molecular mechanics (QM/MM) [9] and semi-empirical methods fall in-between, and QM/MM methods have recently become increasingly popular due to the satisfying computational accuracy (active sites handled by QM methods) and the easy extension to large systems (all atoms except active sites disposed by classical methods). 2013 Noble Prize in Chemistry was awarded to Karplus, Levitt, and Warshel for the development of “multiscale methods for complex systems.”

Figure 1.

The theoretical landscape regarding to the balance between system size and computational accuracy.

Let us look back to ab initio quantum mechanical methods. Strictly speaking, a large portion of DFT methods fall outside this scope; e.g., B3LYP is probably the most popular DFT method and its exchange-correlation functional is written as [10, 11]

$$E_{\mathit{xc}}^{B3\mathit{LYP}}=\left(1-a\right)E_x^{\mathit{\text{LSDA}}}+{\mathit{aE}}_x^{\mathit{HF}}+{\mathit{bE}}_x^{B88}+{\mathit{cE}}_c^{\mathit{LYP}}+\left(1-c\right)E_c^{\mathit{\text{LSDA}}}$$

where a, b, and c are empirical parameters.

There is a five-rung Jacob’s ladder of common density functional approximations, as local spin-density approximation (LSDA), generalized-gradient approximation (GGA), meta-GGA, hyper-GGA, and random phase approximation (PRA-like functionals), see more details in [12]. As far as we know, the most widely density functionals are at GGA (e.g., B3LYP, PBE) and meta-GGA (e.g., M06 L, M06-2X) levels.

Computer simulations, as discussed above, are particularly useful to handle the complex soil systems, and two applications were illustrated in this chapter. Section 2 summarized DFT calculated results of uranyl adsorption onto mineral surfaces and Section 3 elaborately discussed DFT calculated results of SO₂ (one of the main precursors for acid rain) adsorption and formation mechanisms of acid rain. In the end, concluding remarks were provided and some suggestions for DFT applications to the future soil researches were presented as well. A more critical role of computer simulations should have played in soil science, and this chapter aims to arouse the attention of general soil researchers regarding to the applications of computer simulations. More efforts in this regard are greatly beneficial to decipher the experimental results and probe the uncharted principles that will result in a revolutionary for soil science in the near future.

2. Uranyl adsorption onto mineral surfaces

Release of radionuclides into the environment seriously threatens the ecosystem and human health, and adsorption of radionuclides onto mineral surfaces significantly affects their migration and transport into the environment [13]. Accordingly, knowledge about the interaction between radionuclides and minerals is essential for the long-term risk assessment of radioactive waste repositories. Uranium is usually present in the uranyl (UO₂²⁺) form [14, 15] and DFT calculations have been conducted to understand the adsorption of UO₂²⁺ onto mineral surfaces [16, 17, 18, 19, 20, 21, 22, 23, 24, 25].

Perron et al. [16] examined UO₂²⁺ adsorption onto rutile(110) surface, the most stable face of natural rutile reported by Jones and Hockey [26, 27]. UO₂²⁺ forms a bidentate inner sphere complex with three H₂O molecules to fill its first hydration shell, see Figure 2. There exist two types of surface-O atoms, as O_t (terminal-O) and O_b (bridging-O) that are singly and doubly coordinated, respectively. As a result, a total of three adsorption structures are produced: UO₂²⁺ coordinated to two O_b atoms (bb), two O_t atoms (tt) and both of O_b and O_t atoms (bt). The bb mode is the most preferred and has a lower energy than the bt and tt modes by 5.0 and 13.6 kJ/mol, respectively.

Figure 2.

Adsorption structures for uranyl ions on rutile(110) surface: tt (Left), bb (Middle), and bt (Right). U, Ti, O, and H are respectively in yellow, blue, red, and white, and H atoms to saturate surface-O atoms are not shown for clarity.

In addition to the [UO₂(H₂O)_n]²⁺ adsorption complexes, other uranyl species containing anionic ligands such as OH⁻ and CO₃²⁻ also play a critical role in the environmental circumstances. Pan and co-workers [17] investigated the effects of different ligands (H₂O, OH⁻, CO₃²⁻) on the adsorption on rutile(110) surface. Note that only the most stable adsorption mode (bb) has been taken into account therein. The uranyl ion (UO₂²⁺) interacts more strongly with anionic ligands (OH⁻ and CO₃²⁻) than H₂O in the first coordination sphere. The bond lengths of the U⁶⁺ center and two O_b atoms are respectively calculated at 2.248 and 2.333 Å for [UO₂(H₂O)₃]²⁺, 2.808 and 2.930 Å for [UO₂(OH)₃]⁻, 2.852 and 2.960 Å for [UO₂(CO₃)₂]²⁻, suggesting that the adsorption strengths of the uranyl ion with rutile(110) surface are significantly impaired in the latter two cases. Adsorption of these uranyl complexes on partially hydrated surfaces were then studied, which exhibit stronger interactions and are attributed to the formation of H-bonds between H₂O molecules and uranyl ion (UO₂²⁺) that strengthen the adsorption processes.

Kaolinite is the major constituent of sedimentary clay rocks, and adsorption of uranyl species onto its (001) surface [20, 21] and (010) edge surface [22, 23] have been theoretically studied. The adsorption of uranyl species on tetrahedral SiO₄ surfaces is thermodynamically unfavorable, with the adsorption energies of inner- and outer-sphere complexes being, respectively, 239 and 206 kJ/mol [20]; in contrast, the adsorption at octahedral AlO₆ surfaces is preferred due to the presence of upright (perpendicular to the surface) and lying (parallel to the surface) OH groups. The most stable configuration corresponds to an inner-sphere monodentate complex, where UO₂²⁺ bonds directly to the lying surface-OH group and has an adsorption energy of −155 kJ/mol. Martorell et al. [21] continued to study the adsorption of uranyl ion on the bare and solvated octahedral AlO₆ surfaces of kaolinite, where UO₂²⁺ forms two direct bonds with the surface-O atoms of deprotonated OH groups. Two different adsorption structures are assigned according to the coordination of surface-O atoms: two surface-O atoms connected to one and two Al atoms are designated to be AlOO (short-bridge site) and AlO-AlO (long-bridge site), respectively. The formation of adsorption complex at the short-bridge site (AlOO) requires an energy of 195 kJ/mol and is obviously less than at the long-bridge site (AlO-AlO, 261 kJ/mol). Thus, the uranyl ion prefers to adsorption at the short-bridge site; furthermore, similar trends for uranyl adsorption remain when including solvent effects by adding a monolayer of water molecules.

Gibbsite is a primary mineral form of aluminum hydroxide, and adsorption of the uranyl ion onto its (001) surface [18, 24] and edge surface [25] have been investigated, similar to the situation of kaolinite discussed above. Figure 3a depicts three potential adsorption sites for uranyl adsorption: Sites I and III hold almost the same O-O distances (about 2.7 Å), while site II shows an apparently elongated O-O distance (about 3.4 Å). Site I represents the lowest-energy adsorption structure (Figure 3), and the adsorption structures of sites II and III are less stable with relative energies being 22.2 and 46.3 kJ/mol, respectively. Three H-bonds constructed between surface-OH groups and UO₂²⁺ center as shown in Figure 3 stabilize the interactions between uranyl ion and gibbsite surfaces.

Figure 3.

(a) Three sites for uranyl adsorption on gibbsite(001) surface and (b) top and (c) side views of the optimal adsorption complex. U, Al, O, and H are, respectively, in blue, pink, red, and white, and H-bonds are in black dash lines and water molecules are in circles.

3. Acid rain

In 1852, Smith demonstrated the relationship between acid rain and atmospheric pollution in Manchester, and after 20 years (i.e., 1872), he coined the term “acid rain.” Now acid rain has become a popular term and one of the world’s biggest environmental concerns, especially in North America, Europe, and China [28]. Acid rain refers to any form of precipitations with acidic components that fall to the ground from the atmosphere, including rain, snow, flog, hail, or even dust that is acidic. It results mainly from SO₂ and NO_x emissions to the atmosphere and the further transport by wind and air currents, during when SO₂ and NO_x react with water, oxygen, and other substances leading to the formation of sulfuric (H₂SO₄) and nitric (HNO₃) acids. The pH of acid rain is approximately 5.6 and since 1940s, researchers began to recognize its strong impacts on the ecosystem and human health so that soils, freshwaters, forests, and buildings will be damaged. With regard to soils, acid rain inhibits the decomposition of organic matter [29], fixation of nitrogen [30], elution of calcium (Ca²⁺), magnesium (Mg²⁺), potassium (K⁺), and other nutrients [31]. As a result, soil fertility and microbial activity show an obvious reduction [32]. Geochemical modeling indicated that Ca²⁺ leaching in marble due to acid rain neutralization approximates 0.158 mmol/L, in contrast to 10.5 mmol/L by dry deposition, and the corresponding Cu²⁺ losses in bronze are ca. 0.21 and 47.3 mmol/L, respectively [33].

As aforementioned, SO₂ emissions are one of the principal causes of acid rain and also represent a primary source of atmospheric aerosols, which can lead to respiratory diseases, premature deaths, and even climate changes by affecting the properties of clouds and the balance of solar radiation. Therefore, it is of great significance to convert SO₂ to other less contaminated compounds, and a number of measures to control SO₂ emissions have been proposed. During 1982–1999, SO₂ emissions have reduced by approximately 65% in Europe and 40% in the United States, and SO₂ emissions in China decline in the late 1990s while again increase after then. DFT calculations provide useful information about the adsorption of SO₂ onto mineral surfaces as well as reaction mechanisms that seem difficult to capture by current experimental techniques [34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45], which are, however, critical to understand the formation of acid rain at the molecular level and to remediate the ecosystem. Clay minerals, such as alumina (Al₂O₃), iron oxides (Fe_xO_y), are good candidates for the adsorption of acid components from acid rain and then convert them into less hazardous compounds. Lo et al. [35] studied the adsorption of SO₂ on clean (100), dehydrated (110), and hydrated (110) surfaces of γ-Al₂O₃, finding that significant adsorption differences exist for the various surfaces and the calculated adsorption energies (−13 to –85 kcal/mol) are consistent with experimental results.

The γ-Al₂O₃(100) surface is composed by bridging-O and five-coordinated Al atoms, and a total of five stable configurations are produced for SO₂ adsorption (Figure 4). The feeble interaction between S and surface-O atoms results in the physisorption configuration (CM3) with the S-O distance of 2.915 Å, and the corresponding binding energy is very small (−2.0 kcal/mol). The interaction between O@SO₂ and Al atoms leads to a chemisorption state named CM4, and the O-Al bond distance and adsorption energy are 2.123 Å and −23.9 kcal/mol, respectively. The other three configurations are also ascribed to chemisorption. In CM5, one O@SO₂ atom is coordinated to two Al atoms in the vicinity of an octahedral vacancy, and CM7 and CM8 can be considered to generate from CM5 conversion and recombination. As compared to CM5 and CM7, the adsorption configuration CM8 possesses a superior symmetry, and both O@SO₂ atoms participate in the formation of direct bonds with the Al atoms. The adsorption energies of SO₂ are calculated to be −39.2, −15.1, and −45.2 kcal/mol, respectively, for CM5, CM7, and CM8. In consequence, three types of sulfite (SO₃²⁻) are produced during the adsorption of SO₂ onto γ-Al₂O₃(100) surface. For all adsorption configurations including three with positive adsorption energies (CM1: 1.0 kcal/mol, CM2: 2.7 kcal/mol, and CM6: 21.4 kcal/mol), no direct coupling is detected between S@SO₂ and Al atoms. When γ-Al₂O₃(110) surface is hydrated, five stable adsorption configurations arise that are distinct from dehydrated condition: two physisorption modes (HM1 and HM2) and three chemisorption modes (HM3, HM4, and HM5), see Figure 4. HM1 and HM2 are structurally similar in that their S atoms are coordinated to a surface hydroxyl, while the coordination numbers of their Al atoms are different from each other. HM3 is produced by interaction of O@SO₂ atom with five-fold Al sites. HM4 and HM5 contain the sulfite species where the S atom is coordinated with surface-O atoms. The adsorption energies are calculated to be −20.4, −25.3, −31.1, −17.5, and −35.0 kcal/mol, respectively, for HM1, HM2, HM3, HM4, and HM5. In consequence, HM5 with formation of the sulfite species represents the lowest-energy adsorption configuration, which is the same as in dehydrated condition (CM8). Two IR peaks at 1214 and 1349 cm⁻¹ are assigned to the sulfate species, which can be been finely interpreted by DFT calculated results.

Figure 4.

Optimized configurations of SO₂ adsorption on (a) dehydrated (100) and (b) hydrated (110) surfaces of γ-Al₂O₃, where Al, H, O, S, and O attached to S are presented in green, pink, white, blue, and red balls, respectively.

Goethite (α-FeO(OH)), which can be found in soils and other low-temperature environments, is an iron-bearing hydroxide. Because of the considerable adsorption capacity for organic acids and anions, goethite has also been widely used in environmental remediation and protection [35]. Zubieta et al. [37] investigated the adsorption of SO₂ on partially and fully hydrated (110) surfaces of goethite and obtained eight stable products: six sulfite, one bisulfate, and one sulfate (Figure 5). The six adsorption structures containing sulfite species, created only on two types of partially dehydrated goethite surfaces, are further divided into two monodentric mononuclear (MdMn) and four bidentate (Bd) configurations. In the MdMn configurations (I and II), the S-O_Fe distances are elongated as compared to the other S-O distances, and they display two symmetrical stretching modes (OSO and OSO_Fe) centered at ca. 1126 and 976 cm⁻¹. In BdBn configurations (I and II), the two S-O_Fe distances are approximately 1.62 Å, and although with similar geometries and stretching modes, the vibrational frequencies deviate significantly from those of MdMn configurations and fall at around 672 and 661 cm⁻¹. In BdPn configurations (I and II), one S-O_Fe distance is optimized at 1.75 Å and lengthened as compared to those of BdPn configurations, while the other S-O_Fe distance equals 1.55 Å and is obviously contracted. The Bader analyses indicate that all sulfite species carry approximately −1.4 |e| charges.

Figure 5.

Optimized configurations for the adsorption of SO₂ on partially hydrated goethite(110) surfaces forming the sulfite, bisulfite and sulfate species (violet = Fe, red = O, white = H, yellow = S).

The bisulfate species is formed when the S atom constructs direct bonds with μ₁-OH (Figure 5), which is further stabilized by two H-bonds (1.77 and 1.81 Å). A similar bisulfate species was observed by Liu et al. [35] during the adsorption of SO₂ on γ-Al₂O₃ surfaces. The Bader analyses indicate that the HSO₃ species (bisulfate) carries −0.8 |e| charges. The computational sulfite and bisulfate species have IR spectra that are consistent with experimental results [37, 38]. The sulfate species emerges only when two water molecules are deprived from the (110) surface of goethite. In this way, two O atoms are singly coordinated and one O atom is three-fold coordinated with the Fe ions of goethite(110) surface. The sulfate species is corroborated by one H-bond (1.81 Å) and carries −1.6 |e| charges according to the Bader charge analyses. A stable sulfite structure can be produced under identical conditions, whereas its adsorption energy is obviously less, indicating that the formation of the sulfate species is significantly preferred.

The adsorption of SO₂ onto the Cu(100), MgO and carbon surfaces was discussed as well, which may provide insightful clues for resembling processes onto mineral surfaces. It was proposed that SO₂ and H₂O are co-adsorbed onto Cu(100) surfaces [30], through the direct coupling of Cu atom with S and O@H₂O atoms. The Cu-O distances ascend in the order of co-adsorption of SO₂ and H₂O < adsorption of only SO₂ < < adsorption of only H₂O. Accordingly, the interaction between SO₂ and Cu(100) surface is stronger than that of H₂O, and the co-adsorption of SO₂ conduces to the enhanced interaction of H₂O with Cu(100) surface. At the same time, the Cu-S distance of the co-adsorption configuration is optimized at 2.385 Å and is shorter than that with only SO₂ adsorption. That is, water exhibits a promoting effect for the adsorption of SO₂ on Cu(100) surface, as corroborated by the calculated adsorption energies. Eid and collaborators [41] found that as compared to regular MgO surface, the adsorption capacity of SO₂ at MgO(Fs-center) defects is higher, and MgO(Fs-center) corresponds to an enhanced catalytic activity. With regard to pure carbon materials, SO₂ is physisorbed and van der Waals (vdW) is the driving force therein [42]. When carbon materials are modified with functional groups such as carboxyl, lactone, or/and phenolic hydroxyl, the adsorption strengths of SO₂ are enhanced pronouncedly, especially for the sites at edge surfaces. In addition, these functional groups show little effects on SO₂ adsorption, suggesting that the enhanced adsorption is mainly due to regulation of carbon surface properties.

Adsorption of SO₂ is the first step for the formation of acid rain. According to our preliminary studies, the reaction mechanisms of gas phase and mineral surfaces resemble each other, and hence the gas-phase results are beneficial to understand acid rain formation at mineral surfaces. The reaction of SO₂ and H₂O produces two isomers that have close electronic energies [43]; in addition, the two isomers have apparently lower electronic energies than H₂SO₃. Accordingly, the gaseous SO₂ and H₂O mixture is likely to exist as the SO₂·H₂O complex rather than H₂SO₃. The activation barrier of SO₂ reacting with H₂O to form H₂SO₃ is so high (146.7 kJ/mol) that it becomes very difficult to produce the sulfite species (H₂SO₃) in gas phase. Five years later, Stirling [44] investigated the hydrolysis of SO₂ in aqueous solutions, finding that hydrated SO₂ forms the bisulfite anion and hydronium ion after overcoming an energy barrier of about 17 kcal/mol, while the one-step formation of H₂SO₃ has not been detected. The orientation of water molecules in the hydration shell of SO₂ implies a more facile formation of the bisulfite anion rather than H₂SO₃ [45], in line with the results of meta-dynamics calculations [46]. When HO₂ participates in the reaction, the S atom constructs a new bond with O@H₂O atom (S-O: 1.716 Å). The energy barrier reduces considerably and equals 56.6 kJ/mol. HO₂ exists widely in the atmosphere and participates in a variety of chemical reactions [47]. In addition to HO₂, a number of computational studies explicitly indicated that acidic substances in the atmosphere play similar catalytic effects and reduce significantly the energy barriers for the hydrolysis of SO₂ [48, 49, 50]. The effects of H₂SO₃ and H₂O on the hydrolysis of SO₂ were studied by Liu et al. [48], showing that the energy barriers ascend in the order of H₂SO₃ + SO₂·H₂O (26.4 kJ/mol) < H₂O + SO₂·H₂O (84.0 kJ/mol) < SO₂·H₂O (154.6 kJ/mol). The catalytic effect of H₂SO₃ is obviously more pronounced than that of H₂O. As reflected by NBO analyses, the reactant complexes in presence of sulfurous acid have the enhanced second-order stabilization energies as compared to those with addition of only water molecules, which provides strong supports for the pronounced catalytic effect of H₂SO₃. The kinetic calculations further show that the hydrolysis of SO₂ is a nearly autocatalytic reaction.

Alternatively, acid rain can be formed by release of SO₂ into the atmosphere and oxidization of SO₃, which then reacts with H₂O to form sulfuric acid (H₂SO₄). There are many oxidants that are able to convert SO₂ to SO₃, and Calvert [51] have provided the formation mechanism of sulfuric acid (H₂SO₄)

SO₂ + OH → HSO₃

HSO₃ + O₃ → SO₃ + HO₂

SO₃ + H₂O → H₂SO₄

Interaction between SO₂ and the OH radical (OH^•) has been addressed by ab initio electronic structure calculations [52], and the optimized structure of the HOSO₂ radical show differences with that predicted using a small basis set (MP4/6-31G**//HF/3-21G*) [53]: Inclusion of correlation effects results in an increase of S-O distances and OSO angle, whereas the HOS angle shows a substantial decrease. Although geometrically stable at all levels of theory, stability of the HOSO₂ radical is estimated to be 104–110 kJ/mol, suggesting that the direct dissociation to SO₃ is almost infeasible. O₃ and H₂O₂ are two oxidizing substances in atmosphere chemistry, and their reaction mechanisms with SO₂ were studied by Jiang et al. [54].

SO₂ + O₃ → SO₃ + O₂

SO₂ + H₂O₂ → H₂SO₄

SO₂ and O₃ interact mainly through the S and O@O₃ atoms, and their reaction causes the formation of S-O bond (1.708 Å) and the rupture of O-O bond in O₃. This is one-step process and requires to overcome an energy barrier of 46.5 kJ/mol. The reaction of SO₂ with H₂O₂ proceeds via the OH-abstraction mechanism, and two OH radicals generated from the dissociation of H₂O₂ are appended to SO₂ forming H₂SO₄. However, the energy barrier of this reaction is 299.5 kJ/mol, which is extremely difficult to proceed at normal conditions. Chen et al. [46] reported the reaction mechanism between SO₂ and HO₂ and showed that there exist two types of SO₂·HO₂ complexes: one is to combine the terminal O@HO₂ and S atoms. In this reaction, the O-S distance decreases from 2.966 Å (the initial complex: SO₂·HO₂) to 1.598 Å while the O-O bond of HO₂ is gradually elongated until broken. The reaction of SO₂ and H₂O is divided into three stages, and the energy barrier of SO₃·OH formation is so small (1.7 kJ/mol) that can be considered negligible. The final product of this reaction is HSO₄. The other complex is characterized by two types of intermolecular H-bonds that form between the terminal O@HO₂ and S atoms and between the H@HO₂ and O@SO₂ atoms. It is a one-step process and the product for this reaction is HOSO·O₂. When a single water is added, it has a combined effect on the overall process, not only accelerating the reaction of generating HSO₄ by reducing the activation barrier of the second step from 48.8 to 44.8 kJ/mol but also inhibiting the reaction between SO₂ and HO₂ to produce HOSO due to blocking the interaction between the H@HO₂ and O@SO₂ atoms.

It can be seen from the above discussions that the SO₃ and water reaction is an integral section for the formation of acid rain. The direct reaction of SO₃ with one water to produce sulfuric acid (H₂SO₄) requires to overcome a large energy barrier (28.7 kcal/mol) and seems difficult to occur at normal conditions [55, 56, 57]. When the second water participates, the energy barrier reduces substantially and the reaction becomes almost barrierless, which is probably due to the formation of a stable six-membered cyclic transition state, see Figure 6 [58]. The two water molecules transfer their protons in a concerted manner, and a new S-O bond is formed between the newly generated OH and SO₃ fragments. In consequence, multiple reaction pathways may co-exist as illustrated below

Figure 6.

Optimized transition states for the formation of sulfuric acid from SO₃ in presence of one (left) and two (right) water molecules. Arrows indicate the reaction coordinate vectors, and some important distances (Å) and angles (degrees) are shown.

SO₃ + (H₂O)₂ → H₂SO₄ + H₂O

SO₃·H₂O + H₂O → H₂SO₄ + H₂O

SO₃ + (H₂O)₂ → H₂SO₄ + H₂O

The energy barriers for the reaction of SO₃ and water clusters (H₂O)_n to form H₂SO₄·(H₂O)_n-1 are calculated to be 28.7, 11.1, and 4.6 kcal/mol for n = 1, 2, 3, respectively [59]. With increase of water numbers, more water molecules are available to solvate and stabilize the charged transition state complexes, which further reduce the energy barriers. That is, it is favorable to hydrolyze SO₃ to form sulfuric acid in presence of sufficient water vapor. The second molecule acts as a good catalyst that promotes proton transfer from water to SO₃, and according to both experimental and theoretical reports, the rate constant of SO₃ + 2H₂O reaction is approximately 10⁻¹⁵ cm³·molecule⁻¹·s⁻¹ [60, 61, 62]. HO₂ [63], HCOOH [62, 64], H₂SO₄ [49], HNO₃ [65] can replace the role of the second water molecule and exhibit a similar catalytic effect on the formation of sulfuric acid (H₂SO₄).

4. Concluding remarks

In contrast to the rapid development of chemistry, physics, and biology and other disciplines, soil science remains almost stagnant in the past few decades, and to best of our knowledge, no breakthroughs have been reported for rather a long time. Despite that, no one can deny the vital significance of soil to our life, and soil science has been widely acknowledged as the final frontier.

The slow progresses for soil science, we think, should be attributed to two reasons: (1) wrong and outdated perceptions. A majority of soil researchers mistakenly believe that all knowledge worth knowing about soils has already been understood and no revolutionary progresses would take place; (2) complex systems. Soils are very structurally complicated and there are multiple factors to co-function, which makes it very difficult to characterize by experimental techniques. Computer simulations have unique advantages to handle complex systems, while currently its role in soil science is far from being recognized. In this chapter, two examples are elaborately discussed with regard to application of DFT calculations in soil science: one focuses on the adsorption of uranyl onto mineral surfaces, and the other involves the adsorption of SO₂ onto mineral surfaces and reaction mechanisms to form acid rain. It can be seen from these discussions that DFT calculations are able to provide useful and detailed information about the adsorption, interaction and reactions at the atomic level that greatly promote our understanding about soil science.

With advent of high-performance computing platforms, the same DFT calculation tasks of 10 years ago can now finish within a remarkably shorter time, even if you increase model size (periodic model is also an option), consider solvent effects by adding explicit solvent molecules or/and choose more accurate methods. The methodological developments regarding to DFT calculations have also made remarkable progresses over the recent three decades, and as a result, thermodynamics and reaction barriers can now be predicted with nearly chemical accuracy (≤1 kJ/mol). In consequence, computer simulations including DFT are the right key to unravel the complicated phenomena and processes occurring within soils; e.g., with DFT calculations, the aggregation mechanisms of “real” soils and the driving force therein were unveiled at an atomic level [6].

In addition to DFT methods, there are a number of other computational methods, such as QM/MM and Molecular Dynamics (MD). The choice of suitable computational methods is strongly recommended. We are pleased to see the birth of the ClayFF force-field [66] that was developed specially for clay minerals and the capability of the ReaxFF force-field [67] to handle reaction mechanisms. The 2013 Noble Prize in Chemistry was awarded to the development of “multiscale methods for complex systems,” and now it is time to apply these methods to tackle complex soil systems.

Acknowledgments

This work was sponsored by the National Natural Science Foundation of China (21473137), the Fourth Excellent Talents Program of Higher Education in Chongqing (2014-03) and the Natural Science Foundation Project of CQ CSTC, China (cstc2017jcyjAX0145).