Impact of ethylene glycol on ions influencing corrosion in pores between iron oxide and calcium carbonate

ABSTRACT Injection wells for carbon capture and storage are made of iron casings supported by cement. Cement can transform into calcium carbonate and iron into iron oxide. Thus, these wells are subject to corrosion. Monoethylene glycol (MEG), used for hydrate prevention and natural gas dehydration, can adsorb both to calcium carbonate and iron oxide. As a first step towards an atomistic understanding of effects caused by adsorbed MEG, molecular dynamics simulations were applied to estimate free energy changes to hydronium and bicarbonate as MEG was adsorbed to the surface. We found that the global free energy minimum of both hydronium and bicarbonate was moved closer to the surfaces, due to adsorbed MEG, which may be caused by MEG replacing water molecules within the first water layers. This could have an effect on chemical reactions involving hydronium and bicarbonate. Minima in the free energy profiles other than the global one was found to originate from adsorbed water combined with interactions from the surfaces.


Introduction
Carbon capture and storage is one of several attractive methods of reducing the amount of carbon dioxide (CO 2 ) released into the atmosphere. Captured carbon dioxide is transported through pipelines to an appropriate location, which in many cases will be an underground saline aquifer [1] or possibly a methane hydrate reservoir [2][3][4]. Carbon dioxide sequestration to underground aquifers is a wellknown concept and has been commercially available since 1996 [5]. Appropriate saline aquifers will often be located close to abandoned oil and gas wells. Abandoned oil and gas wells consist of iron (Fe) casings both surrounded by and plugged by ordinary Portland cement [6]. However, due to dissolved sodium chloride (NaCl), as well as carbonic acid (H-2 CO 3 ) originating from carbon dioxide, iron will transform into iron oxide (Fe 2 O 3 ), while calcium oxide (CaO) of the cement transforms into calcium carbonate (CaCO 3 ). Thus, nano-to micro-scale pores between the casing and the cement make up ideal environments for corrosion, making the abandoned oil and gas wells a potential source of carbon dioxide leakage [7]. Moreover, operational carbon dioxide injection wells will also contain a similar iron casing, which is supported by cement, thus yielding another potential source of carbon dioxide leakage.
Under appropriate conditions, such as in an underground aquifer or oil reservoir, hydrate (with carbon dioxide as guest molecules) may form [8]. In these cases, hydrate inhibitors must be injected into the system. A common hydrate inhibitor is monoethylene glycol (MEG) [9], which is also used for natural gas dehydration and may thus leave traces of MEG in the dried gas. In some cases MEG is also injected directly into gas processing plants at the critical points of hydrate formation.
Experiments have shown that MEG acts as a corrosion inhibitor of iron in carbon dioxide rich environments. This effect has been attributed to the reduced solubility and diffusivity of carbon dioxide, as well as reduced water activity and solution polarity, due to presence of MEG [10] and is therefore not related to the effects of adsorbed MEG. Thus, for trace amounts of MEG already adsorbed to the iron surface, this effect would not apply. To our knowledge, no investigations have been performed to understand how adsorbed MEG may affect corrosion processes on the surfaces of the pores at the atomistic level.
Carbonic acid will react with water to form hydronium ions (H 3 O + ) and bicarbonate ions (HCO − 3 ), which influence and cause corrosive reactions with iron oxide and calcium carbonate, respectively. A corrosive reaction on calcite can occur via CaCO 3 + H 3 O + Ca 2+ + HCO − 3 + H 2 O, where the calcium ions dissociate from the surface. On iron, the corrosion process, which forms a hematite surface, is driven by electron transport and will be accelerated in the presence of ions, such as sodium, chloride, hydronium, and bicarbonate. As a first step towards understanding the interaction mechanisms that may arise due to adsorbed MEG, we used molecular dynamics simulations to investigate how adsorbed ethylene glycol molecules on both calcite, as well as iron oxide (in its hematite form) affected free energy profiles of hydronium and bicarbonate. Changes in the ions free energy landscape can affect their ability to reach the surfaces and may affect the probability of reactions with calcite, as well as their effect as ions that accelerate corrosion of iron, thus changing the corrosion rate. For example, a large reactant free energy barrier outside the surface, much higher than k B T, where k B is the Boltzmann constant and T the temperature, would prevent any influence from the ions altogether, since the reactant would never reach the surface. Hence, we can utilise molecular dynamics, without reactive force fields, to learn more about the processes that take place before corrosive reactions.
In Section 2, we give a detailed description of our computational experiments. In Section 3, we present our results. Our findings are summarised in Section 4.

