Interaction of Nitrite Ions with Hydrated Portlandite Surfaces: Atomistic Computer Simulation Study

The nitrite admixtures in cement and concrete are used as corrosion inhibitors for steel reinforcement and also as anti-freezing agents. The characterization of the protective properties should account for the decrease in the concentration of free NO2− ions in the pores of cement concretes due to their adsorption. Here we applied the classical molecular dynamics computer simulation approach to quantitatively study the molecular scale mechanisms of nitrite adsorption from NaNO2 aqueous solution on a portlandite surface. We used a new parameterization to model the hydrated NO2− ions in combination with the recently upgraded ClayFF force field (ClayFF-MOH) for the structure of portlandite. The new NO2− parameterization makes it possible to reproduce the properties of hydrated NO2− ions in good agreement with experimental data. In addition, the ClayFF-MOH model improves the description of the portlandite structure by explicitly taking into account the bending of Ca-O-H angles in the crystal and on its surface. The simulations showed that despite the formation of a well-structured water layer on the portlandite (001) crystal surface, NO2− ions can be strongly adsorbed. The nitrite adsorption is primarily due to the formation of hydrogen bonds between the structural hydroxyls on the portlandite surface and both the nitrogen and oxygen atoms of the NO2− ions. Due to that, the ions do not form surface adsorption complexes with a single well-defined structure but can assume various local coordinations. However, in all cases, the adsorbed ions did not show significant surface diffusional mobility. Moreover, we demonstrated that the nitrite ions can be adsorbed both near the previously-adsorbed hydrated Na+ ions as surface ion pairs, but also separately from the cations.

