The Binding of Phosphorus Species at Goethite: A Joint Experimental and Theoretical Study

: Knowledge of the interaction between inorganic and organic phosphates with soil minerals is vital for improving the soil P fertility. To achieve an in-depth understanding we combined adsorption experiments and hybrid ab initio molecular dynamics simulations to analyze the adsorption of common phosphates, i.e. orthophosphate (OP), glycerolphosphate (GP) and inositolhexaphosphate (IHP), onto the 100 surface plane of goethite. Experimental adsorption data per mol P-molecule basis ﬁtted to the Freundlich model show the adsorption strength increases in the order GP < OP < IHP, and IHP adsorption being saturated faster followed by GP and OP. Modeling results show that OP and GP form stable monodentate ( M ) and binuclear bidentate ( B ) motifs with B being more stable than M , whereas IHP forms stable M and 3M motifs. Interfacial water plays an important role through hydrogen bonds and proton transfers with OP/GP/IHP and goethite. It also controls the binding motifs of phosphates with goethite. Combining both experimental and modeling results, we propose that the B motif dominates for OP, whereas GP forms M and IHP forms a combination of M and 3M motifs. The joint approach plausibly explains why IHP is the predominant organically bound P form in soil. This study could be considered as a preliminary step for further studies for understanding the mechanisms of how microbes and plants overcome the strong IHP–mineral binding to implement the phosphate groups into their metabolism.