Software packages
To perform molecular dynamics (MD) simulations we used HOOMD-blue version 1.3.1 including the particle-particle particle-mesh (P 3 M) algorithms [11][12][13][14]. The HOOMD-blue code was extended to support sampling of potential of mean force according to the adaptive biasing force (ABF) scheme (see Supplementary Information), as well as to support both the Buckingham potentials and the 12-8 Lennard-Jones potentials without resorting to tabulated potentials. The 12-8 Lennard-Jones potential was used in combination with the Buckingham potential to create a full potential consisting of r −6 , r −8 and r −12 terms, as well as the exponential term, where the 6th and 8th order terms correspond to the dipoledipole and dipole-quadrupole van der Waals terms. GDIS version 0.99 was used to construct crystal structures based on existing crystallographic information files (CIF) [15]. Visual Molecular Dynamics (VMD) V1.9.1 was used to view and visualise MD trajectories [16], while we used PackMol, in addition to in-house developed software, to create initial system configurations [17]. The Quantum Espresso plane wave package [18,19] (PW) version 6.3 was used to perform quantum mechanical electronic structure calculations. Further details on simulation setup and parameters used within some of these software packages are given in Section 2.5. To calculate solvation free energies the 9 Dec 2014 version of LAMMPS [20,21] with the USER-FEP package was employed. Figure 1 shows equilibrated starting configurations of the four molecular dynamics systems simulated within this work. In the first studied system (Figure 1(a)) we created a {10 14} surface of calcium carbonate in its calcite form. A 24.94 by 24.882 by 34.944 Å calcite slab was constructed based on crystallographic information (at 297 K) found in the literature [22], where the z-axis was positioned perpendicular to the {10 14} surface. Immediately outside the surface we placed a 40 Å long aqueous phase containing 790 water molecules, yielding a density of approximately 1 g/cm 3 at 298 K. Inside the aqueous phase we placed 16 MEG molecules, which were adsorbed through equilibration simulation runs, as well as a single hydronium ion to be controlled by the adaptive biasing force (ABF) sampling scheme (see Supplementary Information). The second system to be studied is shown in Figure 1(b), where we replaced the hydronium ion with a bicarbonate ion. Molecules in the bulk phase were placed in an initial configuration (prior to equilibration) using PackMol to achieve an appropriate initial packing and random distribution of the molecules. In all cases, water was specified to be initially distributed within boxes placed 2 Å outside the bounding-boxes of the surfaces (in the z-direction), thus providing space for the first water layers to form once the simulations were started. MEG was initially distributed within boxes placed 1.5 Å outside the bounding-boxes of the surfaces, stretching to 3.5 Å outside the surfaces. MEG was placed closer to the surfaces than water to ensure proper adsorption of the MEG molecules before execution of the free energy profile simulations for hydronium and bicarbonate.

Model description
It is known from experiments that the {001} hematite surface often occur as a fully hydroxylated surface at low water vapour pressure [23], thus making the hydroxylated surfaces the most commonly studied. However, research has also shown that this surface can be slightly hydrophobic at natural pH levels [24], hence prompting us to study the unhydroxylated surface, which is an important one to study. The 30.214 by 29.08 by 26.892 Å hematite slab was constructed based on crystallographic information found in the literature [25] (where it is referred to as 'crystal 1') and was followed by a 40 Å aqueous phase of 1117 water molecules containing 16 adsorbed MEG molecules and a single hydronium ion (see Figure 1(c)). To study adsorption of bicarbonate we replaced hydronium by a bicarbonate ion (see Figure 1(d)).
To provide a clearer picture of the complete system that is studied, which consists of an aqueous pore between calcite (cement) and hematite (rust), the snapshots in Figure 1 are all drawn with the calcite slab to the left and the hematite slab to the right. Similarly, all of the following plots are drawn with the calcite surface to the left, at z = 0, and the hematite surface to the right, at z = 0, thus yielding a positive z-axis for the plots involving calcite and a negative z-axis for the plots involving hematite.
In all the systems we aimed to study the effect of adsorbed MEG. Thus, we also needed to construct all aforementioned systems without adsorbed MEG, where we wanted to achieve a constant density of the bulk phase, close to 1 g/cm 3 . A bulk phase of 40 by 24 by 24 Å yields a theroretical density of approximately 1.02 g/cm 3 (where the molar masses were m H = 1 g/mol, m O = 15.9994 g/mol and m C = 12.0107 g/mol). Keeping a fixed volume but adding 16 MEG molecules yields a theoretical density of approximately 1, 096 g/cm 3 . Hence, the uncertainties in density due to the gap between the surface and the bulk, as well as density fluctuations near the surfaces, are larger than the differences in density due to added MEG (which will also eventually adsorb to the surfaces). Therefore, we did not attempt to adjust the number of water molecules after adding MEG, to reduce the density. Systems were periodic in all directions.