To simulate the interface of portlandite with NaNO 2 aqueous solution, a system consisting of 11 × 10 × 8 portlandite unit cells was cleaved parallel to the (001) plane and cut in the ab plane of the crystal to a rectangular shape. The (001) interface orientation is chosen because it has one of the lowest cleavage-free energies of all possible portlandite crystal faces [25][26][27][28]. Despite its hydrophilicity, the (001) crystal surface is more resistant to dissolution in water compared to other crystal faces of portlandite and brucite [29][30][31].
The created portlandite surfaces were brought into contact with a layer of NaNO 2 aqueous solution approximately 85 Å thick. According to experimental data [32], this size is comparable to typical capillary pores in hardened cement. After the insertion of the aqueous layer, the dimensions of the simulation supercell were~71.7 × 62.1 × 122.8 Å 3 . The number of H 2 O molecules in the nanopore corresponded to the density of the aqueous solution under ambient conditions (≈1 g/cm 3 ). The Na + and NO 2 − ions were initially uniformly distributed in the nanopore. The resulting nominal NaNO 2 solution molality in the pore space was ≈0.56 mol/kg (12,632 H 2 O molecules and 64 ion pairs). This corresponds to the experimentally measured solution molalities in the pores in cement materials [6]. Figure 1 shows a representative snapshot fragment of the simulation supercell. For the subsequent analysis, the interfacial space occupied by the solution was subdivided into three layers. The first near-surface layer (I) starts from the average positions of hydroxyl hydrogen atoms on the portlandite surface, H h , and has a thickness of 3.5 Å (approximately one molecular diameter of H 2 O. The second layer (II) immediately follows it and has the same thickness of 3.5 Å. In the rest of the pore space (III), the solution is considered to be much less affected by the surface and has properties close to those of the bulk phase.
To simulate the interface of portlandite with NaNO2 aqueous solution, a system co sisting of 11 × 10 × 8 portlandite unit cells was cleaved parallel to the (001) plane and c in the ab plane of the crystal to a rectangular shape. The (001) interface orientation is ch sen because it has one of the lowest cleavage-free energies of all possible portlandite cr tal faces [25][26][27][28]. Despite its hydrophilicity, the (001) crystal surface is more resistant dissolution in water compared to other crystal faces of portlandite and brucite [29][30][31].
The created portlandite surfaces were brought into contact with a layer of NaN aqueous solution approximately 85 Å thick. According to experimental data [32], this s is comparable to typical capillary pores in hardened cement. After the insertion of t aqueous layer, the dimensions of the simulation supercell were ~71.7 × 62.1 × 122.8 Å 3 . T number of H2O molecules in the nanopore corresponded to the density of the aqueo solution under ambient conditions (≈1 g/cm 3 ). The Na + and NO2 -ions were initially u formly distributed in the nanopore. The resulting nominal NaNO2 solution molality in t pore space was ≈0.56 mol/kg (12,632 H2O molecules and 64 ion pairs). This corresponds the experimentally measured solution molalities in the pores in cement materials [6]. Figure 1 shows a representative snapshot fragment of the simulation supercell. F the subsequent analysis, the interfacial space occupied by the solution was subdivid into three layers. The first near-surface layer (I) starts from the average positions of h droxyl hydrogen atoms on the portlandite surface, Hh, and has a thickness of 3.5 Å (a proximately one molecular diameter of H2O. The second layer (II) immediately follow and has the same thickness of 3.5 Å. In the rest of the pore space (III), the solution is co sidered to be much less affected by the surface and has properties close to those of t bulk phase.  The thick black lines schematically indicate the first 3.5 Å (I), second 3.5 Å (II) solution layers next to the surface, and the bulk phase (III). sym is the symmetry plane of this slit-like pore.

Force Field Parameters
A series of preliminary simulations for the bulk Ca(OH) 2 crystals was initially performed with both ClayFF-orig and ClayFF-MOH versions of the force field, in order to assess the relative accuracy of both versions. The MOH variation adds harmonic energy terms for Metal-O-H (M-O-H) angle bending in the crystal structure. The equilibrium angles and stiffness parameters were originally reported for Mg-O-H and Al-O-H angles [17]. They have been fitted to reproduce the lattice structure and dynamics obtained by density functional theory (DFT) calculations of brucite (Mg(OH) 2 )and gibbsite (Al(OH) 3 ). The optimized parameters for Mg-O-H angles are as follows [16]: θ 0,MgOH = 110 • and k MgOH = 6 kcal·mol −1 ·rad −2 . In this work, we use the same parameters for the Ca-O-H angles, based on the close similarity between the Mg(OH) 2 and Ca(OH) 2 crystal structures.
Water molecules were modeled using the flexible SPC/E model [33], and the parameters representing hydrated Na + ions were taken from the updated ClayFF model [16]. As mentioned above, we used here a new parameterization of charges, Lennard-Jones constants, and intramolecular constants for the NO 2 − ions in combination with the ClayFF force field [18]. This parameterization took into account recent results of DFT calculations for the nitrite ion dissolved in water [34,35], some earlier NO 2 − models use in classical MD simulations [36][37][38], and experimental data [39]. All force field parameters used in this work were specifically parametrized to be consistent with the same SPC/E water model, therefore they are expected to be transferable and provide reliable results in the simulations of other systems that are consistent with the same model of water, as has been already demonstrated numerous times in the ClayFF-based simulations of various nanoporous and layered materials and their aqueous interfaces over the last 20 years [16].

Simulation Details
All MD simulations were performed using the LAMMPS simulation package [40]. LAMMPS is one of the most commonly used computer simulation software packages for modeling of cementitious materials and their interfaces with aqueous solutions (e.g., ref. [11]). Periodic boundary conditions in all three dimensions were used in all simulations. A cutoff radius of 12.5 Å was used for short-range forces, and the PPPM method was used to account for long-range electrostatic interactions [41,42]. The GPU accelerator within the LAMMPS package was used to speed up the calculations [42,43]. Standard Lorentz-Berthelot mixing rules (ε ij = √ ε ii ε jj ; σ ij = (σ ii + σ jj )/2) were applied to calculate the Lennard-Jones parameters of interaction between unlike atoms [16,41]. The Newtonian equations of motion were numerically integrated with a time step of 0.5 fs via the velocity-Verlet algorithm [41,44].
The constructed models were initially equilibrated for 1 ns at the temperature of 298 K and pressure of 1 bar using the Nose-Hoover thermo-barostat [41,45,46]. After the equilibration in the NPT statistical ensemble (constant number of particles, constant temperature, and constant pressure), the equilibrium simulation run was performed in the NVT statistical ensemble (constant number of particles, constant temperature, and constant volume) using the Nose-Hoover thermostat [41,45,46]. The equilibrium NVT part was also 1.0 ns long and the generated equilibrium dynamic trajectories of atoms were then used for further statistical analysis.
The quantitative analysis of the structure, dynamics, and energetics of the portlanditesolution interfaces was performed following the standard approaches that are commonly used for such kinds of interfacial simulations (e.g., ref. [47]).  Table 1 shows the equilibrium parameters of the crystallographic unit cell of portlandite calculated with the ClayFF-orig and ClayFF-MOH versions of the force field, in comparison with experimental data [24,48], DFT calculations [26,49,50], and earlier results of classical MD simulations [11,[51][52][53][54]. In our simulations, all unit cell parameters were optimized in the NPT ensemble MD runs without any symmetry constraints. The results clearly reproduced the hexagonal crystal structure, in agreement with experimental data, and the discrepancy of lattice parameters did not exceed 2.4% for the ClayFF-MOH version. Both versions of the ClayFF force field also showed a good agreement with the DFT results. There were no significant differences in the crystallographic parameters between the ClayFF-orig and ClayFF-MOH models, indicating that introducing the explicit M-O-H angle bending terms in portlandite does not affect the symmetry of the crystal. One has to note, however, that this agreement is not particularly surprising because portlandite was one of the model structures used to parametrize the original ClayFF force field [20,21].

Results and Discussion
The small differences in the crystallographic parameters between the ClayFF-orig and ClayFF-MOH models are due to the much stronger localization of the hydrogen atoms and the orientation of hydroxyl groups for the latter case ( Figure 2). In the case of ClayFF-orig, the hydroxyl groups are freely oriented with their hydrogen atoms, H h , towards one of the three hydroxyl oxygen atoms, O h , of the opposite crystal layer, forming a strong hydrogen bond. Due to the unconstrained Ca-O-H angle in the ClayFF-orig, the hydroxyl groups can frequently and freely rotate, and those hydrogen bonds switch between three O h atoms from the opposite layer. Such poorly constrained rotations of hydroxyl groups lead to a slight expansion of the interlayer space. Consequently, the parameter c of the portlandite unit cell increased, and the parameters a and b decreased compared to the ClayFF-MOH. As    Figure 3 shows the distributions of angles of hydroxyl groups to the ab plane in the original and ClayFF-MOH versions, which also reflect a more disordered behavior of hydroxyls in ClayFF-orig. In contrast, the hydroxyl O-H bonds are more strongly localized around the normal to the ab plane in the ClayFF-MOH model. According to the DFT results, such localized behavior of hydroxyls is more realistic of brucite-like materials [17] and is in better agreement with experimental data [22].

Power Spectra of Atomic Vibrations in Portlandite Crystals
Analysis of atomic vibration spectra in the crystal gives further insights into the dy namics of portlandite hydroxyls. Figure 4 shows smoothed power spectra reflecting the vibrations in the ab crystallographic plane for all portlandite atoms and, separately, the contribution of Hh-atoms calculated in both versions of ClayFF. The smoothing of the spectra was carried out by applying the exponential damping with the characteristic time τdamp = 0.25 ps to the velocity autocorrelation functions, as proposed by Szczerba et al. [55] (a) (b)

Power Spectra of Atomic Vibrations in Portlandite Crystals
Analysis of atomic vibration spectra in the crystal gives further insights into the dynamics of portlandite hydroxyls. Figure 4 shows smoothed power spectra reflecting the vibrations in the ab crystallographic plane for all portlandite atoms and, separately, the contribution of H h -atoms calculated in both versions of ClayFF. The smoothing of the spectra was carried out by applying the exponential damping with the characteristic time τ damp = 0.25 ps to the velocity autocorrelation functions, as proposed by Szczerba et al. [55].

Power Spectra of Atomic Vibrations in Portlandite Crystals
Analysis of atomic vibration spectra in the crystal gives further insights into the dynamics of portlandite hydroxyls. Figure 4 shows smoothed power spectra reflecting the vibrations in the ab crystallographic plane for all portlandite atoms and, separately, the contribution of Hh-atoms calculated in both versions of ClayFF. The smoothing of the spectra was carried out by applying the exponential damping with the characteristic time τdamp = 0.25 ps to the velocity autocorrelation functions, as proposed by Szczerba et al. [55].   The low-frequency region of the spectrum corresponds to the librational modes of the hydroxyl groups and to the lattice vibrations. In the ClayFF-orig model, the individual peaks were unresolved and all merged into one single wide spectral band. In the ClayFF-MOH model, the individual peaks were better resolved for both libration and lattice vibration modes, which is in a better qualitative agreement with the vibrational behavior of portlandite crystal observed in experiments [56] and DFT calculations [50], compared to the ClayFF-orig model.
The most noticeable differences between the two versions of ClayFF in Figure 4a,b are in the high-frequency range of the spectra. The peaks corresponding to the O-H stretching modes of hydroxyls are located at positions 3624 cm -1 and 3676 cm -1 for the ClayFForig and ClayFF-MOH models, respectively. For the ClayFF-MOH, the peak frequency is slightly higher than the experimental data [56][57][58], where the O-H stretching peaks are located at 3620 cm -1 (A 1g mode) and 3646 cm -1 (A 2u mode) for Raman and IR spectroscopy, respectively. The lower frequency of the O-H stretching peak in the ClayFF-orig compared to the ClayFF-MOH indicates that the interlayer hydrogen bonding is stronger in the former model. As mentioned above, in the ClayFF-orig, H h atoms instantaneously bond with only one opposite-layer O h atom. In ClayFF-MOH, an H h atom can participate in up to 3 weaker hydrogen bonds at the same time, due to the higher localization of H h atoms in the structure. For instance, in brucite under normal conditions, the preferred orientation of the hydroxyl groups is along the normal to the magnesium octahedral layer [17,59]. The reorientation of hydroxyl groups in brucite and portlandite towards the oppositelayer O h atoms occurs only at high temperatures and pressures [23,60]. This exaggerated reorientation of hydroxyls under normal conditions in the ClayFF-orig model is especially noticeable in the isolated spectrum of H h -atoms in the ab crystallographic plane, where a very strong peak is observed in the high-frequency region, as seen in Figure 4b.

Structural Properties of Interfacial Solutions
First, we discuss the structural differences in the simulated aqueous solution between the two ClayFF versions. Figure 5 shows atomic density profiles of H w and O w atoms of the H 2 O molecules along the direction normal to the portlandite surface for both versions of ClayFF. As expected for the ClayFF-orig version, the profiles for H w and O w were similar to the earlier MD simulation results with this force field [20,51]. However, slightly different profiles are observed with the ClayFF-MOH model, indicating some changes in the local interfacial solution structure. Below, we argue that the reason for the changes in the density profiles in the first two solution layers, (I) and (II) (as defined in Figure 1) are directly related to very specific interfacial molecular coordinations.
bration modes, which is in a better qualitative agreement with the vibrational behavior portlandite crystal observed in experiments [56] and DFT calculations [50], compared the ClayFF-orig model.
The most noticeable differences between the two versions of ClayFF in Figure 4 are in the high-frequency range of the spectra. The peaks corresponding to the O stretching modes of hydroxyls are located at positions 3624 cm -1 and 3676 cm -1 for t ClayFF-orig and ClayFF-MOH models, respectively. For the ClayFF-MOH, the peak f quency is slightly higher than the experimental data [56][57][58], where the O-H stretchi peaks are located at 3620 cm -1 (A1g mode) and 3646 cm -1 (A2u mode) for Raman and spectroscopy, respectively. The lower frequency of the O-H stretching peak in the ClayF orig compared to the ClayFF-MOH indicates that the interlayer hydrogen bonding stronger in the former model. As mentioned above, in the ClayFF-orig, Hh atoms insta taneously bond with only one opposite-layer Oh atom. In ClayFF-MOH, an Hh atom c participate in up to 3 weaker hydrogen bonds at the same time, due to the higher loca zation of Hh atoms in the structure. For instance, in brucite under normal conditions, t preferred orientation of the hydroxyl groups is along the normal to the magnesium oc hedral layer [17,59]. The reorientation of hydroxyl groups in brucite and portlandite wards the opposite-layer Oh atoms occurs only at high temperatures and pressures [23,6 This exaggerated reorientation of hydroxyls under normal conditions in the ClayFF-o model is especially noticeable in the isolated spectrum of Hh-atoms in the ab crystal graphic plane, where a very strong peak is observed in the high-frequency region, as se in Figure 4b.

Structural Properties of Interfacial Solutions
First, we discuss the structural differences in the simulated aqueous solution betwe the two ClayFF versions. Figure 5 shows atomic density profiles of Hw and Ow atoms the H2O molecules along the direction normal to the portlandite surface for both versio of ClayFF. As expected for the ClayFF-orig version, the profiles for Hw and Ow were simi to the earlier MD simulation results with this force field [20,51]. However, slightly diff ent profiles are observed with the ClayFF-MOH model, indicating some changes in t local interfacial solution structure. Below, we argue that the reason for the changes in t density profiles in the first two solution layers, (I) and (II) (as defined in Figure 1) a directly related to very specific interfacial molecular coordinations.  [26,51]. With this orientation, the H 2 O molecule donates two hydrogen bonds to the O h atoms of the hydroxyls on the surface of the crystalline phase and this orientation of the H 2 O molecule corresponds to the bidentate coordination towards the surface [55]. In addition, there may be a monodentate coordination with the dipole vector oriented ≈90 • to the surface normal [55]. Such orientations of H 2 O molecules on the (001) surface of portlandite have been reported in other recent MD simulations [61]. Due to the more accurately constrained vertical orientations of the surface hydroxyl groups with the ClayFF-MOH model, H 2 O molecules in those two orientations are adsorbed closer to the portlandite surface. That, in turn, results in stronger hydrogen bonds in the donor-acceptor pairs H w ···O h and H h ···O w . The stronger bonding of H 2 O molecules in mono-and bidentate orientations is reflected in the shift of the first peak of the H w profile closer to the portlandite surface (in layer (I)) compared to the ClayFF-orig version.
However, the O w and H w interfacial density profiles for both versions of the force field clearly reproduce the distinct layering of water at least within the first two molecular layers, which is fully consistent with the results of previous simulations [20,51,61] and justifies our definition of the solution layers in Figure 1.
The Thus, the ClayFF-MOH model provides a better qualitative and quantitative description of the portlandite crystal and the portlandite-water interface. Therefore, in the rest of the paper, the behavior of dissolved ions and H 2 O molecules on the portlandite surface are presented and discussed only for the ClayFF-MOH model. Figure 6 shows the density profiles of Na + ions and N and O n atoms of NO 2 − in the sodium nitrite aqueous solution next to the portlandite surface. The density profile of Na + is similar to the previously reported profiles [20,51] with Na + is mainly adsorbed with inner-sphere surface coordination. Inner-sphere coordination is also observed for adsorbed and On atoms towards the Hh hydroxyls. In both orientations NO2 -ions accept hydrogen bonds from surface hydroxyls.  Figure 1).
Visual analysis of the trajectories with the OVITO package [62] has shown that the adsorbed NO2 -ions demonstrate four main types of surface coordination as illustrated in Figure 7a-d. In the first coordination (Figure 7a), all three atoms of the NO2 -ion are in a plane oriented almost parallel to the portlandite surface and coordinated by surface hydroxyls. The other three orientations have only the N atom, one On atom, or both On atoms coordinated by the surface hydroxyls, as shown in Figure 7b, 7c and 7d, respectively. These types of coordination can be seen for the anions adsorbed next to Na + ions on the surface as well as for the anions located far from adsorbed cations. As mentioned above, the first minima of the N and On density profiles are located within the solution layer (II), so that the density of NO2 -in that layer is lower than in layer (I) and in layer (III) that reflects bulk solution properties. This behavior of NO2 -in layer (II) is caused by the repulsion from the first NO2 -adsorption layer in layer (I). The surface distribution of molecules visualized by atomic density contour maps for the solution layer (I) shows that Na + ions are adsorbed at the centers of triangles formed by Oh atoms, as illustrated in Figure 8. The inner-sphere adsorption of Na + ions and the coordination of the H2O molecules around them by Ow atoms agrees well with the previous MD results [20,51]. Unlike Hu et al. [63], we did not observe the preference of NO2 -adsorption surrounding the cations and forming surface clusters. Instead, NO2 -ions can be adsorbed both near Na + and far from them. In both cases, the adsorbed anions are quite strongly bound to the surface and do not change their binding site within the time scale of the simulation. Nevertheless, the nitrite anions frequently change their orientation while staying at the same surface adsorption site. Due to the frequent changes of adsorbed NO2 -ions orientation, the H2O molecules that form hydrogen bonds with these NO2 -ions These types of coordination can be seen for the anions adsorbed next to Na + ions on the surface as well as for the anions located far from adsorbed cations. As mentioned above, the first minima of the N and O n density profiles are located within the solution layer (II), so that the density of NO 2 − in that layer is lower than in layer (I) and in layer (III) that reflects bulk solution properties. This behavior of NO 2 − in layer (II) is caused by the repulsion from the first NO 2 − adsorption layer in layer (I). and On atoms towards the Hh hydroxyls. In both orientations NO2 -ions accept hydrogen bonds from surface hydroxyls.  Figure 1).
Visual analysis of the trajectories with the OVITO package [62] has shown that the adsorbed NO2 -ions demonstrate four main types of surface coordination as illustrated in Figure 7a-d. In the first coordination (Figure 7a), all three atoms of the NO2 -ion are in a plane oriented almost parallel to the portlandite surface and coordinated by surface hydroxyls. The other three orientations have only the N atom, one On atom, or both On atoms coordinated by the surface hydroxyls, as shown in Figure 7b, 7c and 7d, respectively. These types of coordination can be seen for the anions adsorbed next to Na + ions on the surface as well as for the anions located far from adsorbed cations. As mentioned above, the first minima of the N and On density profiles are located within the solution layer (II), so that the density of NO2 -in that layer is lower than in layer (I) and in layer (III) that reflects bulk solution properties. This behavior of NO2 -in layer (II) is caused by the repulsion from the first NO2 -adsorption layer in layer (I). The surface distribution of molecules visualized by atomic density contour maps for the solution layer (I) shows that Na + ions are adsorbed at the centers of triangles formed by Oh atoms, as illustrated in Figure 8. The inner-sphere adsorption of Na + ions and the coordination of the H2O molecules around them by Ow atoms agrees well with the previous MD results [20,51]. Unlike Hu et al. [63], we did not observe the preference of NO2 -adsorption surrounding the cations and forming surface clusters. Instead, NO2 -ions can be adsorbed both near Na + and far from them. In both cases, the adsorbed anions are quite strongly bound to the surface and do not change their binding site within the time scale of the simulation. Nevertheless, the nitrite anions frequently change their orientation while staying at the same surface adsorption site. Due to the frequent changes of adsorbed NO2 -ions orientation, the H2O molecules that form hydrogen bonds with these NO2 -ions did not have any predominant coordination type, as seen in Figure 8b. The surface distribution of molecules visualized by atomic density contour maps for the solution layer (I) shows that Na + ions are adsorbed at the centers of triangles formed by O h atoms, as illustrated in Figure 8. The inner-sphere adsorption of Na + ions and the coordination of the H 2 O molecules around them by O w atoms agrees well with the previous MD results [20,51] ions did not have any predominant coordination type, as seen in Figure 8b.

Power Spectra of the Interfacial Dynamics
Power spectra reflecting the vibrations of Hh atoms at the surface and in the bulk structure of portlandite are shown in Figure 9a. As mentioned in Section 3.1.2, the lowfrequency region of the bulk portlandite spectrum has a band where several peaks can be resolved. In contrast to that, only two peaks can be resolved in that band for surface Hh, which indicates a superposition of normal vibration modes for surface hydroxyls. In the high-frequency range, the peak corresponding to the stretching of O-H bonds is broader, less intense, and significantly blue-shifted (3676 cm -1 vs. 3743 cm -1 ) for surface hydroxyls compared to the bulk phase, which is typical for surfaces of brucite-like minerals [59]. The shift of the O-H stretching vibrations peak for surface Hh at the crystal-solution interface is significantly smaller than with the crystal-vacuum interface [59]. The reduced peak shift is due to the formation of hydrogen bonds of the surface hydroxyls with the interfacial H2O molecules and NO2 -ions. Figure 9b shows that the power spectra of the H2O solution molecules at all distances from the surface are identical. This may merely reflect the fact that simple harmonic potentials used here to model the intramolecular vibrations in the SPC/E water model are not accurate enough, and more complex inharmonic models are necessary (see, e.g., ref. [55]). However, the libration, bending, and stretching modes of H2O molecules are reproduced with only small deviations from the experimental data [64,65].
Unlike the case with H2O molecules, the spectra reflecting the surface dynamics of Na + ions show some dependence on the distance from the surface, as seen in Figure 9c, which confirms the previous MD results [20], where such spectra have been analyzed in detail. However, the NO2 -power spectra in different solution layers have peaks at the same frequencies, differing only in magnitudes, as shown in Figure 9d. The positions of the spectral peaks are in good agreement with the results of Raman spectroscopy [39]: the calculated frequency of the bending mode, δ, is at 779 cm -1 , compared to 817 cm -1 from the Raman spectroscopy. The calculated asymmetric, νas, and symmetric, νs, N-On stretching modes peaks are at 1250 cm -1 and 1371 cm -1 , respectively, compared to the experimentally determined values of 1242 cm -1 and 1331 cm -1 .

Power Spectra of the Interfacial Dynamics
Power spectra reflecting the vibrations of H h atoms at the surface and in the bulk structure of portlandite are shown in Figure 9a. As mentioned in Section 3.1.2, the lowfrequency region of the bulk portlandite spectrum has a band where several peaks can be resolved. In contrast to that, only two peaks can be resolved in that band for surface H h , which indicates a superposition of normal vibration modes for surface hydroxyls. In the high-frequency range, the peak corresponding to the stretching of O-H bonds is broader, less intense, and significantly blue-shifted (3676 cm -1 vs. 3743 cm -1 ) for surface hydroxyls compared to the bulk phase, which is typical for surfaces of brucite-like minerals [59]. The shift of the O-H stretching vibrations peak for surface H h at the crystal-solution interface is significantly smaller than with the crystal-vacuum interface [59]. The reduced peak shift is due to the formation of hydrogen bonds of the surface hydroxyls with the Interfacial H 2 O molecules and NO 2 − ions. Figure 9b shows that the power spectra of the H 2 O solution molecules at all distances from the surface are identical. This may merely reflect the fact that simple harmonic potentials used here to model the intramolecular vibrations in the SPC/E water model are not accurate enough, and more complex inharmonic models are necessary (see, e.g., ref. [55]). However, the libration, bending, and stretching modes of H 2 O molecules are reproduced with only small deviations from the experimental data [64,65].
Unlike the case with H 2 O molecules, the spectra reflecting the surface dynamics of Na + ions show some dependence on the distance from the surface, as seen in Figure 9c, which confirms the previous MD results [20], where such spectra have been analyzed in detail. However, the NO 2 − power spectra in different solution layers have peaks at the same frequencies, differing only in magnitudes, as shown in Figure 9d. The positions of the spectral peaks are in good agreement with the results of Raman spectroscopy [39]: the calculated frequency of the bending mode, δ, is at 779 cm -1 , compared to 817 cm -1 from the Raman spectroscopy. The calculated asymmetric, ν as , and symmetric, ν s , N-O n stretching modes peaks are at 1250 cm -1 and 1371 cm -1 , respectively, compared to the experimentally determined values of 1242 cm -1 and 1331 cm -1 .   Figure 1).

Structure and Dynamics of the Interfacial Hydrogen Bonding
To better characterize the structural and dynamic properties of the interfacial solution, we have also studied the dynamics of hydrogen bonds (H-bonds), which were formally defined using a common geometric criterion: an H-bond between a hydrogen Hw (donor) and oxygen Ow (acceptor) atoms of two H2O molecules was assumed to exist if the Hw···Ow distance is less than 2.45 Å and the angle between the H-bond direction and the vector connecting the Ow atoms of the donor and acceptor molecules is less than 30° (e.g., ref. [66]). For H-bonds between H2O molecules and NO2 -(Hw···N and Hw···On pairs), the distance thresholds were taken according to [35]: RHw-N ≤ 2.25 Å and RHw-On ≤ 2.35 Å. The geometrical definitions for H-bonds between portlandite surface hydroxyls and H2O (or NO2 -) were taken to be the same as for H2O molecules (or H2O and NO2 -) in solution. The average lifetimes τHB for various types of H-bonds were then assessed by integrating the so-called continuous-time autocorrelation functions [35,[66][67][68] as illustrated in Figure  10. In addition, the average number of hydrogen bonds nHB per H2O molecule and acceptors N and On were also calculated.
It was found that the portlandite hydroxyls form the most stable H-bonds with the H2O molecules of the solution (pairs Hh···Ow and Hw···Oh). The lifetime of the Hh···Ow bonds is longer than that for Hw···Oh, despite the fact that in the ClayFF model, the charge of Oh is more negative than that of Ow, while the charges of Hh atoms and Hw atoms in SPC/E are almost equal. In the bidentate orientation of H2O molecules, which is most preferable on the surface of portlandite [26], H2O molecules are also coordinated by neighboring hydroxyls (Hh···Ow pairs) without forming hydrogen bonds between Hw and Oh of this hydroxyls.  Figure 1).

Structure and Dynamics of the Interfacial Hydrogen Bonding
To better characterize the structural and dynamic properties of the interfacial solution, we have also studied the dynamics of hydrogen bonds (H-bonds), which were formally defined using a common geometric criterion: an H-bond between a hydrogen H w (donor) and oxygen O w (acceptor) atoms of two H 2 O molecules was assumed to exist if the H w ···O w distance is less than 2.45 Å and the angle between the H-bond direction and the vector connecting the O w atoms of the donor and acceptor molecules is less than 30 • (e.g., ref. [66]  (III) denotes the bulk phase of the solution (see Figure 1).
As mentioned above, monodentate orientation can also occur, in which H-bonds between the Hh···Ow and Hw···Oh pairs are formed. Therefore, a significant role on the hydrated (001) portlandite surface is played by the H-bonds donated by the hydrogen atoms of the surface hydroxyl groups to the oxygen atoms of H2O molecules (pair Hh···Ow), which inhibits the penetration of H2O molecules closer to the portlandite surface and the dissolution of calcium ions [29]. Table 2 shows that the lifetime of H-bonds between H2O molecules in the solution layer (I) is longer than in the layer (II), while the number of H-bonds per H2O molecule remains about the same for both near-surface layers. This means that the interaction of H2O molecules with surface hydroxyls orients them in a way that favors the formation of a more stable H-bond network between themselves in layer (I). That H-bonding network of the H2O molecules nearest to the surface (in layer (I)) creates an ice-like twodimensional layer on the surface of portlandite [29]. In the bulk solution (layer (III)), the number of H-bonds per H2O molecule and the bond lifetime have typical values for lowconcentration aqueous solutions [69].   Figure 1).
As mentioned above, monodentate orientation can also occur, in which H-bonds between the H h ···O w and H w ···O h pairs are formed. Therefore, a significant role on the hydrated (001) portlandite surface is played by the H-bonds donated by the hydrogen atoms of the surface hydroxyl groups to the oxygen atoms of H 2 O molecules (pair H h ···O w ), which inhibits the penetration of H 2 O molecules closer to the portlandite surface and the dissolution of calcium ions [29]. Table 2 shows that the lifetime of H-bonds between H 2 O molecules in the solution layer (I) is longer than in the layer (II), while the number of H-bonds per H 2 O molecule remains about the same for both near-surface layers. This means that the interaction of H 2 O molecules with surface hydroxyls orients them in a way that favors the formation of a more stable H-bond network between themselves in layer (I). That H-bonding network of the H 2 O molecules nearest to the surface (in layer (I)) creates an ice-like two-dimensional layer on the surface of portlandite [29]. In the bulk solution (layer (III)), the number of H-bonds per H 2 O molecule and the bond lifetime have typical values for low-concentration aqueous solutions [69].
For interactions between water molecules and nitrite ions, the H-bonds H w ···O n are stronger than H w ···N, which is consistent with the DFT results for hydrated NO 2 − [35]. The H-bonds for these pairs have approximately the same lifetimes in all layers of the aqueous solution. However, the average number of H-bonds per acceptor atom in H w ···N pairs in layer (I) is lower than in layer (II), because some of the H 2 O molecules in layer (I) bind with surface hydroxyls instead of NO 2 − nitrogen atoms. The lifetime of H h ···N bonds between hydroxyls and NO 2 − is close to the lifetime of the H w ···N bond, however, the H h ···O n bonds have shorter lifetimes than Hw···On. Such a decrease in the lifetime is due to the vibrational behavior of the surface hydroxyls of the crystalline phase, which is also noticeable when comparing the lifetimes of the H h ···O w and H w ···O w pairs. We can also see that the H-bonds between surface hydroxyls and adsorbed NO 2 − ions are weaker than the ones between hydroxyls and H 2 O molecules. This is consistent with the absence of a strongly preferred NO 2 − orientation on the portlandite surface.

Water Molecules and Nitrite Ions Orientational Relaxation
The orientational relaxation of NO 2 − ions and H 2 O molecules can be characterized using the autocorrelation function of the unit vectors directed along the N-O n or O w -H w bonds [35,[70][71][72], which was calculated using the second-rank Legendre polynomial: where cos θ(t) = d T (t 0 )d(t 0 + t) is the linear correlation between the unit vector of bond orientation d at time moments t 0 and t 0 + t, and angular brackets . . . denote the averaging over all molecules and over all possible initial moments t 0 along the equilibrium MD trajectory. The correlation functions of NO 2 − and H 2 O molecules in three layers of the aqueous solution are shown in Figure 11. The orientational correlation times were calculated by fitting C 2 with bi-exponential functions [35]: where the shorter time constant τ 1 reflects the rotational inertia of an ion or molecule in the libration mode, which is affected by H-bonding around the ion/molecule [71]. The longer time constant τ 2 characterizes the reorientation of an ion or molecule due to the rotational Brownian motion and to large angular jumps associated with the breaking of H-bonds [35,71]. Table 3 shows the orientational relaxation times for NO 2 − and H 2 O, depending on the molecule location in the pore. For H 2 O molecules, the time τ 1 for all layers of the solution have similar values, indicating low sensitivity of librational molecular motions to the variation of local H-bonding environment between bulk solution and near-surface layers. The power spectra in Figure 9b also indicate the same behavior in the low-frequency region. On the other hand, the time τ 2 for H 2 O molecules increases in the near-surface layers. Such an increase in rotational relaxation time confirms the aforementioned formation of a more stable H-bonding network at the surface of crystalline portlandite due to hydroxyl group-water interactions.  Figure 1).
The increase in τ2 near the surface compared to their bulk solution values is seen for NO2 -ions as well, which is also due to the participation of the nitrite ions in the more stable H-bonding network at the surface of portlandite. In the closest near-surface layer, we also observe an increase in time τ1, indicating that the NO2 -ion adsorption affects its libration modes. The translational dynamics of NO2 -, Na + , and H2O molecules in the near-surface layers and in the bulk of the solution were studied by calculating the self-diffusion coefficients of these species using the Einstein-Smoluchowski relation [41,72,73]. For H2O molecules, we calculated the 3-dimensional (3d) diffusion coefficients separately for the solution layers (I), (II), and (III), as defined in Figure 1. For the ions, due to their relatively small number in the simulated systems, only bulk (layer (III)) 3-dimensional diffusion coefficients were calculated. Table 4 shows the calculated self-diffusion coefficients. The average self-diffusion coefficients for the ions are found to be = (0.84 ± 0.17) × 10 -5 cm 2 /s, = (1.46 ± 0.16) × 10 -5 cm 2 /s. The calculated diffusivities of ions in water are lower than the experimental values = 1.33 × 10 -5 cm 2 /s and = 1.91 × 10 -5 cm 2 /s [74]. The diffusional mobility of water molecules decreases as they approach the surface, with = (1.81 ± 0.16) × 10 -5 cm 2 /s, = (1.47 ± 0.18) × 10 -5 cm 2 /s, and = (1.36 ± 0.21) × 10 -5 cm 2 /s. The diffusivity of water in the bulk layer is underestimated compared to the experimental data = 2.33 × 10 -5 cm 2 /s [75]. The reduced calculated diffusivities of water and ions reflect the effect of their interactions with the portlandite surface but are also partly due to the known underestimation of the molecular mobility by the flexible SPC/E model and due to the finite size effects of the periodic boundary conditions [76]. (III) denotes the bulk phase of the solution (see Figure 1). The increase in τ 2 near the surface compared to their bulk solution values is seen for NO 2 − ions as well, which is also due to the participation of the nitrite ions in the more stable H-bonding network at the surface of portlandite. In the closest near-surface layer, we also observe an increase in time τ 1 , indicating that the NO 2 − ion adsorption affects its libration modes.

Diffusional Mobility of Water and Ions at the Portlandite Surface
The translational dynamics of NO 2 − , Na + , and H 2 O molecules in the near-surface layers and in the bulk of the solution were studied by calculating the self-diffusion coefficients of these species using the Einstein-Smoluchowski relation [41,72,73]. For H 2 O molecules, we calculated the 3-dimensional (3d) diffusion coefficients separately for the solution layers (I), (II), and (III), as defined in Figure 1. For the ions, due to their relatively small number in the simulated systems, only bulk (layer (III)) 3-dimensional diffusion coefficients were calculated. Table 4 shows the calculated self-diffusion coefficients. The average self-diffusion coefficients for the ions are found to be D Na = (0.84 ± 0.17) × 10 -5 cm 2 /s, D NO 2 = (1.46 ± 0.16) × 10 -5 cm 2 /s. The calculated diffusivities of ions in water are lower than the experimental values D exp Na = 1.33 × 10 -5 cm 2 /s and D exp NO 2 = 1.91 × 10 -5 cm 2 /s [74]. The diffusional mobility of water molecules decreases as they approach the surface, with D I I I H 2 O = (1.81 ± 0.16) × 10 -5 cm 2 /s, D I I H 2 O = (1.47 ± 0.18) × 10 -5 cm 2 /s, and D I H 2 O = (1.36 ± 0.21) × 10 -5 cm 2 /s. The diffusivity of water in the bulk layer is underestimated compared to the experimental data D exp H 2 O = 2.33 × 10 -5 cm 2 /s [75]. The reduced calculated diffusivities of water and ions reflect the effect of their interactions with the portlandite surface but are also partly due to the known underestimation of the molecular mobility by the flexible SPC/E model and due to the finite size effects of the periodic boundary conditions [76].

Conclusions
The recently modified ClayFF force field, ClayFF-MOH [10], with added Ca-O-H angular bending terms, together with a new parameterization for NO 2 − ions [18] were used to study the structure and dynamics of sodium nitrite solutions at the surface of portlandite using classical MD simulations. The (001) crystal surface of portlandite was selected for these simulations because it is the least reactive when interacting with water compared to other surfaces. We show that the inclusion of the Ca-O-H bending term does not affect the simulation results for the crystal structure but does affect the vibrational properties of the crystal, bringing them closer to the experimental data. This, in turn, affects the structural properties of the aqueous solutions in the near-surface layers due to the increased localization of surface hydroxyl groups of portlandite.
The better-constrained orientation of the hydroxyls on the (001) surface of portlandite leads to a more realistic description of the hydrogen bonds formed between surface hydrogen atoms and water molecules (H h ···O w donor-acceptor pair), facilitating the formation of a network of H-bonds between water molecules themselves, which form a surface layer inhibiting the dissolution of calcium ions into water. The formation of the interfacial Hbonding network manifests itself in the increase in the rotational relaxation time for water molecules near the portlandite surface.
Despite the formation of the water-water H-bonding network on the surface of the crystalline phase, the surface can adsorb dissolved nitrite and sodium ions. We show that the adsorption of nitrite ions involves the formation of H-bonds donated by surface hydroxyls not only to the oxygen atoms of NO 2 − , but also to the nitrogen atoms, i.e., both H h ···O n and H h ···N donor-acceptor pairs are present. Although the H-bonds H h ···O n and H h ···N are weaker than the H-bonds H h ···O w donated to the water molecules, NO 2 − ions infiltrate the H-bonding network on the surface of portlandite and are quite strongly adsorbed at it with no visible translational diffusion over the timescale of the entire MD simulation (~1 ns). On the other hand, the orientation of the adsorbed ions can change over time, so that there is no single preferred coordination of the adsorbed NO 2 − to the surface, meaning that the H-bonds in the H h ···O n and H h ···N pairs are less stable and can frequently break and re-form. The NO 2 − ions can be adsorbed onto the surface of portlandite next to the previously adsorbed Na + ions as well as separately without any strong association with the cations.