Introduction
Phosphorus (P) adsorption capacity of soil is a crucial factor affecting the P immobilization process and thus soil fertility, and fate of P in natural environments [1][2][3][4]. This adsorption capacity, that is defined as the amount of adsorbed substance (adsorbate) reached in a saturated solution for a specific adsorbent [5], is fundamentally influenced by the diverse interactions of P with soil constituents involving free metal ions [6,7], soil minerals [8][9][10][11], and soil organic matter [12,13]. In particular, the strong interaction between phosphates and soil mineral surfaces, and especially Fe-and Al-(oxyhydr)oxides, plays a very important role in controlling this process [8,[14][15][16][17][18]. The Fe-oxides are widespread in surface environments and constitute a major component of highly weathered soils and sediments. They have attracted considerable attention due to their high P-adsorption capacity [19][20][21]. Phosphate adsorption on iron oxide surfaces shows a biphasic behavior consisting of a rapid and strong ligand exchange step followed by a slower step [22][23][24][25]. The latter phase was described previously by the formation of monodentate (M) complexes and conversion into bidentate (B) complexes, the competition with other anions, and/or precipitation events [26]. Alternatively, the slower phase has been assigned to a diffusion process of phosphate ions from an outer-sphere complex to the surface [16,23,27,28]. 3 of 19 resulting ferrihydrite was aged for 120 h at 55 • C to be converted to goethite. Then the suspension was subjected to pH adjustment to pH 6 by adding 0.1 M HCl and followed by centrifugation-washing cycles with distilled water until the electrical conductivity was lower than 10 µS cm −1 . The so-prepared goethite was freeze-dried and stored as powder for further analysis and experiments.
The results of X-ray diffraction (XRD) analysis confirmed the well-crystallinity of goethite samples, see Fig. ??. The measured specific surface area (SAA) of goethite by Brunauer-Emmett-Teller (BET, N 2 adsorption; Nova 4000e, Quantachrome, USA) method is 64.5 m 2 g −1 . This value lies in the reported common range values of SSA of goethite in the literature [27,33,43]. The mean hydrodynamic diameter of goethite particles is 1.6±0.2 µm as determined by a particle sizer (Zetapals, Brookhaven, USA) on a 100 mgL −1 suspension treated with a dispersant agent (0.01 mM Na 4 P 2 O 7 , pH 8) and sonicated for 30s (Labsonic M, Sartorius Stedim). The measured point of zero charge (PZC) of goethite particles is 6.3, see Fig. ??.
In the present study, adsorption of three different phosphate compounds involving inorganic orthophosphate (OP, potassium dihydrogen phosphate, CAS number: 7778-77-0) and organic phosphates (α-glycerol phosphate (GP, CAS number: 17603-42-8) and myo-inositol hexakisphosphate (IHP, CAS number: 83-86-3)) on goethite was performed. Here, a stock solution of 10 mM P of each P compound was prepared in 0.01 M CaCl 2 solution. For each experiment, a 200 mg sample of goethite was weighed into a 50 mL centrifuge tube and equilibrated with 40 mL of each initial P concentration.
Here twelve different initial P concentrations (0, 0.05, 0.1, 0.15, 0.2, 0.3, 0.5, 1, 2, 3, 5, and 10 mM P) were considered. Each initial concentration was prepared by dilution of the P stock in 0.01 M CaCl 2 background electrolyte solution adjusted at pH 5. After an equilibrium period of 24 h at 25 • C under end-over-end shaking at 20 rpm, samples were centrifuged at 4500g for 15 min. The P content in the filtered supernatant was quantified by inductively coupled plasma-optical emission spectroscopy (ICP-OES) method.
It is noteworthy that all adsorption experiments were performed in triplicate and data were presented as the means of three repeats. Moreover, all used P compounds (myo-inositol hexakisphosphate: C 6 H 18 O 24 P 6 ; α-glycerol phosphate: C 3 H 9 O 6 P; potassium dihydrogen phosphate: KH 2 PO 4 ) in the present experiments were of analytical grade chemicals and purchased from Carl Roth GmbH and Sigma-Aldrich. Working solutions were prepared fresh daily by adding accurate quantities of the prepared P stocks into 0.01 M CaCl 2 solution.
To describe the adsorption behaviors of P compounds on goethite, the adsorption data were fitted to the Freundlich [44], Langmuir [45], and Temkin [46] models. The Freundlich model provides a two-parameter equation that describes the relationship between the equilibrium concentration and the adsorbed one onto heterogeneous surfaces through the following equation: that is rewritten in a linear form as follows: where Q ads is the amount of adsorbate adsorbed per unit mass (or surface area) of adsorbent, C e is the adsorbate equilibrium concentration in solution, K f is Freundlich adsorption constant (sometimes it is defined as Freundlich unit adsorption capacity), and n f is Freundlich exponent [44]. In the present study, Q ads , C e , and K f are expressed in µmol m −2 , µmol L −1 and µmol 1−n f L n f m −2 , respectively. Notice that Freundlich model assumes that the adsorption enthalpy depends on the amount of adsorbate. In the limit of small Q ads where the adsorption enthalpy should not depend on Q ads one could describe the isotherm by a Langmuir model as well. The Langmuir adsorption theory assumes that the adsorbate forms a monolayer on a homogeneous adsorbent surface. The following equation expresses the Langmuir isotherm: that is rewritten in a linear form as follows: where Q max is the maximum amount of the adsorbate which is required to form a monolayer by complete saturation of all binding sites (it is defined also as the maximum monolayer coverage capacity or simply as monolayer capacity [5], µmol m −2 ), K l is the Langmuir adsorption constant (L µmol −1 ) which is mainly related to the adsorption energy [45]. The Temkin model exhibits a factor considering the adsorbent-adsorbate interaction with a uniform distribution of binding energies. The model is expressed by the following equation: where R is the universal gas constant (8.314 J K mol −1 ), T is the absolute temperature (K), b T is Temkin isotherm constant (J mol −1 ), A T is Temkin isotherm equilibrium binding constant (L µmol −1 ), and B T is constant related to the heat of adsorption. For better understanding for the P adsorption behaviour, the calculated adsorption coefficients are intrepreted on three different bases (per mol P (i.e. per mole of P element), per mol molecule (i.e. per mole of P-containing molecule (OP, GP, and IHP)) , and per mass (i.e. per mass (mg) of P-containing molecule)).

Molecular Modeling and Computational Details
The current molecular modeling approach for the binding process of phosphates at the goethite-water interface is illustrated in Fig. 1. Here, three different abundant phosphates are considered, i.e. OP (H 3 PO 4 ), GP (C 3 H 9 O 6 P), and IHP (C 6 H 18 O 24 P 6 ), to study the complexation reaction of each phosphate with the goethite surface in the presence of water (see Fig. 1a-c). Each constructed goethite-phosphate-water complex model (for example see Fig. 1d) consists of 1-the goethite surface, 2-a phosphate compound with a specific binding motif to the surface, and 3-water molecules surrounding the goethite-phosphate complex. Here, the 100 goethite surface is considered which is one of the most abundant goethite surface planes in soil systems [33,47]. This goethite surface plane is modeled by repetition of the goethite unit cell (lattice constants : a = 9.95, b = 3.01, c = 4.62 Å) as 1a × 8b × 5c which consists of 640 atoms (160 Fe, 160 H, and 320 O atoms). The binding motifs formed between these phosphates and the goethite surface are modeled based on previous literature studies [9][10][11]41,[48][49][50]. The initial motifs of OP and GP with the goethite surface include monodentate (M, 1Fe+1O one covalent bond between phosphate non-protonated O atom and a surface Fe atom) and bidentate (B, 2Fe+2O i.e. two covalent bonds between two phosphate O atoms, one protonated and other non-protonated and two surface Fe atoms), see Fig. 1e-f. For IHP, additional, 4M (4Fe+4O i.e. four covalent bonds between four separate phosphates non-protonated oxygens and four surface Fe atoms) is considered since IHP is known to interact with goethite through multiple phosphate groups, see Fig. 1g. To simulate the effect of water on the goethite-OP/GP/IHP complexes, each modeled binding motif is solvated with water molecules at a density of ≈ 1g cm −3 perpendicular to the studied surface plane with height of ≈ 18 Å by using the visual molecular dynamics (VMD) package (see Fig. 1d) [51]. To ease the discussion about the interactions, oxygen atoms bonded surface Fe atoms are denoted as O p .
Due to the large size for the modeled goethite-OP/GP/IHP-water complexes (each model consists of more than 1200 atoms), simulation of these systems with pure ab initio methods is computationally expensive and thus not feasible [52,53]. Therefore, the present models will be studied by molecular dynamics (MD) simulations based on the electrostatic embedding QM/MM hybrid approach. Here, the reactive part (system of interest) is described at a QM level of theory, while the remaining part is described using MM. The QM part includes 1-OP/GP/IHP, 2-the goethite top surface layer (160 atoms), and 3-water molecules within a layer of about 10 perpendicular to the surface (≈ 55 molecules). The total QM box size is 22 × 8b × 5c i.e., 22 × 24.08 × 23.1 Å. The QM part is simulated with density functional theory (DFT) as implemented in the quickstep code in CP2K [54]. A hybrid Gaussian and plane wave (GPW) dual basis method is used with ionic cores described by Goedecker-Teter-Hutter (GTH) pseudopotentials [55] in combination with the Perdew-Burke-Ernzerhof (PBE) [56] exchange correlation functional and the Grimme D3 empirical dispersion correction [57]. The valence electrons of all atoms are defined with the double-ζ valence polarized MOLOPT (DZVP-MOLOPT-SR-GTH) basis set except water for which the single-ζ valence (SZV-MOLOPT-SR-GTH) basis set is used to reduce computational cost [58]. The SCF convergence threshold was chosen to be 10 −4 hartree. The MM part is described by classical force fields (FF) with FIST module in CP2K [59]. The goethite surface is modeled with the CLAYFF FF [60] while water is modeled with the single point charge (SPC) water model [61] and OP/GP/IHP with CHARMM FF using the SwissParm tool [62]. Both CLAYFF and CHARMM FFs are compatible with the SPC water model. The interaction between the QM and MM parts in CP2K is implemented using the Gaussian expansion of the electrostatic potential method (GEEP) [63], wherein the MM charge is distributed by defining it using Gaussians instead of point charges to avoid electron spilling. QM/MM based canonical (NVT i.e. constant number of atoms N, volume V and temperature T) MD simulations are performed for 25 ps with a 0.5 fs time step while the temperature was maintained at 300K using velocity rescaling thermostat (CSVR) [64]. The first 10 ps of each trajectory are assigned for the equilibration and the remaining 15 ps (production trajectory) are used for analysis. The interaction energy between each phosphate and the goethite surface is calculated for every 100 fs (i.e. 150 snapshots) along the production trajectory by using: Here, E G−P complex , E P , and E G denote the total electronic energies of the goethite-phosphate complex, the phosphate (OP/GP/IHP), and the goethite surface, respectively. The interaction energies involving water are divided by the number of water molecules in the simulation box for better comparison. The basis set superposition error (BSSE) in interaction energies is corrected using counterpoise scheme [65]. Despite a few data scatter in OP adsorption, Fig. 2 shows that almost the maximum adsorption capacities are reached for goethite and the selected range of P concentrations is sufficient to achieve the equilibrium. Here, the sequence of adsorbed P (Q ads ) values per mol P is IHP > OP > GP (see Fig. 2). For most cases, the coefficients of determination (R 2 ) values indicated that the Freundlich equation fitted the adsorption data better than the Langmuir and Temkin models, see Table 1. The Freundlich adsorption constant (capacity) K f ranged from 0.22 to 4.79 µmol 1−n f L n f m −2 . The order of K f values suggests that IHP exhibits the strongest adsorption and the highest capacity at the goethite surface, followed by OP and GP (GP < OP < IHP). The magnitude of the Freundlich exponent n f , that ranged from 0.07 to 0.29 (see Table 1), gives an indication that the sorption mechanism is dominated by adsorption and not absorption [66,67]. Furthermore, the exponent points to the diversity of the energies associated with adsorption of P compounds on the goethite surface. Moreover, n < 1 for all cases indicates that upon increasing the P concentration/loading the binding energy between the surfaces and P compounds is reduced. The order of the n f values (IHP < GP < OP) suggests that the binding energy decreases with increasing the P loading in the order OP < GP < IHP. This means that affinity of the goethite surface to adsorb/bind a P molecule, with increasing the P concentration, increases in the order IHP < GP < OP. According to the Langmuir model, the maximum monolayer adsorption capacities (Q max ) are 1.2, 7.64 and 8.35 µmol m −2 , respectively, for GP, OP and IHP (see Table 1). This shows that the order of saturation of the goethite surface with P per mol P is GP < OP < IHP. This trend fits well with that one observed based on the Freundlich K f values. The Langmuir parameter K l increased in the order 0.001 (OP) < 0.003 (GP) < 0.06 (IHP) L µmol −1 . This constant is mainly related to the adsorption energy and could give information on how strong (i.e. strength) the goethite-P interaction/binding process. Therefore, based on Langmuir model, one expects that the Preprints (www.preprints.org) | NOT PEER-REVIEWED | Posted: 8 February 2021 doi:10.20944/preprints202102.0171.v1

Adsorption Isotherms
goethite-P interaction increases in the order OP < GP < IHP. The same order but with different values was obtained from Temkin binding constant A T which is also related to the binding strength, see Table  1. The Temkin B T values for OP (0.75), GP (0.13) and IHP (0.52) suggests that the heat of adsorption increases in the order GP < IHP < OP. Regardless the Temkin constant B T , all other parameters from the represented isotherm models in the present contribution refer to stronger adsorption and higher capacity for the IHP case compared to OP and GP cases by considering the number of P moles. Considering the size and chemical structure of IHP (six phosphate groups), GP and OP (only one phosphate group), the adsorption data were presented as well in terms of whole P molecular system. Specifically, both Q ads and C e concentrations and all related fitted adsorption isotherm parameters are represented per mole and per mass of the whole P-containing molecule (OP, GP, and IHP), see Table  1. It should be noted that n f and R 2 values for all P compounds do not change up on changing the interpretation/representation of P concentration (Q ads and C e ). In addition, the order of all adsorption parameters (k f , n f , k l , Q max , B T , A T ) does not change for OP and GP by different interpretations for the P concentration (i.e. per mol P verus per mol molecule and per mass bases). Here, the only change was observed for the IHP case. By considering the number of moles of the P adsorbed molecules (i.e. per mol molecule), the K f , Q max and B T values for IHP decreased by a factor of six, while the corresponding K l and A T values increased by a factor of six (see Table 1). In principle, this means that the goethite adsorption capacity will decrease to one sixth by considering the whole IHP molecule compared to considering six P moles of each IHP mole. Consequently, the binding/adsorption strength between the goethite surface and each IHP mole will increase by a factor of six compared to the case of considering only one mole of P. In addition, the orders of K f (GP < OP < IHP), K l (OP < GP < IHP) and A T (OP < GP < IHP) do not change compared to the mol P basis. In contrast, the orders of Q max (GP < IHP < OP) and B T (IHP < GP < OP) changed comparing with considering the mol P cases. Same orders were obtained by considering the P concentration as mass of the adsorbed P molecule. With respect to the K l order, one can suggest that GP has stronger interaction than OP with the goethite surface. In contrast, the B T order suggests that OP has stronger interaction than IHP. In a study evaluating adsorption properties of goethite, Celi et al. [68] reported the adsorption ratio of 3:2 between IHP and OP at pH 4.5. Later, Martin et al. [69] confirmed the greater affinity of IHP than of OP for goethite, indicating that the maximum amount of adsorbed P were 3.6 and 2.4 µmol P m −2 for IHP and OP, respectively. Likewise, Celi and Barberis [70] found that goethite shows a higher affinity for IHP (deduced from Langmuir K l values) than for other organic and inorganic phosphates. Li et al. [71] reported a maximum GP adsorption of 1.95 µmol P m −2 at pH 5 which is close to the value of 1.20 µmol P m −2 observed here. However, the lower adsorption of GP as compared to that reported by Li et al. [71] may be due to the differences in experimental conditions. Specifically, Li et al. [71] conducted their adsorption experiment by β-glycerophosphate at a goethite sample with SAA of 46.4 m 2 g −1 and PZC of 9.2 in 0.1 M KCl as background electrolyte. Overall, both organic and inorganic phosphates exhibited affinity for goethite. Parfitt [29] compared the adsorption of OP on goethite in CaCl 2 solution and reported that larger amounts of specific adsorbed P reflect more reactive and available adsorption sites on the goethite. This was explained by thermal gravimetric analyses that indicated that goethite contains 24% structural OH [71]. Similarly, Guzman et al. [72] showed the affinity for phosphate (exponent coefficient of Freundlich equation) for goethite on a surface area basis. Details of IHP adsorption on goethite have been extensively studied [68,69,73,74]. Ognalaga et al. [48] investigated the adsorption of IHP on synthetic goethite and attributed it solely to the inner-sphere complexation of phosphate groups with reactive OH groups. Celi et al. [68] and Ognalaga et al. [48] suggested that four of the six phosphate groups of each IHP molecule were involved in bonding to goethite while the other two remained free.

Orthophosphate
The MD simulations showed that OP maintained the initial motifs (see Fig. 1e-f) and formed stable M and B motifs throughout the simulation trajectory, see Fig. 3a-d and Fig. 3e-h, respectively. The time averaged Fe-O p bond length and Fe-P distance observed for the B motif are shorter than those for the M motif which is probably due to higher stability in the former case, see Table 2. Here, the total interaction energy between OP and the goethite surface of the goethite-OP complex for the B and M motifs are -82 and -35 kcal mol −1 , respectively. The interaction energy per bond of goethite-OP in the B motif being higher than in the M motif suggests that the alignment of OP with the goethite surface favors the formation of an additional Fe-O p covalent bond in the B motif. For both M and B motifs, two protons are transferred from OP to water (see Fig. 3a-b and e-f), forming an average of five HBs each. Consequently, the calculated OP-water interaction energy per water molecule for M (-1.7 kcal mol −1 ) and B (-1.6 kcal mol −1 ) motifs is essentially the same. In addition, proton transfers are observed from goethite to water, see Fig. 3a and e-f. The proton transfers observed here occurred only during the equilibration phase of the complexes and OP remained twice deprotonated in the production trajectory for both cases. The goethite-water interaction energy per single water molecule (-4.8 kcal mol −1 ) is lower than the goethite-OP interaction energy for both M and B motifs. The goethite-water interaction energy is higher than OP-water interaction energies in both motif cases because of goethite's larger solvent accessible surface area (SASA), proton transfer events from goethite to water (see Fig. 3), formation of multiple Fe-O H 2 O M motifs (on average 15 of 40 surface Fe atoms, see Fig. ??) and HBs between goethite and water.

Glycerolphosphate
GP is aligned perpendicular to the goethite surface to form initial M and B motifs as shown in Fig. 1d and Figs. ??c-d. The MD simulations show that these initial motifs are stable during the whole trajectory (see Fig. 3i-l and Fig. 3m-p, respectively)  see Table 2. GP forms an average of six and eight HBs with water in the M and B motifs, respectively. It also exhibits multiple proton transfers to water, two in the M motif case and three in the B motif. The proton transfer events observed for both cases occurred during the equilibration phase. The total interaction energy between GP and goethite for the B motif (-76 kcal mol −1 ) is higher than for the M motif (-48 kcal mol −1 ) due to additional covalent bond in former motif. The higher overall goethite-GP interaction energy of the B motif indicates that it is a more favorable motif than the M one. But the goethite-GP interaction energy per bond for the B motif is lower than for the M motif, suggesting comparatively unfavorable conditions for Fe-O p bonds in former motif. This is in contrast to the case of OP. The difference comes from the glycerol group in GP. Both OP and GP interact through the phosphate, but through the glycerol group, GP has a higher SASA and exhibits stronger interaction with water than OP (OP/GP-water interaction energies will be discussed below). In the GP B motif case, two oxygens of GP's phosphate group are bound to two individual Fe atoms and when the glycerol group interacts with flexible water atoms, GP fluctuates with water while the surface Fe atoms competitively pull O p oxygens towards the surface as a counteraction. This may induce strain in the Fe-O p bonds for the B motif and, therefore, the goethite-GP interaction energy per bond for the B motif is less compared to the M motif. The stretch in the Fe-O p bonds found in the GP B motif is larger than in the OP B motif, see Fig. ??, which signifies a strain in the Fe-O p bonds in the former case. In addition, GP is found to oscillate occasionally in a see-saw type of motion over the surface which randomly stretches one of the Fe-O p bonds to its extreme (maximum Fe-O p bond length found ≈ 2.4 Å), see Fig. ??. A similar observation was made for GP on 100 diaspore surface [41] which is isomorphous with goethite [33]. The study by Xu et al. [75] showed that GP induces more negative charge onto hematite mineral surface than OP which might also influence the strength of the Fe-O p covalent bonds. Regarding GP-water interaction energies, GP exhibited a slightly stronger interaction with water for the B motif (-2.8 kcal mol −1 ) case than for the M motif one (-2.5 kcal mol −1 ). This is due to additional proton transfer observed from GP to water for the B motif case (see Fig. 3m-o). The GP-water interaction energy here is higher than that of OP-water due to additional glycerol group in GP which contains two polar OH groups that interact strongly with water. Similar to goethite-OP-water complexes, proton transfer is observed from goethite to water in both motifs here, see Fig. 3i and n.

Inositolhexaphosphate
IHP initially was aligned perpendicular to the goethite surface to form M and B motifs (see Fig. ??a-b), and parallel to the goethite surface to form the proposed 4M motif (see Fig. 1g). In contrast to goethite-OP-water and goethite-GP-water cases, transformation of initial binding motifs between IHP and the goethite surface with time is observed. The simulation results show that the initial M  Table 2). The M(1) motif's goethite-IHP interaction energy is higher than that for the M(2) one (see Table 2). This is probably because of the strain in the Fe-O p bond in M(2) motif due to intramolecular HBs between the phosphate group bound to Fe and its adjacent phosphate groups (see Fig. 4h). The 3M motif has a higher total interaction energy than both M motifs due to additional covalent bonds and proton transfer from IHP to the goethite surface. But the interaction energy per bond between IHP and goethite for the 3M binding motif is less than for the M(1) motif. Thus, in spite of the 3M motif being the more favorable motif compared to the M motifs, the alignment of IHP with goethite as in the 3M motif weakens the overall strength of the individual Fe-O p bonds. Even though an additional proton transfer is observed from IHP to water for the M(1) motif, the IHP-water interaction energy (-7.3 kcal mol −1 ) is very close to that for the M(2) motif case (-7.2 kcal mol −1 ). This could be probably due to higher number of time averaged HBs observed for the M(2) motif case (25) compared to the M(1) case (23). The IHP-water interaction energy observed for the 3M motif case (-3.1 kcal mol −1 ) is lower than those for both M cases due to fewer proton transfer processes (2) and HBs (15) observed in the 3M motif case.
The stronger interaction of IHP with water is vital for understanding the reason for transformation of initial B and 4M motifs to M(2) and 3M final motifs, respectively. For initial B motif case, when IHP interacts with flexible water molecules it fluctuates with respect to goethite surface which might induce strain in one of the Fe-O p covalent bonds leading to its dissociation. In addition, an intramolecular HB is observed between the Fe bonded phosphate group and its adjacent phosphate group, see Fig. 4l. In our previous work on IHP interaction at the diaspore(100)-water interface (diaspore isostructural with goethite), [41] the initial 4M motif was transformed into a 2M motif. The reason is that the Al-O p covalent bonds are inclined and restricted between diaspore's surface hydroxyl groups. Moreover, when IHP interacts with water and fluctuates during the equilibration phase, the restricted Al-O p covalent bonds under unfavorable conditions do not have an alternative but to dissociate and move towards water. The same reason holds true here for the transformation process of the initial 4M motif to a stable 3M binding motif.