Force field description
For the calcite crystal we employed the force field parameters of Xiao et al. [26], with long-range forces modelled using Coulomb potentials of the form where r ij is the distance between atoms i and j, e is the elementary electron charge, q i , is a partial charge in units of |e| and e 0 is the vacuum permittivity. Short-range potentials were described using Lennard-Jones 12-6 potentials of the form where e ij and s ij are fitting parameters of the model. Angular forces of the oxygen-carbon-oxygen angles of the carbonate were represented by harmonic potentials where u ijk is the angle between atoms i, j and k, u 0 ijk is the equilibrium angle and K u ijk is a fitting parameter of the model. To keep the carbon atoms of the carbonates in the same plane as the oxygen atoms of the carbonate a dihedral of the form was used, where f ijkl is the torsional angle between the planes i −j−k and k−j−l and K f,h ijkl is the fitting parameter. Short-range interactions between calcite and water were taken from the force field of Xiao et al. [26], while short-range interactions between calcite, bicarbonate ions, hydronium ions and MEG where calculated using geometric combination rules (s ij = s i s j √ and e ij = e i e j √ ) based on pure parameters for calcite from the same force field. For hematite we used the force field parameters from Kerisit [23], where we chose the force field parameters based on a modified CLAYFF force field.
As is the case for the standard CLAYFF model [27], all interatomic interactions in the crystal are described by Equations (1) and (2). Short-range interactions between hematite and water were primarily modelled using the Lennard-Jones parameters of Kerisit [23], except for interactions between oxygen of water and iron, where we used the potential More information about the form of this potential can be found in the literature [28,29]. Pure parameters of the hematite model were used with geometric combination rules to calculate short-range interactions between hematite, bicarbonate ions, hydronium ions and MEG. We used hydronium parameters from Jang et al. [30], which in-turn are based on the DREIDING force-field [31], where long and short-range forces were modelled using Equations (1) and (2), respectively. Both bonds and angles were simulated using harmonic potentials and the potential of Equation (3), respectively. Short-range interactions between hydronium ions, water, bicarbonate ions and MEG were estimated using geometric combination rules.
Interactions between bicarbonate ions and water were parameterised using the force field of Zeebe [32], where long and short-range potentials are described by Equations (1) and (2). For this force field the parameters e Ow−X and e Ow were known, where Ow is oxygen of water and X is some atom in bicarbonate. This enabled us to calculate pure parameters e X , and similarly s X using geometric combination rules. The pure parameters were used between bicarbonate molecules, as well as between bicarbonate and MEG.
The original bicarbonate force field is rigid [32], but was made flexible using bonds, angles, dihedrals and out-of-plane bending parameters from Demichelis et al. [33]. Making the molecule flexible was done to achieve more realistic mechanical behaviour when moving between narrow cavities formed by adsorbed MEG.
Bonds, angles and dihedrals were modelled using the potentials of Equations (5), (3) and (4), respectively. To incorporate the out of plane bending potential [33], we were required to convert an out-of-plane bending potential of the form U a = k 1 d 2 + k 2 d 4 (where d is the distance to the plane) to a potential of the form U b = k(1 − cos f). To achieve this, we first note that U a ≈ k 1 d 2 since d ≪ 1. Taylor expanding U b (f(d)) to second order in d around d = 0 we found a coefficient that could be calculated with respect to k 1 , which was given in the force field from Demichelis et al. [33] (also see Supplementary information).
All force field parameters used are listed in Supplementary information. Short and long-range potentials were turned off between atoms connected by one bond (1-2 interactions), atoms connected by two bonds (1-3 interactions) and atoms connected by three bonds (1-4 interactions), except when modelling MEG, where 1-4 interactions were scaled by 0.5. Reinstating part of the 1-4 interactions was achieved by extending the HOOMD-blue software package to include a special type of Lennard-Jones plus Coulomb bond type.
Our choice of forcefields have been justified in Supplementary Information.

Hematite surface depolarisation
Once the hematite surface was cut it showed a large dipole moment, which caused a long range tail in the free energy profiles. This could be observed for the simulations performed for low sampling of the ABF simulations (free energy plots shown in Supplementary information) before adjustments were made to depolarise the surface. To achieve a dipole moment close to zero, we decreased the outermost positive iron (Fes) charges and increased the outermost oxygen (Ors) charges until we reached a net dipole moment of near zero. The new partial charges were 0.5025641 for outermost iron and −0.70145833 for outermost oxygen.
A second method of providing a hematite surface with vanishing dipole moment was also applied, where the original surface charges were applied and depolarisation was achieved by moving ions from one side of the surface to an equivalent site on the opposite side. This is referred to as reconstruction of the surface. If nothing is stated, the free energy profiles are generated using hematite slabs with surface charge scaling, while other plots are generate using an unmodified surface. If reconstructed hematite slabs were used, this is stated explicitly. The dipole moment was reduced to 0.126 Å using reconstruction.

Simulation details
Eight systems were simulated, where variations of those simulations are listed in Table 1. Four systems contained adsorbed MEG and four systems did not. Each of the systems were first equilibrated, where initially all particle velocities were assigned randomly from a Gaussian distribution with a vanishing mean value and standard deviation v, where v is the solution to mv 2 = k B T (where k B is the Boltzmann constant). Subsequently, the average momentum was subtracted from each particle to achieve a vanishing total momentum. Figure 1 shows the equilibrated systems containing adsorbed MEG. All simulations were performed in the NVT ensemble at 298 K, where temperature was controlled using the Nosé-Hoover thermostat [48] with a relaxation time of 100 times a single time step. P 3 M Ewald summation was used to handle electrostatic forces across periodic boundaries, where we adjusted the grid spacing to approximately 2 Å in all directions applying a second order interpolation between the grid points. All cutoff radii were set to 12 Å.
After equilibration, free energy profiles were estimated using the ABF method on hydronium and bicarbonate ions of the systems, where a force constant of 1000 kJ·mol −1 Å −2 was used for the moving windows. We applied seven windows over a range of 20 Å outside the surfaces, where each window was sampled for 10 ns with 1 fs time step. Within each window a value N fs = 500 was applied (see Supplementary Information) to assure a proper average for the biasing force calculations. The range of 20 Å was divided into 400 bins. All of the eight ABF simulations thus had a 70 ns timespan. A second set of ABF simulations were also performed to investigate how well the systems converged. In that case, a force constant of 200 kJ·mol −1 Å −2 was applied and N fs was set to 1000. Moreover, the simulations were run with 20 windows over a range of 20 Å, where each window was sampled for 20 and 17 ns for calcite and hematite systems, respectively.
For the second set of simulations, we included two ABF molecules to double the number of samples, as well as to ensure a larger coverage of the (x, y)-plane. Each molecule was equipped with their own biasing force, thus resulting in two independent ABF simulations that were combined by cumulatively adding the bins of average forces from each ABF molecule into a single set of bins of average forces. In cases where the interaction between the two molecules were strong enough to significantly disturb the resulting free energy profiles, this would show up as one unnatural peak within each window of the free energy profile. Each such peak originated from each ABF molecule having to overcome the cluster of all other ABF molecules within the current window. For the calcite systems, such behaviour was only observed if 10 or more ABF molecules were included within the simulation and it was not observed in the case of two ABF molecules.
As for the first set of ABF simulations, the number of bins, over 20 Å, was set to 400, thus yielding 20 bins per window. The second set of ABF simulations will be referred to as the 'high sampling' scheme and the first set of simulations as the 'low sampling scheme'. It should be noted that the first and the second set of simulations were started with random positioning of MEG, hydronium and bicarbonate molecules in the (x, y)-plane.
To estimate solvation free energies of hydronium, bicarbonate and MEG, three systems, consisting of 1000 water molecules and the three molecules, respectively, were constructed, using the same force fields as for the systems described above. An initial cubic system with sides of 32 Å were used. First 20,000 steps of 0.1 fs were executed using a Langevin thermostat (NVE ensemble), followed by a 250,000 step equilibration run with 1 fs time steps. Each production run was set to 8·10 6 steps of 1 fs. Temperature and pressure were set to 298 K and 1 bar, respectively.
To arrive at the solvation free energies, long-and short-range interactions between water and the respective molecules were scaled by a factor λ, where λ was varied from 0 to 1 throughout the production run. The production run was divided into 100, where λ was increased evenly for each 80,000 step. At the end of each 80,000 step, states at 16,000 steps separated 5 steps apart, were used to calculate 〈U(l i+1 ) − U(l i )〉, 〈V exp { − [U(l i+1 ) − U(l i )]/[k B T]}〉 and 〈V〉, where U is the potential, V is the system volume, T is the temperature and k B is the Boltzmann constant. These values were then used to calculate the free energies using thermodynamic integration (TI) and free energy perturbation (FEP). All other simulation details were identical to those used by Olsen et al. [49], where a more detailed description of the calculations are also given. The TI simulations are listed in Table 1.

Results and discussions
3.1. Convergence of ABF simulations Figure 2(a-d) shows the average forces sampled over 400 bins for both sets of ABF simulations (both for the low sampling scheme and for the high sampling scheme). The best convergence could be seen for hydronium on calcite, where the results of high and low sampling schemes overlap with a maximum error at the first peak of approximately 20% (with near overlap of the remaining peaks) and a deviation along z of approximately 0.3 Å. Convergence was more difficult to achieve for bicarbonate since this molecule has a more complex shape and requires more samples than the simpler hydronium molecule. However, in both cases it can be observed that the same basic features are present both at high and low sampling schemes. That is, we can observe the same number of maxima and minima and the roots of the plots occur at approximately the same places along the z-axis. Hence, there will be a limitation related to comparing free energy values resulting from these curves, unless the value differences are qualitatively different (e.g. a sign shift or similar). Systems with adsorbed MEG experienced the highest convergence problems with up to 0.5 Å displacements along the z-axis and additional features in the force profiles.
Free energy is obtained by integration of the mean force profiles shown in Figure 2(a-d). Therefore, vanishing forces indicate local or global extremums of the free energy. It is the rate of change of the force that governs the magnitude of the free energy extremums. Since the force profiles, to a large extent, retain transitions from positive to negative values, the changes of the free energy profiles remain similar, despite changes made to the parameterisation. However, the absolute adsorption free energies are expected to be relatively dependent on parameterisation. This can also be seen by comparison with free energy profiles obtained using the low sampling scheme, shown in Supplementary information. Hence, changes in the free energy are likely to be more reliable than absolute free energies.