Discussion of Experimental and Modeling Results
In what follows we will correlate the present modeling results and the adsorption isotherms as per mole of P-compound. First, let us recollect that K f , K l and B T relate to the phosphate's binding energy and that n f relates to the binding energy of the next incoming P molecule to surface (i.e. binding affinity). Since the Freundlich isotherm provided the best fit to the experimental data for all cases (see Table 1), we will discuss the theoretical results in Table 2 in the context of K f and n f . Analyzing the strength of phosphate interaction shows that the K f values are in line with the order of overall binding energies observed here, i.e. GP B < OP B < IHP 3M. The K f value for OP adsorption is more than twice that of the GP adsorption, suggesting a significant difference between OP and GP interaction with goethite. In the light of these results combined with the theoretical binding energy values (see Table 2), one could suggest that OP might predominantly form B motif and GP might form M motif as the predominant binding motif. Therefore, the binding strength order could be updated as GP M < OP B < IHP 3M. Comparing K f values for IHP and OP, the K f for IHP adsorption is about 1.58 times of that for OP adsorption which is very close to the reported adsorption ratio by Celi et al. [68] and Martin et al. [69]. Correlating this ratio to the present calculated binding energies would suggest that one should rather use the average of all three bonding motifs, i.e. both M and 3M binding motifs for the adsorbed IHP molecule are present at the goethite surface. Therefore, based on both experimental and theoretical binding strength values, one could suggest that the overall binding strength increases in the order GP (M) < OP (B) < IHP (M + 3M).
The order of n f values (IHP < GP < OP, see Table 1) suggests that IHP adsorption leads to a faster saturation of the goethite surface than GP and OP. Even though IHP binds to goethite through a few phosphate groups, the remaining phosphate groups (often deprotonated) would induce more negative charge to the surface compared to OP and GP cases [68,75]. Moreover, it was found that GP leads to a more negative charged surface than OP [75]. This explains the reason why the n f value for GP is lower than for OP.
Let us compare the present modeling results with previous experimental and theoretical studies. In case of OP experimental [36,76,77] and theoretical [10,40] resutls pointed to a formation of both, i.e. M and B, motifs with goethite. Specifically, for the 100 goethite surface studied here, Kubicki et al. [40] and Ahmed et al. [10] demonstrated that OP often exists in a doubly deprotonated state forming both M and B binding motifs but with the predominance of the B one. This is confirmed by the present MD results, i.e. OP has formed stable M and B motifs and remained twice deprotonated in the production Preprints (www.preprints.org) | NOT PEER-REVIEWED | Posted: 8 February 2021 doi:10.20944/preprints202102.0171.v1 trajectory with the B motif exhibiting higher interaction energy with the goethite surface compared to the M motif. Even though the B motif is the more stable motif under the present conditions, it is not always the dominant one, a topic which has been highly debated in literature [39,40,76,77]. For instance, Persson et al. [77] interpreted their FTIR spectra of the OP adsorption at goethite and hematite by the formation of only M binding motifs with different protonation states. This analysis was done as function of total phosphate concentration, pH, and time. In contrast, Tejedor-Tejedor and Anderson [36] proposed that the B motif is the dominant motif, for low surface coverage of OP on goethite over the pH range of 3.6-8.0. Further, they proposed that the M motif exists at low surface coverage but as a non-dominant motif. By using DFT calculations, Kwon and Kubicki [39] concluded that the B motif is dominant between pH 4-6. Recently, a joint approach by Ahmed et al. [49] involving adsorption experiments and DFT simulations showed that the M motif is dominant at both extremely low and high pH values, while the B motif is dominant in the intermediate pH range. Keeping in mind that the present simulations correspond to a low surface coverage scenario and acidic pH, one can conclude that the stable M and B motifs with abundance of the B one here are in line with the studies of Tejedor-Tejedor and Anderson [36], Ahmed et al. [49], and Abdala et al. [76]. Moreover, the time averaged Fe-P distance observed in both motifs are within the range of values reported in literature, thus giving further support for the present model. For instance, an extended X-ray absorption fine structure study of OP on goethite [76] showed that Fe-P distances (Å) for M and B motifs are 3.6 and 3.28, respectively. The Fe-P distances from other studies are in the range of 3.48-3.55 Å for the M motif and 3.13-3.37 Å for the B motif [10,40,78,79].
Compared to OP, GP adsorption on goethite has not been extensively studied and hence, the information about its binding motifs with goethite is limited. Li et al. [71] showed that GP forms inner-sphere complexes with goethite by replacing water and hydroxyl groups from surface active sites. The study also proposes that GP forms only M motifs at the goethite surface based on FTIR spectra analysis. The B motif is sterically hindered by the organic moiety. Also, Hartree-Fock simulations and FTIR spectra studies by Persson et al. [80] of monomethyl phosphate (CH 3 -H 2 PO 4 ), an organic P with a single phosphate group like GP, showed that it forms mostly M motifs with goethite. Adsorption isotherms and FTIR spectra by Sheals et al. [34] proposed that glyphosate binds predominantly through an M motif but B motifs might be formed at low P concentration and neutral pH. By using periodic DFT based MD simulations, Ahmed et al. [9] studied the glyphosate binding process at the goethite-water interface by considering three goethite surface planes. The results indicated that the M binding motif is the most dominant one at the 100 goethite surface plane, although both M and B motifs exist at the 100, 010 and 001 surface planes. Summing up these literature data, the M motif is the dominant motif for GP like molecules at the 100 goethite surface. Regarding the Fe-P distances observed in both GP motifs here, the distances are within the range of Fe-P distances observed for OP on goethite in literature [10,40,76,78,79]. This suggests that the GP's phosphate group interacts with goethite surface in a similar way to the OP case.
There is no consensus in the literature about the number of IHP's phosphate groups that bind to goethite surface or about dominant binding motifs. The binding motifs observed for goethite-IHP complexes in the current study (M and 3M) have been found in literature before [11,41] but it is not clear which motif is the dominant one. For instance, 3M [11,81], 2M [11,41] and 1M [11] are different binding motifs observed for IHP on minerals. In contrast, Johnson et al. [74] proposed that IHP interacts with goethite by forming outer-sphere complexes. Interestingly, none of the above mentioned studies for IHP suggest that it forms a bidentate binuclear (B) motif. Similar to our previous studies [11,41], we found the initial B motif of IHP being not stable due to the strong interaction of IHP with water and intramolecular HBs. Therefore, we propose that IHP phosphate groups form M motif with goethite surface.
De Groot and Golterman [82] showed that IHP adsorption onto goethite had an inhibitory effect on OP adsorption and it can release adsorbed OP from the goethite surface. The same could be inferred from binding energies here, where the goethite-IHP binding energies are larger than for goethite-OP complexes. The goethite-GP binding energies are also less than goethite-IHP and one might suggest that IHP might replace GP as well from goethite surface. The same could be inferred from the adsorption strengths order calculated from experiments here. The Fe-P distances observed here show that they are close to the Fe-P distances for OP on goethite [10,40,76,78,79] which suggest that IHP's organic moiety might not influence the individual phosphate groups interaction with goethite but only the conformational flexibility of the overall binding motif [11,41].This comes in accord with the observation by Celi et al. [68] that phosphate groups of IHP react with goethite surface similar to OP.