Solvation free energies
Using combination rules to obtain intermolecular force field parameters can only provide a rough estimate of the real potential energy surface. To provide an idea of how well the thermodynamics are reproduced by the chosen force fields, solvation free energies were estimated for hydronium, bicarbonate and MEG. The thermodynamic integration curves, 〈∂U(l)/∂l〉, obtained from the simulations, together with the free energy perturbation steps, are shown in Figure 3.
The solvation free energy of MEG was found to be −28 + 2.2 and −31 + 2.1 kJ/mol using TI and FEP, respectively, where the errors are based on statistical fluctuation that were removed by smoothing, resulting in the curves of Figure  3. These values correspond well with the ones obtained by Olsen et al. [49], with a different force field for water.
For hydronium the solvation free energies were estimated to −168 + 2.2 kJ/mol and −172 + 2.2 kJ/mol, using TI and FEP, respectively. Reported experimental results yield a solvation free energy of −405 kJ/mol [50], where we needed to convert from conventional free energies, using the expression DG = DG conv. − 1057 kJ/mol [51].
For bicarbonate solvation free energies of −221 + 2.6 and −227 + 2.6 kJ/mol were found using TI and FEP, respectively. Reported experimental results have measured a solvation free energy of −335 kJ/mol.  ) and with (c, d) adsorbed MEG, and in vacuum without adsorbed MEG (i, k). Free energy profiles of hydronium (e, f) and bicarbonate (g, h) in the vicinity of calcite and hematite surfaces in pure aqueous environment with and without adsorbed MEG, as well as in vacuum (i-l). Calcium at calcite and iron at hematite are located at approximately z = 0, respectively. Plots coloured in grey were produced with low sampling, while plots in blue were produced with high sampling (longer simulations per ABF step, more ABF steps and two simultaneous ABF molecules).

Hydronium and bicarbonate free energies
In the following, we define the magnitude of free energy of adsorption as |A bound − A unbound |, the desorption free energy barrier as A barrier − A bound and the adsorption free energy barrier as A barrier − A unbound . The definitions of the three energy barriers are illustrated in Supplementary information. Figure 2(e-h) shows the free energy profiles of hydronium and bicarbonate in the vicinity of the {10 14} calcite surface, as well as in the vicinity of the {0001} hematite surface.
In case of MEG adsorbed to calcite the magnitude of free energy of adsorption for hydronium varied by approximately 5 kJ/mol, which based on the discussion of convergence can be considered to be insignificant. Successive minima in the free energy curves can be observed both with and without adsorbed MEG. However, the free energy curves in the presence of adsorbed MEG are displaced closer to the calcite surface by 0.8 Å. Such a displacement towards the surface could also be observed in the low sampling scheme (see Supplementary information).
The free energy profile of bicarbonate on calcite experienced a doubling of the desorption free energy barrier, accompanied by a displacement of the free energy minimum towards the surface, due to adsorbed MEG. Thus, adsorbed MEG would enhance any reaction involving bicarbonate ions with calcite. Changes of the free energy profile of bicarbonate in the vicinity of hematite, due to adsorbed MEG, were of similar proportions as hydronium near calcite, as well as in the same direction, thus being insignificant compared to the accuracies of the simulations.
For hydronium in the vicinity of hematite, a three fold increase of the magnitude of free energy of adsorption was observed, due to adsorbed MEG, as well as a displacement of the free energy profile in the direction of the surface.
Within the accuracies of the computational procedures that were performed, the free energy profiles of both hydronium and bicarbonate either experienced a displacement towards the surface or remained unchanged, regardless of the chosen surface. Moreover, it could be observed that the magnitude of free energy of adsorption either increased or remained approximately the same, due to adsorbed MEG. Successive minima in the free energy profiles close to the free energy global minima were present, regardless if adsorbed MEG was added to the surface or not. It could also be observed that the free energy mnima remained approximately the same both with and without adsorbed MEG, except for being displaced in the same direction as the free energy minima. Figure 4 shows the free energy profiles for hydronium and bicarbonate in the vicinity of the reconstructed hematite slab. Comparing with Figure 2(f,h), it can be seen that in both cases MEG results in hydronium and bicarbonate free energy minima closer to the surface. However, in the case of using the reconstructed hematite slab, this effect is stronger. Moreover, due to the higher surface charges, in the case of reconstruction, the adsorption free energy magnitudes were higher.  Figure 5 shows the adsorption geometry of MEG and its distribution at the two surfaces after equilibration of the systems. At calcite we could observe a tendency towards hydrogen atoms of MEG adsorbing to oxygen, while oxygen atoms of MEG adsorbed to calcium. This is not unnatural due to the sign of the partial charges involved. At hematite, hydrogen atoms of MEG preferred to adsorb to oxygen, while the oxygen atoms of MEG preferred to adsorb to the positively charged iron atoms.