Summarizing Discussion
The current study contributes to the efforts in understanding the interaction of phosphates with soil minerals. Here, both experimental and theoretical approaches are adopted to characterize inorganic (OP) and organic (GP, IHP) phosphates interaction with abundant and reactive goethite mineral. The goethite-OP/GP/IHP-water complexes are simulated with the multiscale QM/MM method which provides molecular level insights into adsorption experiments performed for OP/GP/IHP on goethite. The binding energies and interaction mechanisms of phosphates adsorption on goethite from the modeling study are correlated to adsorption data fitted to the Freundlich isotherm, which provided a uniformly better fit than the Langmuir and Temkin models. The model coefficients provided an overview about the pattern of phosphates interaction with goethite and the QM/MM simulations zoom into the molecular level attributes that build these pattern.
The modeling results show that, OP forms stable M and B motifs and OP B motif has higher interaction energy per bond than OP M motif. This suggests that goethite-OP interaction favors the additional covalent bond which makes the OP B motif more stable. For goethite-GP complexes, the order of interaction energies is same as OP case, except that GP B motif interaction energy per bond is lower than GP M motif. Compared to OP, the GP B motif is weaker than OP B motif due to a strain in Fe-O p bonds in former case. The strain is due to GP interaction with water molecules. In addition, Xu et al. [75] study shows that adsorption of GP induces more negative charge on hematite surface than OP which might strain the Fe-O p bonds further. Therefore, literature [34,71] and the calculated binding energies in the present contribution show that GP's dominant motif might be M or B motif depending on the GP interaction with environmental molecules. IHP forms M and 3M motifs with multiple intramolecular HBs between adjacent phosphate groups [74]. IHP's 3M motif exhibits the strongest binding energy with the goethite surface among all goethite-OP/GP complexes here. Therefore, we propose that the 3M motif might be the most dominant motif at low surface loading.
The adsorption data from experiment fitted to Freundlich model show that the order of adsorption strength (K f ) is GP < OP < IHP whereas the order of incoming P binding strength with surface (n f ) is IHP < GP < OP. Based on the mangnitude and order of experimental K f values combined with theoretical binding strength values, one could suggest that the overall binding strength increases in the order GP (M) < OP (B) < IHP (M + 3M). This shows that GP might often form M motif whereas OP might form B and IHP form both M and 3M motifs. The n f order suggests that IHP adsorption on goethite saturates its surface faster than for the GP and OP cases. The fact that GP adsorption induces more negative charge on goethite than OP supports our suggestions that GP would often form M and not B motifs and eventually its interaction energy would be less than OP.
Eventually, the present experimental and theoretical results, in agreement with the pertinent literature, plausibly explain why IHP is the predominant organically bound P form in soil. A challenge for further studies and application is to understand and explore the mechanisms by which microbes and plants can overcome the strong IHP-mineral binding and incorporate the phosphate groups into their metabolism.