Adsorption geometry of MEG
Full adsorption (i.e. all 16 MEG molecules of the simulation) on the surfaces seen in Figure 5 corresponds to only 0.027 and 0.019 g/cm 2 surface coverage of MEG adsorbed to calcite and hematite, respectively.
In Figure 6 the density profiles of adsorbed MEG were plotted, based on the ABF simulations of hydronium on calcite and hematite. However, the first 35 ns were excluded to avoid interference from the presence of hydronium in the vicinity of the region where density was to be estimated (we know that during the three last ABF windows, hydronium would not be in this region).
In addition to plotting the density based on all atoms of MEG, the densities for hydroxyl groups of MEG and ethylene groups (i.e. remaining carbon and hydrogen atoms) of MEG were plotted separately in the same figures.
At calcite, hydroxyl groups of MEG had density peaks both at 2 and at 3.5 Å with approximately the same density.
Moreover, the ethylene groups had a density peak at approximately the same location as the second hydroxyl density peak. This implies that an adsorption geometry consisted of one adsorbed hydroxyl, while the rest of the ethylene group was horizontally oriented with respect to the surface, thus leaving the other hydroxyl group at the same position along the zaxis as the ethylene group or adsorbed at the same distance from the surface as the first hydroxyl group.
The results from MEG adsorption on calcite can be compared with adsorption of ethanol (C 2 H 5 OH), which has a similar structure as that of MEG. Molecular dynamics simulations, together with experimental observations, for ethanol   adsorption to calcite is reported by Cooke et al. [52]. There, it was observed that the OH-group of ethanol oriented towards the calcite surface, as was also the case for MEG. In the case of OH of MEG we observed a density peak approximately 2 Å from the surface, while for ethanol it was reported to be located 2.2 Å from the surface. For ethanol, the centre of mass was reported to be located 3.2 Å away from the surface, while for MEG a centre of mass density peak was found at approximately 3.5 Å. The desorption free energy barrier for ethanol was reported to be approximately 23 kJ/mol. As can be seen from the comparison with reported results, the adsorption behaviour observed for one of the MEG adsorption modes is similar to that observed for ethanol.
A different adsorption behaviour of MEG could be observed at hematite (compared to adsorption at calcite). Density peaks at hematite were in general closer to the surface, thus indicating a stronger adsorption. Furthermore, at hematite, for the majority of adsorbed MEG, both hydroxyl groups attached to the surface, but at two different locations relative to the surface (at 1.7 and 1.3 Å from the surface), while the ethylene groups of MEG pointed towards bulk.
In Figure 6, density profiles produced for hematite slabs that were reconstructed to achieve a vanishing dipole moment are also shown. It can be seen that the adsorption behaviour on the two types of slabs are qualitatively similar. For the reconstructed surface, two peaks are present at approximately 4 Å from the surface, one for the hydroxyl group and one for the ethylene groups. Visual inspections of the molecular trajectories revealed that these originate from MEG molecules where one hydroxyl group points towards the surface. The ethylene group of these molecules point towards bulk and the other hydroxyl group is bent away such that it lines up besides the outermost carbon atom (along z). Hence, the ethylene group peak at 4 Å originates from the outermost parts of the MEG ethylene groups. Figure 7 shows pair correlation functions of MEG, pair correlation functions of bicarbonate without adsorbed MEG, and pair correlation functions of bicarbonate with adsorbed MEG. For each surface, only pair correlations with peaks closer than 2.5 Å from the surfaces were included in the plots. I.e. if a pair correlation had a peak close enough only with adsorbed MEG, then it was included also for the case without MEG and vice versa.
Peaks further away from the surfaces were considered to have relatively weak interactions with the surfaces. It can be seen that, in all cases, hydrogen and oxygen atoms of both MEG and bicarbonate appear at similar distances from the surfaces, implying that there may be a competition between the two types of molecules for adsorption. In the case of adsorbed MEG, the hydrogen atom of bicarbonate oriented itself significantly closer to the calcite surface carbon atoms. Adsorbed MEG also resulted in the hydroxyl oxygen of bicarbonate (Ob1) orienting itself closer to surface calcium atoms. Furthermore, a small peak for bicarbonate carbon atoms at a distance of approximately 2.5 Å from the surface oxygen atoms could be observed to form due to adsorbed MEG. These are all effects that support a stronger attraction between bicarbonate and calcite when MEG is adsorbed to the surface. The situation is different at hematite, in which case the pair correlation peak, at approximately 2.2 Å, between bicarbonate hydroxyl oxygen and surface iron atoms vanishes, due to adsorbed MEG. It can also be observed that the first peak between bicarbonate carbon and surface oxygen atoms was displaced in the direction of the surface, due to adsorbed MEG.

Density profiles of water
To investigate the impact of adsorbed MEG on adsorbed water the density profiles of water were plotted, with and without adsorbed MEG, where both hydrogen and oxygen of water were included when estimating the densities. The density profiles were generated from an 8.5 Å region outside each surface, which was divided into 100 bins. As basis we used the last 35 ns of trajectories from ABF simulations of hydronium. The results are shown in Figure 8.
As can be seen, the presence of even small amounts of adsorbed MEG on calcite resulted in a significant attenuation of the first three water layers, effectively displacing water into the fifth water layer and bulk.
On hematite, no attenuation of the first water density peak was observed since water molecules forming this density peak resided closer to the surface than the larger MEG molecules. However, a significant displacement of water away from the second density peak into the third and fourth density peaks and bulk could be seen. On Figure 8, we also plotted the density profile for the hematite slab that was reconstructed to achieve a vanishing dipole moment. It can be seen that the density peak close to the surface was lowered, hence indicating a weaker water layer close to the surface. It can also be observed that other peaks were moved closer to the surface peak.

Origins of free energy changes
It can be observed that for all free energy profiles of Figure 2(e, f) involving both water and hydronium, there were non-zero free energy minima other than the global one. The free energy profile for hydronium, with and without adsorbed MEG, had similarities with regards to shape and free energy minima. This implies that the origin of the free energy minima outside the global one is other than from adsorbed MEG. Thus, the origin of the barrier must have been from water or the surface, where the effect from passing through the water layers is well known from previous literature and should be a primary contributor to the free energy minima. Figure 2(j,l) shows that all free energy profiles of hydronium and bicarbonate in vacuum lack the aforementioned barriers. However, a closer look at the underlying mean force profiles, shown for hydronium and bicarbonate near calcite in Figure  2(i,k), show two distinct peaks in the forces that manifest themselves as plateaus with non-vanishing derivatives in the free energy profiles. The forces seen in these plots become reshaped once water is introduced, thus implying that both water and the surface may contribute to the observed minima.
It should be noted that correlations between minima in the free energy and water density profiles of Figure 8 do not correlate perfectly, which is not to be expected, since both hydronium and bicarbonate have different spatial shapes and will therefore be affected differently by the layers of water close to the surface. Moreover, the different layers of water are defined by distinct orientations of the water molecules, resulting in distinct dipole moments that will have different effects on the two ions. Figure 9 shows a correlation between edges of the water density peaks (where the derivatives of the density profile are largest between a maximum and minimum) and the free energy profiles of hydronium near calcite. Positive density derivatives are shown in red and negative ones in blue. It can be seen that the range of the water layers correlate well with the range of the free energy minima.
Comparing free energies of hydronium and bicarbonate in vacuum of Figure 2(j,l) with those in water without adsorbed MEG from Figure 2(e,g), it can be seen that the free energy profiles in vacuum have much larger values (in general). Moreover, this enhancement of free energies in vacuum reaches between 5 to 10 Å outside the surfaces, implying that the cause is of a long-range nature.
For all free energy profiles shown in Figure 2(e-h), it could be observed that free energy minimum appeared closer to the surfaces with adsorbed MEG and water, compared to surfaces with pure adsorbed water phases.
An effect that may have contributed to the displacements of the free energy minima towards the surfaces (due to adsorbed MEG) is the lowering of the first density peaks of water, thus  yielding lower barriers for hydronium and bicarbonate to cross to reach the respective surfaces.

Charge densities
In the same manner as the mass density profiles seen in Figure  8 were sampled, we sampled the charge density fluctuations in proximity to calcite and hematite. Charge densities originating from water, with and without adsorbed MEG, are shown in Figure 10 together with charge density originating from adsorbed water and adsorbed MEG combined.
Adsorbed MEG caused a significant attenuation of charge density originating from water near calcite (approximately same attenuation for both positive and negative charge densities). Total charge density, originating from both adsorbed water and adsorbed MEG was attenuated in the same manner, but less. Since the charge density profiles with and without MEG had a similar shape in the vicinity of calcite, Coulomb forces from water on ions traversing the adsorbed water were in approximately the same directions (with and without adsorbed MEG) but were weaker with adsorbed MEG.
As can be seen from Figure 10, at hematite, adsorbed MEG caused a noticeable reduction, as well as qualitative change, of the charge density fluctuations in the region between 1.5 and 4 Å from the surface. However, the adsorption geometry of MEG was such that the total charge density of water remained approximately unchanged after MEG adsorption. As was the case for calcite, the forces originating from Coulomb interactions of water on the traversing ions were weaker in case of adsorbed MEG.

Robustness of the results
The models applied to obtain our results are simplified compared to the physical systems we simulate. For example, realworld surfaces will often contain defects and irregularities that are not accounted for in the simulations. Moreover, charges can differ between atoms at a surface and the corresponding bulk atoms of the surface slab in ways that are difficult to mimic using point charges. Using readily available force-fields we also need to make assumptions about their transferability.
In Figure 11, we have compared the water density profile near calcite from our simulations (see Figure 8) with experimentally obtained values from the literature [53]. Our calculations show five distinct peaks, where the first two peaks and the last peak are no more than 0.3 Å from the three peaks measured from X-ray reflectivity (XR) experiments [53]. Thus, the applied water and calcite models show good agreement with experiments without introducing surface defects or irregularities.
To further investigate the robustness of our models, we performed a series of simulations where selected parameters were altered. Figure 12(a) shows the effect on free energy profiles of MEG near calcite at 298 K and 1 atm, due to changing between a cutoff radius of 9 and 11 Å, as well as changing between a 26×26×51 and 96×96×192 k-space grid for the Ewald summation spline interpolation. Figure 12(b-c) shows the effect on free energy profiles of hydronium ions near calcite and hematite, respectively, due to varying the ABF sampling from using a single ABF molecule with seven windows of 10 ns to using two ABF molecules with 20 windows of approximately 20 ns. In Figure 12(e,f), the same change in ABF sampling was done for bicarbonate ions near calcite and hematite, respectively.
We also made small perturbations to the hematite surface geometry by moving ions from one side of the surface to the other (surface reconstruction) resulting in reduced dipole moment, where the effect on free energy profiles for hydronium and bicarbonate ions near hematite is shown in Figure  12(d,g), respectively. Note that the applied changes are very similar to creating surface defects and irregularities on the surface.
The forcefield parameters between bicarbonate and calcite, as well as between bicarbonate and hematite, were exchanged with ones obtained using quantum chemical calculations on simplified systems (see Supplementary Information). The observed effects of perturbations to the surface force field on the corresponding free energy profiles are shown in Figure  12(h,i), respectively.
It can be seen that nearly all positions of global and local free energy maxima and minima survived changes in cutoff  radius, Ewald summation parameters, force-fields, as well as minor changes in surface geometry (such as introduction of defects and irregularities). Moreover, the same invariance seems to be present for hydronium, bicarbonate, as well as for MEG free energy profiles. We therefore believe that the shape of the free energy profiles, not including the relative values, is determined largely by the overall geometry of the systems (i.e. the shape of the free energy profiles is only affected by larger changes in the surface geometries), as well as the primary features of the applied force fields. Hence, our models can be used to predict the overall shapes of the free energy profiles, which is also the primary goal of the performed simulations.

Summary and conclusion
To investigate the effects of adsorbed MEG on corrosion of the walls of pores formed between pipelines and cement, we studied changes in free energy of hydronium and bicarbonate, due to adsorbed MEG, on calcite and hematite in an aqueous Figure 12. (Colour online) Free energy profiles change in response to variations in force field cutoff radius, Ewald summation k-space grid resolution, variation in simulation time (i.e. sampling time) and force field adjustments (in this case a fit to a similar system where forces are calculated using quantum mechanics). Thick lines were calculated by taking the average (along the y-axes) of the shaded areas with the same colour. Green and blue denote MEG, purple denotes H 3 O + , red denotes HCO − 3 and yellow denotes HCO − 3 compared to a similar system where forces were fit using quantum mechanical calculations. To better visualise the range of free energy values between compared profiles at each position, the space between each compared profile is filled with a colour representing the molecule for which the profile belongs. environment. Origins of the free energy changes were studied through density profiles, radial distribution functions, as well as charge density fluctuations the first nanometer outside the respective surfaces.
Density profiles showed that adsorbed MEG displaced water molecules from the first three water layers to bulk and a fifth layer of water, which effectively reduced water-hydronium and water-bicarbonate interactions close to the surface, which can contribute to the free energy minimum moving closer to the surfaces.
Free energy minima of hydronium and bicarbonate ions, other than the global one, were observed between the free energy minima and the free energy in bulk. These were found to originate from each layer of water, possibly combined with the potential originating from the surfaces. How peaks and valleys of free energy profiles of the ions correlate to peaks and valleys of water density is dependent on geometry, charge distributions and short-range behaviour of the ions. Similar atoms competed for a finite number of adsorption sites, which plays a significant role in determining the effect that adsorbed MEG can have on the desorption free energy.
It should be noted that, even if a given free energy profile is moved closer to the surface or gets a higher adsorption free energy magnitude, a series of contiguous free energy barriers with maximum heights of 5 kJ/mol (resulting in an adsorption free energy barrier) would require a temperature of approximately 400 K to be overcome by an ion before that ion could reach the surface (i.e. using E = 3/2k B T). Hence, it is not certain that ions will adsorb easily, even with large magnitude of free energy of adsorption.
From this work, we see that the changes in free energy profiles of a molecule close to a surface, due to other adsorbed species, is extremely complex and is dependent on numerous factors, such as adsorption geometry, geometry of the molecules, charge distribution, hydrogen bond sites, as well as surface geometry.
We believe that this work may provide valuable information for further development of theoretical models for determining corrosion processes taking place between iron pipees and cement casings in aqueous MEG environments. In future works it could also be useful to investigate the effects of CO 2 within the studied systems, as well as the effect adsorbed MEG would have on CO 2 .