First Steps towards Understanding the Non ‐ Linear Impact of Mg on Calcite Solubility: A Molecular Dynamics Study

: Magnesium (Mg 2+ ) is one of the most common impurities in calcite and is known to have a non ‐ linear impact on the solubility of magnesian calcites. Using molecular dynamics (MD), we observed that Mg 2+ impacts overall surface energies, local free energy profiles, interfacial water density, structure and dynamics and, at higher concentrations, it also causes crystal surface deformation. Low Mg concentrations did not alter the overall crystal structure, but stabilised Ca 2+ locally and tended to increase the etch pit nucleation energy. As a result, Ca ‐ extraction energies over a wide range of 39 kJ/mol were observed. Calcite surfaces with an island were less stable compared to flat surfaces, and the incorporation of Mg 2+ destabilised the island surface further, increasing the surface energy and the calcium extraction energies. In general, Ca 2+ is less stable in islands of high Mg 2+ concentrations. The local variation in free energies depends on the amount and distance to nearest Mg in addition to local disruption of interfacial water and the flexibility of surface carbonate ions to rotate. The result is a complex interplay of these characteristics that cause variability in local dissolution energies. Taken together, these results illustrate molecular scale processes behind the non ‐ linear impact of Mg 2+ concentration on the solubility of magnesium ‐ bearing calcites.


Introduction
Calcium carbonate minerals, and in particular calcite, are extensively found in nature, for example as the main component in marble, limestone and chalk, and it is therefore one of the most abundant non-silicate minerals at the Earth's surface.Furthermore, calcite has proven its functionality and usability in different contexts in both natural and industrial processes [1].To enhance the understanding of the mechanisms behind observations on the formation or dissolution of calcite, theoretical models describe calcite growth and dissolution processes at many different scales.Pore-scale and continuum scale models are becoming more and more sophisticated, yet there are still important discrepancies between model outcomes and experimental observations, e.g., [2][3][4][5].Similarly, statistical (e.g., [6]), kinetic (e.g., [7]) and ion-by-ion (e.g., [8][9][10][11]) models have also revealed discrepancies between model and experimental observations.Possible reasons for those discrepancies are further elucidated using atomistic-scale approaches and ab initio models, for example in [12][13][14].
Impurities are more the rule than the exception in natural calcite and are commonly found not only in different concentrations but also in a range of different types of impurities [15].One of the most common impurities found in calcite crystals is the magnesium ion, Mg 2+ .Magnesium carbonate (magnesite) is isomorphic with calcite and a wide range of solid solutions are therefore formed between these end-members.Generally, two categories are employed to differentiate between low and high Mg 2+  concentrations, with arbitrary limits [16][17][18][19], i.e., low-Mg 2+ calcite in which a maximum of 2 2% of the Ca 2+ is replaced by Mg 2+ , and high-Mg 2+ calcite with 8% to 40% of Mg 2+ substituted for Ca 2+ [20].Beyond 40% we cannot speak of calcite anymore, but should consider the mineral as dolomite or huntite instead.
According to computer simulations, the incorporation of magnesium impurities into the step edges of a calcite surface is initially thermodynamically very favourable and is independent of the growth site's topography (obtuse or acute edges) [21,22].It has been shown experimentally that during calcite growth Mg 2+ incorporates more into acute edges than obtuse steps [23], although in more recent computer simulations it was found that Mg 2+ is randomly distributed over the calcite crystal [24].Surface studies have also shown that high concentrations of impurities in the growth edge will poison the edge and thereby inhibit further crystal growth [21,22,25].Furthermore, it was shown that it is thermodynamically unfavourable to sequentially incorporate CaCO3 units next to a MgCO3 edge corner in a growing calcite crystal [21,22].In contrast, a mechanism that was shown to be thermodynamically feasible is to incorporate magnesium via diffusion of Mg 2+ towards the surface, followed by Ca 2+ -Mg 2+ exchange with an existing surface calcium site [25].Experimental work has shown an increase in initial growth rate when a metallic (Mn 2+ , Zn 2+ , Co 2+ and Cd 2+ ) ion was adsorbed or exchanged at the calcite surface [26].Other studies concluded that calcium units were dissolving faster than Mg 2+ units [27,28], resulting in non-stoichiometric (incongruent) dissolution.Other, more recent, experimental work has concluded that the dissolution rate of calcite is strongly but nonlinearly dependent on the Mg 2+ concentration in the crystal [29].
Arguably one of the most difficult variables to capture with dissolution rate laws is the impact of impurity ions such as Mg 2+ in the solution and in the dissolving calcite crystal.Experimental work using atomic force microscopy (AFM) [30][31][32][33] has shown alterations in crystal geometries as a result of variations in dissolution velocities of calcite due to the incorporated magnesium and its distinct preference for (interactions with) different surface sites.Additionally, the effect of salt ions on water structure and dynamics [34] and the consequence for calcite dissolution, with an increase in water viscosity and a resulting decrease in diffusion coefficients [12,35], has been widely suggested as an explanation for observed variations in dissolution behaviour.
Molecular dynamics (MD) simulations have been used extensively to study reactivity of calcite crystals under different conditions.Among others, these studies have investigated the importance of variability in reactivity of edges and kink sites on a 101 4 surface and defined rate equations in up-scaled models [36][37][38][39][40] or the impact of confinement and solution composition [12,[41][42][43].Moreover, MD was used to study the free energy landscape of adsorption and dissolution of metal and carbonate ions on calcite surfaces [38,39,44].As in experiments [10,11,45] and unbiased MD simulations [22,36,37,46], free energy studies do not show consensus in their observations on the acute and obtuse edges [38,40,47].
The aim of this work is to provide atomistic insight into the non-linear impact of Mg 2+ on the solubility of calcite, and to contribute to the understanding of the mechanism(s) in place and responsible for the variable reactivity with Mg concentration in calcite.In the literature, the focus has been on explaining the reactivity of Mg 2+ and its poisoning effect on the reactivity [21,22,[48][49][50].Using molecular dynamics, simulating Mg 2+ rich systems is complicated due to the slow exchange rates of water around this cation, necessitating timescales beyond what is currently feasible to gather statistically meaningful data.Consequently, we describe the reactivity of the crystal surface in terms of stability of the surface calcium ions instead.We have studied the residence time of water coordinated to surface Ca 2+ , interfacial water structure, the water diffusion, the deformation of the crystal, the surface energies and eventually the free energy profiles of the extraction of surface calcium ions.Our results reveal that Mg 2+ incorporation at a low concentration does lead to a local stabilisation of surface Ca 2+ , whereas a higher Mg concentration disturbs the crystal structure, leading to local variabilities in the crystal surface, which are expressed in local destabilisation of Ca 2+ , thereby facilitating the dissolution of an adjacent CO3 2− and enhancing dissolution [39] compared to pure calcite.At the macroscopic scale, these observations imply a lower solubility than pure calcite when the concentration of Mg in calcite is low (<2%) and a higher solubility with increasing Mg, in agreement with the nonlinear solubility observed experimentally.

Methods
We simulated two different calcite 101 4 crystal surfaces in contact with pure water.One of the crystals had an island of 16 crystal cation-anion units adsorbed onto the flat 101 4 surface and an edge pit of the same size on the opposite side of the crystal slab, whereas the other surface was atomically flat.Since the edges of the small pit are prone to interact [36], we have focused on the edges of an island that can also represent the edges of a bigger pit.Both crystals were simulated with and without 30% magnesium randomly distributed over the entire crystal slab, representing the high Mg 2+ calcite system.In addition, we simulated a flat 101 4 surface with one magnesium impurity in the surface, representing the low Mg 2+ calcite.In total we thus had five different molecular systems.The classical molecular dynamics (MD) simulations were performed in the large-scale atomic/molecular massively parallel simulator (LAMMPS) [51].We used the velocity-Verlet integrator [52,53] and the Nosé-Hoover thermostat [52,53] to integrate the equations of motion with a time step of 1 fs and to maintain the temperature at 300 K.The simulations were carried out in constant number of particles, volume and temperature (NVT) ensemble with a relaxation time for the thermostat of 0.1 ps.The unbiased simulation time was 10 ns with 1 ns of equilibration period.
System details.The calcite slab was 3.3 nm thick, which corresponds to 12 atomic layers of material, and was previously shown to be thick enough to have no interactions between the surfaces through the crystal slab [14].A surface of 4.8 nm by 4.0 nm was exposed to a water column of 6.5 nm ensuring bulk-like water behaviour in the centre of the column.After an initial equilibration period, eight out of the twelve layers in the centre of the crystal were immobilised and kept in the same position, for two reasons: (i) to speed up the calculations; and (ii) to simulate the bulk crystal and prevent the box from deviating from its reference point, so comparison of the different Mg-doped systems was facilitated.Initially, we exchanged a randomly selected single Ca 2+ by a Mg 2+ in the surface of the material, to be used as a reference of low Mg 2+ concentrations (LMg).Based on previous computer simulations, there is no preference position for Mg 2+ in a calcite crystal [24].For the high-Mg calcite, we therefore randomly distributed Mg 2+ over the entire crystal, replacing 30% of the Ca 2+ ions with Mg 2+ (from here on referred as high Mg 2+ or HMg).In the slab with an island on top, the cations in the acute and obtuse corners were kept as Ca 2+ , and Mg 2+ was disseminated on the island around the corners.The flat surface under the island maintained a 30% random distribution of Mg 2+ .An image of our defective surface model can be seen in Figure 1.
Force field.As described previously [34,54], the simple point-charge flexible water (SPC/fw) in combination with a Buckingham potential [44,55] provided the best water dynamics around calcium in solution when compared to ab initio calculations.For comparison with more recent work on calcite surfaces [12,37,38,41,56,57], the force field to describe the inter-and intra-molecular interactions was taken from Raiteri et al. [58], including the SPC/fw water description [59].An 0.9 nm cut-off was used for the van der Waals forces, except for the SPC/fw-tail force field, in which the cut-off was defined at 0.9 nm, but with a tail from 0.6 nm [58].The exact values used in our simulations can be found in Table S1.
Coordination number.The number of molecules directly coordinated to a central ion was reported as the coordination number.The coordination number is determined from the integral of the radial distribution function (RDF).To include the whole first coordination shell, the cut-off distance was set at 0.35 nm as this distance included the full first peak in the RDF (whether or not split in the case of oxygen in the carbonate (Oc) in the high Mg 2+ calcite), unless stated otherwise.Surface energy.An estimation of the surface energy was obtained by simulating bulk pure calcite, LMg and HMg at zero K (a "dry run" using T = 10 K using the energy of the minimised structure).The energy of the system was then subtracted from the energy of a zero K simulation of the surfaces, where each surface system was first relaxed with water on top at 300 K for 1 ns before the dry zero K run was conducted to minimise the energy.The energy in J/m 2 is the estimated surface energy of the water-relaxed-dry surfaces: Us is the energy of the crystal with expressed 101 4 surface after relaxation at 0 K, Ub is the energy of the bulk crystal, A is the surface area of one side of the crystal.The surface energy indicates the stability of the surface, with a low positive value indicating a stable surface, whereas the more positive the value, the less stable is the surface.An estimation of the "wet" surface energy was calculated after hydrating the surface by placing a water molecule approximately every 1 nm 2 .Only the first four layers could relax, as the rest of the slab was fixed according to its energetic minimum of the dry surface.The energy of one water molecule multiplied by the number of water molecules on the surface was subtracted from Us according to Equation (2): where  is the energy of one liquid water molecule in bulk water, Us is the energy of the configuration excluding the water and Ub is the energy of the bulk crystal.
We have also used Equation (1) to calculate the surface energy, employing the final configuration resulting from the dynamic simulation including water, but excluding the water in the estimation of Us.The standard deviation was calculated to demonstrate the variance of the surface energy and thus, indirectly, the mobility of the ions in each system studied.
Free energy profiles.The free energy calculations were performed to extract information about the reactivity of Mg-calcite surfaces and to compare them to the pure calcite surface.We simulated various single calcium ion extractions from different positions in a crystal surface (with or without an island on top and with different amounts of Mg in the vicinity) into the pure liquid water layer 1.2 nm away from the surface.The PLUMED [60] plug-in for LAMMPS was used to perform all free energy calculations.In a single steered MD simulation, we extracted a specific calcium out of the surface and dissolved it.To have a converged simulation for construction of the energy profile, the frames of the steered MD were used to perform biased umbrella sampling at different heights above the initial surface position.In our exploration of the free energy, we selected the distance perpendicular to the surface as the collective variable (CV), based on the assumption that this distance would affect the free energy the most.Due to convergence issues in the directions parallel to the plane, the energy profile was not stable after 120 ns of simulation and we therefore decided to limit the energy exploration only in the direction orthogonal to the surface, by restricting the movement in the plane parallel to the surface with a strong harmonic spring (spring constant of 5.0 eV/Å 2 ).With this restriction in place, the free energy calculations were carried out for 10 ns.The calculation consisted of more than 61 harmonic umbrella potentials with spring constants of 0.5 eV/Å 2 , or 5.0 eV/Å 2 in positions where necessary to obtain a satisfactory overlap in the sampling along the CV of neighbouring windows (see Section 2 in the supplementary materials and Figure S1 for more details).The unbiased free energy profiles were obtained through self-consistent histogram re-weighting, using the weighted histogram analysis method (WHAM) code [61].For the interpretation of the energy profiles, the ΔGextraction was used, which is defined as the Helmholtz free energy of extraction of a calcium ion from bulk water (~1.3 nm) to its position in the crystal structure.The positive energy indicates that the cation is favoured in the surface position.
Diffusion.The mean square displacement (MSD) was calculated with a modified version of the Fortran code written by W. Smith from the UK Daresbury Laboratory [62].Thereafter, the self-diffusion coefficient was derived as the slope of the MSD as a function of time divided by the number of dimensions times two (for both positive and negative direction).The self-diffusion coefficient of interfacial water parallel to the calcite surface was calculated by adjusting this code to only account for water molecules present in the interface (e.g., below 0.35 nm from the calcite surface) at the scanned time.Due to the small size of the interfacial water layer, water molecules leave the specified layers within short times.To improve the statistics, we modified the code to consider a linear extrapolation of the MSD when the water molecule entered the interfacial water layer between the time-origin and the remaining simulation time, and in the input file we set the time interval between MSD origins to one frame.In this way, water molecules that stayed longer in the specified interface area had more weight in the MSD, and short visiting molecules could be included without wrongly interfering with the MSD [54].
Vibrational spectrum.The vibrational spectrum at different distances from the surface was built from the sum of all the vibrational density of states (VDOS) of each atom involved [63].The VDOS can be calculated using the Fourier transformation of the velocity-autocorrelation function (VACF) of the individual atoms in the water molecule.We used the following equation to calculate the VACF: where Ntsteps and Natm are the number of timesteps and the number of atoms, oxygen and hydrogen in this case; νi is the velocity vector of O or H atoms in the i th water molecule.
To produce the vibrational spectrum of interfacial water (<0.35 nm from the surface), the water molecules in the interface were flagged in the starting configuration, which was compared with all configurations (the final 40 ps of the simulation) to make sure that the flagging was still valid (that is, that the water molecules still remained in the interface).

6
Water density.The water density profile perpendicular to the surface was determined to visualise the layering of water at the interface.The density is calculated based on the distance from the surface to the centre of mass of the water.
where dnx is the number of particles within a water layer with volume Vx.To visualise the water z-density profile above a certain surface position, we only considered the water molecules in the column on top of the position within a rectangular box with width of 0.2 nm by 0.2 nm.Number of hydrogen bonds between water molecules.Another feature to characterise the interfacial water structure is the number of hydrogen-bonds that a water molecule, coordinated to the surface, forms with other water molecules.This interaction stabilises the configuration and is fundamental to create an intermolecular network at the solid-water interface, as well as in the bulk solution.The coordination distances were defined based on the first minimum in the RDF and was set at 0.28 and 0.40 nm for the cation -oxygen and anion-oxygen, respectively [64].A hydrogen-bond was defined with a proton-oxygen distance of 0.245 nm and an angle of less than 30 degrees between the oxygen acceptor, oxygen donor and proton donor [64].
Water exchange frequency.To investigate the water dynamics around dissolved calcium ions, the water exchange frequency (Nex) of water molecules in the first hydration shell of the cation was determined, using the "direct method" as described in previous work [54].The average exchange frequency for the whole surface was calculated by dividing the total number of exchanges by the number of Ca ions in the surface.

Water Density, Structure and Dynamics in the Interface
The water z-density profile on top of the pure, low Mg 2+ flat surface (LMgF) and high Mg 2+ flat surface (HMgF) calcite 101 4 surfaces revealed the layering of the interfacial water (Figure 2).The water profile of LMgF was identical to the profile of water on top of pure calcite, but there was a clear difference observed in the water layering above the HMgF surface compared to the pure and LMgF calcite surfaces.For HMgF, the first peak had a lower intensity and showed a shoulder at shorter distance, indicating higher densities closer to the surface.
The hydrogen-bond network varied with the level of Mg incorporation, while the total number of hydrogen-bonds with water coordinated to Ca was comparable between the different interfaces (Table 1).Our analysis of the interactions between oxygens in water and carbonate molecules showed that water molecules coordinated to Mg 2+ form on average 1.37 hydrogen-bonds with Oc, compared to 1.02 for water molecules coordinated to Ca 2+ .The latter value is irrespective of whether the surface Ca 2+ is located in a pure or HMgF calcite.Almost half of the water molecules coordinated to surface Mg (40.8%) showed more than one H-bond to at least one of the neighbouring surface carbonate oxygen (Oc) (Table 1), in contrast to water molecules coordinated to Ca 2+ , where this was the case for only 15-16%.The lifetime of the formed H-bond on top of Ca 2+ on the HMgF calcite was 6.5% longer than on the pure calcite.On top of Mg 2+ , this was 19.4% shorter compared to a water molecule on top of pure calcite.Furthermore, the water dipole angle with the surface is slightly smaller (0.24°) for water in the interface with HMgF (Section S3 of the supplementary materials for further details).The vibrational spectrum (VDOS) showed a difference in the behaviour of the water molecules present in the HMgF interface compared to the water molecules present in the pure calcite interface.The lower intensity of the peaks for intermolecular stretching and intermolecular O-O-O bonding motion (Figures S2 and S3) indicates that there are fewer hydrogen-bonds that vibrate with the same frequency as seen in the pure calcite interface.The water dynamics were also altered when a high amount of Mg was incorporated.The exchange frequencies for water molecules coordinated to Ca 2+ showed a broader range in HMgF than in LMgF and pure calcite (Figures S4 and S5), and the average frequency was slightly higher (31.4 versus 25.7 ns −1 for HMgF versus pure calcite, Table S2), whereas LMgF showed similar dynamics to the pure surface.The translational diffusion coefficient parallel to the interface (Dx/y) was also comparable for LMgF and pure calcite, in contrast to HMgF, which showed an increase of 38.9%.Further details can be found in Section S5 until Section S8 in the supplementary materials.

Surface Energies and Structural Relaxation
The energies of the surfaces are shown in Table 2 (see alternative way of measuring this in Table S3 and literature values in Table S4).LMgF calcite crystals had a slightly higher surface energy compared to the 101 4 surface of pure calcite for both dry and wet surfaces.In contrast HMgF calcite had a slightly lower surface energy for the dry surface and a decrease of 0.1 J/m 2 on the wet surface was observed compared to the pure wet calcite surface.Moreover, the surface energy of a surface with an island (PureI and HMgI), and edge pit in the opposite side of the slab, was higher for the dry island.This trend was in contrast with the wet HMgI surface which showed the lowest surface energy (Table 2).
To summarise, the overall dry 101 4 surface energy varied according to HMgF < Pure < LMgF < PureI < HMgI and the overall wet surface energy was ordered as HMgF < HMgI < Pure < LMgF < PureI.Evidence of the degree of relaxation in the surfaces can be seen in the distances between the constituent ions, depicted in Figure 3 (and Figure S6).The average cationcarbon distances (Cc) in the surfaces, observed during the simulations, showed variations depending on the presence of Mg (Figure 3a), in contrast to the cation-Oc distances (Figure S6).In Figure 3a, the second peak (~0.41-0.43nm) represents the distance between the cation and carbon of the first carbonate along the c-axis of calcite, as illustrated in Figure 3b.This distance was shortened upon Mg incorporation.The RDF of LMg had similar Ca-Cc distances compared to pure calcite and the Mg-Cc distances were the same as in HMg.S5 and Figure 3a); calcium ions are in green, magnesium ions in purple, carbonate oxygen and carbon in red and grey, respectively.

Free Energy Profiles
The free energy difference between a fully solvated Ca 2+ and a Ca 2+ in the crystal position, ΔGextraction, was calculated for eight different environments on a flat 101 4 surface, in pure, LMgF and HMgF calcite.The values were in the range ~191-230 kJ/mol, where the ΔGextraction for different crystal positions with respect to Mg 2+ ranged from ~197-223 kJ/mol and the ΔGextraction of a Ca 2+ from the 101 4 flat calcite crystal surface was 204 kJ/mol (Table 3).The profiles of the different flat surfaces showed a steep increase in free energy until 0.5 nm and converged at around 0.9 nm, corresponding to the second interfacial water layer and the start of bulk water behaviour, respectively.In Figure 4, we present free energy profiles obtained from the extraction of Ca 2+ from a pure and LMgF calcite 101 4 surface.All calcium ions in the nearest surroundings of Mg 2+ , i.e., without another cation in-between, had a higher ΔGextraction compared to the ΔGextraction of Ca 2+ in pure calcite.For HMgF, the ΔGextraction for surface calcium varies with the local environment of the crystal between 191 to 230 kJ/mol.Higher energy barriers were found when a Mg 2+ ion was in the direct vicinity of the target Ca 2+ .Two out of three Ca 2+ gave a ΔGextraction higher than the energy needed to extract one Ca 2+ from a pure system (Figure 5).To describe in more detail the conditions leading to a difference in ΔGextraction, we used the labels in Table 3 to distinguish between the different Ca 2+ environments, both with respect to the number of Mg neighbours and the surface topography.

Low Mg 2+ Calcite
The ΔGextraction profiles for the LMgF system showed narrower and shallower local minima (Figure 4b) than in the pure calcite (green line, Figure 4b).Furthermore, while the average distances in LMgF were the same as in pure calcite, locally there were small but visible variations (Table S6). Figure 4b illustrates the small local variations in the free energy profiles, water densities and distances for four Ca 2+ ions located at a unique distance from the surface Mg 2+ .Note that the water densities (right axis in Figure 4b) represent the local water z-density of the unbiased simulation, i.e., the water density before calcium extraction.CaLMgF1 had the shortest distance to the Mg 2+ and was the only Ca 2+ for which the two neighbouring Ca 2+ (at ~0.4-0.5 nm) were not at the same distance (the closest at 0.4525 nm and second nearest at 0.4925 nm).The ΔGextraction profile for this surface calcium shows a small dent at the distance where the highest water density is observed (Figure 4b), although there is no strong correlation between the water density and free energy profiles.For CaLMgF2 and CaLMgF3, the free energy profiles, water density profiles and distances to neighbours were comparable (Figure 4b).The three Ca 2+ nearest to Mg 2+ had ΔGextraction which was higher than pure calcite.For CaLMgF4, the second Ca 2+ distance was slightly longer than the other distances and it was the only site with lower ΔGextraction than pure calcite. in small green or colour-coded (cf.Table 3) larger spheres to link them to the free energy profiles in turquoise (▼), pink (•), orange (no marker) and purple ().(b) Overlay of free energy profile (continuous line, left y-axis) of Ca 2+ extracted from its equilibrium position in the crystal (0.0 nm) with the water z-density (dashed line, right y-axis) along the same axis.The colours are referring to the corresponding Ca 2+ equilibrium position in the left image and Ca 2+ in pure calcite in green (+).

High Mg 2+ Calcite
For the free energy profiles determined in the HMg system, three calcium sites were selected that could be differentiated by the number of Mg 2+ in the vicinity.The first site (CaHMgF0-0; Figure 5a) showed average coordination distances to Oc (0.232 nm) close to the average coordination distance in pure calcite.The energy profile showed a shallower energy minimum at its relaxed position in the surface, resulting in a total ΔGextraction that was 13 kJ/mol lower compared to the extraction of Ca 2+ from a pure system.The first and second peak in the water density profile correlated to moderation of the slope of the energy profile, but no other local minima were observed upon extraction.The second Ca 2+ (CaHMgF1-2) had one Mg 2+ next to it and another at a slightly larger distance (Figure 5a).Furthermore, two Mg 2+ ions were present nearby in the second crystal layer.The average coordination distance to Oc (0.282 nm) was somewhat shorter than for pure calcite (Table S5).In the free energy profile of the second Ca 2+ (CaHMgF1-2; Figure 5a), shallow local minima on both sides of the second peak of higher water density (at 0.45 nm) were observed.With a total ΔGextraction that was 26 kJ/mol higher than the pure calcite, this Ca showed the 11 highest energy needed for extraction.Also, the water column on top showed the highest water density at the first peak in the density profile of all calcium sites.This high water density was comparable to CaLMgF1, which also had a Mg 2+ at 0.39 nm (cf.Figures 4a and  5a).Moreover, in both cases the second peak in the water density profile was slightly shifted away from the flat surface.The last investigated site on HMgF (CaHMgF3-5; Figure 5) had the most Mg 2+ ions in its vicinity and a total of five Mg 2+ in the surface layer.It also had one Mg 2+ in its vicinity in the second layer.It is worth noting that the distance of the closest Mg 2+ and CaHMgF3-5 of 0.385 nm was similar to the distance between Ca 2+ and Mg 2+ in dolomite (0.387 nm).The average distance of the other two Mg 2+ and one Mg 2+ in the second layer was 0.486 nm.The average coordination distance to Oc (0.234 nm) is comparable with pure calcite.The Ca site was the only site that showed a clear shallow local minimum at <0.2 nm.The total ΔGextraction had an intermediate value of 5 kJ/mol, which was the closest to pure calcite.This Ca-site showed the lowest water density in the first water layer.

Calcite with an Island
In the Mg-doped island topography on the 101 4 surface (HMgI, Figure 6a), we focused on four different sites of the 4×8 atoms island.Note that initially the Mg 2+ impurities were randomly distributed, causing a 50% impurity level in the island.Subsequently, we moved a few Mg 2+ closer to the corners to make them have equal amounts of impurities before starting the simulations.The Mg-rich island showed relaxation of the surface structure (Table S7), where the local symmetry of the island changed and the distances and positions were altered relative to the pure island.Consequently, the edges in HMgI did not have the same symmetry in the crystal as observed in the PureI edges.The RDFs between island-surface Ca 2+ and Oc of high HMg (Figure 6) showed more distinguishable peaks compared to island-surface Ca 2+ in pure calcite (in Figure 7).
In general, the energy profiles of PureI and HMgI (Figure 6) showed similar shapes with local minima and maxima at the same positions, yet the minima in the HMgI profiles were deeper, for example creating a more substantial minimum around the second peak of water density (~0.5 nm) on the obtuse edge.The total ΔGextraction (Table 3) for the obtuse edge in PureI is 16 kJ/mol higher than the same edge in HMgI.Furthermore, the extraction of Ca 2+ is less favourable from the obtuse edge compared to the acute edge.For the acute edge, the PureI had a total ΔGextraction that was 16 kJ/mol lower than HMgI and was therefore the most easily extracted (dissolved).For the corner calcium sites in HMgI, the free energy profile (Figure 6b) showed several shallow minima and clear local minima on both sides of the water density peaks.Exceptionally, for the calcium site at the obtuse corner of HMgI, these minima were lower in free energy than the minimum for Ca 2+ in the crystal structure.Although the difference in extraction energy (Table 3) between the position in the crystal and fully solvated was only 7 kJ/mol, the total ΔGextraction for this calcium was 24 kJ/mol, due to the low second minimum, yielding a free energy of extraction 15 kJ/mol lower than its pure calcite equivalent.The difference between total ΔGextraction of HMgI and PureI at the acute corner site is 6 kJ/mol and thereby the smallest difference observed on the island.
The HMgI edges showed a similar density for the obtuse and the acute edge, both lower than the PureI edges.The main differences in the water density profiles were observed in the first peak, i.e., the first layer of interfacial water.The water density on top of the obtuse corner in PureI had the highest water density of the island corners.The lowest density was observed on the acute corner of PureI, which also showed a small shift of the peak towards the surface.As in the trend for the corners, the obtuse calcium edge sites CaPOE showed the highest water density in the first peak.

Impact of Surface Carbonate on Energy Profiles
When determining the free energy profiles for calcium extraction from the island surfaces, we observed relaxation of carbonate ions in the vicinity of the extracted calcium sites.To study the impact on ΔGextraction, we froze all carbonate ions in the island and compared the free energy profiles (i.e., the carbonate groups were fixed in position and could not rotate or vibrate, Figure 8).The second energy value for the HMgI surface in Table 3 corresponds to the ΔGextraction for frozen CO3 2− .Due to the frozen carbonate groups, the energy profiles changed significantly, the first minimum was deeper and the profiles are more similar in shape to the profiles obtained from the flat surface.In addition, the shapes and locations of the local minima showed little or no correlation to the unfrozen profiles.When the CO3 2− groups were frozen in their position, the ΔGextraction for the extraction from the corners were calculated at 90 and 97 kJ/mol for the obtuse and the acute corners respectively.These energies are 66 and 60 kJ/mol higher, respectively, than in the system where the carbonate ions were not frozen.To remove calcium from the edges, we calculated the total ΔGextraction to be 175 kJ/mol for the obtuse and 206 kJ/mol for the acute edges, respectively.As such, the free energy required to remove calcium from the acute edge was within the range of free energies obtained for the flat HMg surface.

Discussion
Before we focus on the impact of Mg, we compared our observations for the pure system with those previously published.Note that an extensive study of the force field used here, and the resulting 101 4 surface interaction with water, was compared to Xray data and our results compare well with these previously published values [14,43].The surface energy we obtained for the atomically flat dry 101 4 surface of a pure calcite compares well with previously published values (0.23-0.86 J/m 2 , Table S4) [65], including the value reported for the same forcefield (0.71 J/m 2 ) [66].The difference with this last value is possibly due to the extra relaxation by the MD simulation with water prior to the dry optimization of the slab in our methodology.Experimentally dry surfaces are never completely free of water, and hence lower surface energies were measured.The surface energies for wet surfaces are more comparable to the experimental values for the dry surfaces (0.23-0.34 J/m 2 ); this discrepancy can be explained by the fact that experimental surfaces have defects (covered with water) that are not reproduced in the atomistic models [67].

Influence of Magnesium on Calcite Solubility and Dissolution
It is important to note that we do not see a systematic trend in a single parameter with changing Mg concentration in the solid, which suggests that the impact of increasing Mg content is a complex interplay between the different system characteristics.Consequently, upscaling to, for example, the meso-scale leads to numerous issues (see also [68,69]).This complexity is inherent in the microscopic scale, since there are many parameters that can play a role (simultaneously).The non-linear correlation of solubility with Mg concentration may suggest that the controlling parameters might change with increasing Mg.
The overall crystal and surface structures, as well as the interfacial water, for LMgF was very similar to pure calcite, suggesting that the local environment controls the stability of specific Ca-sites in the surface.For HMgF, the explanation becomes more complex, since there are more parameters in the interfacial water and the crystal structure that are distinct from pure calcite.In the sections below we discuss these parameters consecutively, showing trends and possible explanations for why some calcium sites are more prone to dissolve than others.

Interfacial Water Structure
The interfacial water is key in the growth and dissolution of any material (e.g., [70,71]).The structure and dynamics of the water in direct contact with the material can be described by the hydrogen-bond structure [72].In this section, we focus on the flat surfaces to discuss the impact of Mg on the interfacial water, to avoid the added interplay of roughness.A clear layering of the water until 0.8 nm from the surface was observed for our pure calcite system, which has been seen before on calcite and other minerals in both theoretical [73,74] and experimental studies [75].While the interfacial water structure in LMgF was indistinguishable from that in pure system, the water structure around the Mg 2+ in the HMgF surface showed clear differences.As a result of the higher magnesium content, HMgF attracts water closer to the calcite surface, thereby disturbing the water layering (Figure 2).This behaviour agrees with that shown in previous computer simulations, where Mg 2+ was adsorbed onto the surface [44,55,76].With this disruption in the water layer, the hydrophilic and hydrophobic zones in the layering are both shifted 0.01 nm towards the surface.
The hydrogen-bond network for interfacial water is highly affected by the interaction with the crystal surface, with all vibrational modes shifting to higher energies (Figure S2).The higher intensities in our VDOS spectra indicate a stronger H-bond network directly over the surfaces compared to bulk liquid water for pure, LMg and HMg systems.For HMgF, the somewhat lower vibrational intensities around 110 and 400 cm −1 suggest that Mg 2+ disrupts the H-bond network.This agrees with the lower number in H-bonds on top of Mg 2+ (Table 1) and the altered dipole angle, indicating that the plane of water molecules is positioned more parallel to the surface (Section S3 in the supplementary materials and Table 1).Combined with the shorter Mg 2+ to water oxygen distances (Figure 2, see c), this results in a lower tendency to form H-bonds with other water molecules from, for example, the second water layer.Instead, water molecules coordinated to surface Mg 2+ are more likely to have proton-oxygen interactions with neighbouring carbonates.As a result, water molecules are librating towards neighbouring water molecules and Oc that are at the H-bond distance.Furthermore, the increase in librations (Figure S3, ~600 cm −1 ) for HMgF interfacial water suggests that there are more possibilities to form H-bonds, in agreement with the higher average water density (Figure 2) and shorter H-bond lifetimes.Due to the disruptions in the water layer, surface Ca 2+ sites experience a weaker interaction with water molecules, resulting in higher water exchange frequencies and higher water diffusion coefficients to and from the surface calcium sites (Table S2).
The differences in intensities in the water z-density profiles do not correlate with the energy profiles or the ΔGextraction.The effect of water on the stability of the Ca 2+ in the crystal is therefore minor.This observation is in agreement with experiments [48], although the disruption of the interfacial water structure in HMgF, resulting in higher water exchange frequencies and diffusion rates near surface calcium sites, are likely to contribute towards the release of Ca 2+ and thereby enhance dissolution of the crystal.

Surface Energies and Structural Relaxation
Since higher Mg-bearing calcites are generally observed to be more soluble than pure calcite, the surface energy for HMgF was expected to be higher than for pure calcite [77].Counter-intuitively, the lower surface energy for HMgF obtained indicates a more stable flat 101 4 surface compared to pure calcite.This is most likely due to the strong relaxation observed in the crystal surface, in combination with the re-positioning of water molecules in the interface.Alternatively, in the atomically flat low-Mg calcite, no alteration of the global crystal structure and interfacial water was observed, and a small local deformation at the surface resulted in a 0.01 J/m 2 increase of the surface energy relative to the pure system.Again, this was counter-intuitive, given that low-Mg calcites generally show lower solubility than pure calcite.Despite the differing methodologies, the higher surface energy for LMgF compared well to the trend seen in the literature, where a calcite slab with one Mg 2+ had a slightly higher energy compared to a pure calcite 101 4 surface [78].
In the case of a 101 4 surface with a small island, which may be considered more comparable to realistic surface topographies, we obtained dry surface energies indicating that HMgI calcite is less stable (Table 2).Due to the random Mg distribution, the island contains higher concentrations of Mg 2+ than the average HMgI slab, causing local alterations in the crystal structure (Table S7) that cannot be compensated fully, causing higher stress at the interface.
As can be seen from the surface energies of the wet surfaces (Table 2), water is stabilizing the interface as it interacts with the surface and in particular with surface Mg 2+ .This agrees with previous work where dolomite and magnesite have a lower wet surface energy than pure calcite and the water showed a greater impact on the surface energies for the Mg-bearing minerals [25].Due to the stabilizing effect of water on the surface energy, there is no indication that the HMgF surface is more prone to dissolution than the pure and LMgF.We therefore explored the variability of the surface energy during the simulation relative to the optimised bulk configuration (cf.[79,80]).Note, that the average dynamic surface energies and standard deviations obtained (Table S3) are significantly higher than the surface energies in Table 2, due to the absence of energy minimization.Yet, the order of the values is the same as for the dry surface energies, implying again that the HMgI is the least stable system.The larger variance in surface energy for the surfaces with an island (Table S3), indicates a more mobile surface, with the largest variance for the HMgI indicating a very dynamic surface.As expected, the structure further into the slab and the interfacial water structure and dynamics were affected, which may imply that dissolving a surface island on a HMgI calcite is easier compared to a pure calcite and may be related to the deformed crystal template and the difficulty to find an island configuration that matches the crystal layer below.
The deformation of the crystal due to Mg 2+ incorporation has beenn mentioned before in experiments, and is generally thought to be the reason for the increase in solubility of high-Mg calcite [81].Furthermore, anhydrous experiments found that discrepancies in reactivity can be explained by stress in the crystal caused by the presence of Mg 2+ , independent of the water-material interaction [48].Ab initio simulations supported both experiments and found a carbonate ion tilt by 6% due to the introduction of Mg 2+ into the crystal [82].We also observed some minor tilting of the carbonate ions in the surface, although we did not quantify this behaviour (Figure 3b), while the compression observed for HMgF calcite along the c-axis (Figure 3a) is related to local alterations of Ca 2+ -C distances (Figure 3b).
To summarise, the atomically flat surfaces show trends in overall surface energies that are opposite to those expected from known solubility trends, while the variability in surface energy with time does follow the solubility trend, suggesting that a more dynamic surface reflects a more soluble surface.Increasing structural deformation with increasing Mg content confirms previously reported structural destabilization of the crystal structure.

Local Variability in Surface Free Energies
The trend in the surface energies of the flat surfaces, showing HMgF being more stable due to the structural rearrangement, appears to contradict experimental observations.However, when also considering more realistic surface structures such as edges and kinks/corners, the surface energies, crystal deformation and interfacial water distortion suggest that once the surface contains imperfections, the system may become thermodynamically less stable, increasing the solubility product for HMgF as measured in experiments [83][84][85][86][87].In order to investigate the impact of magnesium impurities, it is worth investigating the influence on the actual mechanism of dissolution, i.e., extraction of calcium ions from surface edge and kink sites as well as from atomically flat surfaces.
In general, calcite crystals will dissolve faster from terrace corners, edges and kink sites (e.g., [30,88,89]) than from a flat surface.The extraction of a cation from an edge or corner generates new kinks (strongly under-coordinated surface sites).The free energy profile for such calcium edge or kink site extraction was investigated in HMgI and pure calcite with an island (PureI).With the notable exception of the acute edge calcium site, the overall ΔGextraction was consistently lower for HMg island surface sites than for the same sites in the pure material (Table 3 and Figure 6b).Generally, it may therefore be concluded that rough HMg surfaces will dissolve more easily via calcium detachment than pure rough surfaces.The first minimum in the energy profiles corresponds to the equilibrium position in the crystal surface, the second minimum is at the same height as the first water layer near the surface (Figure 6b), which shows that the interfacial water molecules slightly stabilise extracted calcium ions.Note that there are subtle differences in interfacial water structure, density and local energy minima (Figure 6; and potentially carbonate molecule flexibility, see below) and that this interplay may contribute to defining the overall ΔGextraction for surface calcium ions.For extraction of PureI and HMgI surface corner and edge calcium sites, the local energy minima within <0.5 nm from the surface suggest that there is stabilisation of the extracted calcium, most likely as an inner sphere complex.For the obtuse corner on HMgI, this inner sphere complex is energetically more stable than its position in the crystal.
The dissolution of flat surfaces first needs new etch pit nucleation events, which is an important (rate limiting) dissolution mechanism at near equilibrium conditions [90,91].The extraction of a single ion from a pure flat calcite face represents an unassisted etch pit nucleation event (i.e., without related sub-surface structural defect), which needs 204 kJ/mol according to our simulations (Table 3).The extracted Ca 2+ from the flat surfaces shows a similar energy profile as published before, although the shallow energy 18 minimum observed around 0.3 nm is more pronounced in our simulation and the total ΔGextraction at ~250 kJ/mol is higher [39].A possible explanation can be the use of a different water force field, whereas our water force field also gave a slightly lower density in the coordination shell of Ca 2+ [44,54].In both our results and the literature, the free energies converge at >0.5 nm away from the surface.The extracted ion can still interact with the vacant site, even though there are layers of water in between.The lack of a local minimum closer than 0.5 nm to the surface suggests that there is no stabilisation of an inner sphere complex at the flat surface.The formation of outer sphere complexes (i.e., the energy minima at >0.5 nm from the surface) initiates the convergence of the free energy profiles.This observation is also seen in the impure flat surfaces, where the shallow minima/maxima are positioned between the first and second water layer.Structural imperfections, such as (local) deformation due to Mg impurities, can assist or inhibit nucleation of new etch pits.Accordingly, in LMgF it is expected that calcium extraction generally needs more energy than from pure calcite.Depending on the local environment of the extracted calcium, we observed higher ΔGextraction for LMgF than the pure mineral, albeit not for all calcium extractions (Table 3).
The more favourable ΔGextraction for one of the LMgF sites may be linked to the flexibility of the coordinated carbonate molecules, since we observed that the substitution of Ca by smaller Mg ions generates more space in the crystal structure and therefore more flexibility in the carbonate ions coordinated to that substitutional ion.This flexibility may enhance calcium extraction, and we observed a major increase in the ΔGextraction by 300% when interfacial CO3 2− ions were frozen (Table 3).Such flexibility of surface carbonate groups has been found previously in X-ray reflective and diffractive experiments, where the flexibility is mainly expressed in the standard deviation of the tilting angle of carbonate, which was reportedly 13 times larger for the surface layer compared to the bulk [92,93].Although this variation could be due to the two main orientation directions of CO3 2− observed in AFM [45] and MD (e.g., Figure 1) or the surface roughness expected in every non-ideal, cleaved surface, there was no evidence found in the X-ray data for the existence of two distinguishable groups [92].Furthermore, based on the 'roughness factor' introduced by Fenter et al., the impact of the surface roughness is negligible for X-ray reflectivity data interpretation [92].It is very likely that the observed CO3 2− flexibility in X-Ray Diffraction (XRD) is due to imperfections and thereby indicates the presence of rotating carbonate sites on the calcite surface at geometrically more open surface sites (i.e., edge and kink sites).This is supported by the small flexibility observed in MD calculations of a perfectly flat calcite surface here and in previous investigations [14], indicating that the increased coordination of corner and edge sites with water molecules might not be the only reason for more favourable dissolution energetics and higher reactivities of such sites, but that more opportunity for carbonate molecule rotation and flexibility plays a (potential key) role as well.When Mg 2+ is in the direct vicinity of Ca 2+ , the shared coordinated carbonate molecules are held in place more strongly by the stronger bond with Mg 2+ .Consequently, carbonate is less flexible and less likely to assist in the extraction of calcium.When distortion occurs in the local crystal environment, for example due to the longer Ca 2+ -distance (Table S6), the carbonate molecules are more flexible, resulting in a ΔGextraction that is slightly below the energy value of Ca 2+ in pure calcite.
The lowest value for assisted etch pit nucleation (191 kJ/mol) was observed for the flat HMg surface, although not all calcium extractions led to lower ΔGextraction than in pure calcite (Table 3).These results indicate that it may be energetically more favourable to nucleate new etch pits on HMg calcite than pure calcite surfaces, although this strongly depends on the local magnesium distribution.Two out of three ΔGextraction derived from the HMgF agree with the calculated surface energy, indicating a more stable surface.In one of these sites, smaller distances between Ca and Oc are observed, which causes the Ca 2+ to be better charge-compensated and increases ΔGextraction by 26 kJ/mol.In the other site, the local environment is structurally more similar to dolomite than high-Mg calcite and shows an increase in ΔGextraction by 5 kJ/mol, despite the larger Ca-Oc distances and higher water-exchange frequencies (Table S2).The calcium site with the most favourable ΔGextraction showed unaltered bond distances compared to pure calcite, so it appears that here there is no extra compensation from the crystal, unless the nearby presence of Mg ions alters the local charge distribution.Note that classical MD is not capable of revealing such differences at the electronic level and previous ab initio calculations have not reported alterations in charge distributions upon Mg incorporation, but is found to stiffen the calcite structure [82,94] Alternatively, it may be that this calcium site is less stable than Ca 2+ in pure calcite, due to the characteristics of interfacial water: with a higher diffusion (Table S2) and less structured water (Figure 2 and Figure S3), Ca 2+ is more easily extracted and solvated.
To conclude, calcium extraction is generally more favourable from edge, kink and flat surface sites in a high-magnesium surface than from pure calcite, depending on the subtle interplay between crystal surface and interfacial water structure and density that affects local energy minima.Moreover, the flexibility and opportunity for rotation of surface carbonate molecules may contribute to this complex interplay.This behaviour is particularly clear on atomically flat surfaces, where magnesium impurities limit this flexibility, in particular in low-Mg calcite, thereby potentially playing a key role in the non-linear impact of Mg on calcite solubility.

Implications for the Influence of Magnesium on Calcite Growth
The forced MD may also be considered to represent the molecular scale energetics of ion approach and attachment/adsorption, when regarding the reversed process of ion extraction.For example, in contrast to the enhanced dissolution of a high-magnesium island surface, it will be more difficult to grow such a feature on a HMgF surface compared to growth on pure calcite.In particular, the deformed crystal template makes it more complicated to match the structure and propagate a step edge.AFM experiments have revealed similar behaviour; the first monolayer of Mg 2+ calcite grew normally on a pure calcite template crystal, but subsequent monolayers grew significantly more slowly [31].In addition, AFM experiments have shown that Mg 2+ is found to be a kink poisoner, by occupying a growth site, which drastically slows down growth [50,[95][96][97].In case of growth, our results strongly suggest Ca 2+ is less stable on the high-magnesium surface (island) and less prone to approach and attach compared to Ca 2+ on a pure calcite crystal.Moreover, we observed that the formation of adsorbed inner sphere complexes of calcium likely precedes incorporation of that calcium in acute corners and edge sites of PureI and HMgI, as well as the obtuse edges of the pure calcite island.In contrast, on the HMgI obtuse edge, the formation of such an inner sphere complex may be energetically more favourable than incorporation into the edge, potentially inhibiting the growth of the obtuse edge on high-Mg calcites.

Conclusion
We have performed molecular dynamics simulations of 101 4 calcite surfaces with and without different Mg 2+ concentrations and with and without an island on the crystal surface.With these simulations, the impact of Mg 2+ was determined on the overall and local solid-water interface structure, energetics, as well as the free energy profiles for the onset of calcite dissolution via Ca 2+ removal.
The overall 101 4 surface energy varies according to HMgF < Pure < LMgF < PureI < HMgI.Strong relaxation of the crystal structure was observed in HMgF and HMgI.LMg calcite relaxed its structure to a lesser extent than HMg calcite, showing a crystal structure very similar to pure calcite.
The average free energies for calcium ion extraction from the pure and Mg-doped surfaces followed roughly the opposite trend as the overall surface energy differences.The average energy needed to remove calcium ions increased from HMgI PureI < Pure < LMgF = HMgF, although very large local differences (from −13 kJ/mol to + 26 kJ/mol) were observed.The large variations and the resulting appearance of lower ΔGextraction, 20 compared to pure calcite, lowers the threshold for unassisted nucleation of new etch pits locally, in particular for HMgF.
The local variation in free energy (minima) depends on the amount and distance to the nearest Mg, in addition to local disruption of interfacial water and the flexibility for carbonate ions to rotate.Local configurations were observed to be less stable when Mg 2+ was nearby in the surface, supporting experimental data in which calcite with higher percentages of Mg 2+ has a higher solubility compared to pure calcite up until the dolomite ratio (i.e, Ca:Mg = 1).Some of the free energy profiles showed a local energy minimum where the highest interfacial water density was observed (reflecting the hydrophilic nature of calcium).
Based on the interfacial water structure and dynamics, the surface energies and the ΔGextraction, the low-Mg calcite surface is comparable to pure calcite, although locally Mg 2+ induces stabilization of neighbouring Ca 2+ , which results in slightly unfavourable new etch pit nucleation energetics.We can conclude that low concentrations of Mg 2+ tend to stabilise Ca 2+ and increase etch pit nucleation energies, and thereby decrease the calcite solubility, whereas higher concentrations of Mg 2+ lead to deformation of the surface crystal and interfacial water structure, which leads to local variabilities.The result is a thermodynamically less stable crystal and a higher solubility of high-Mg calcite compared to pure calcite, including more facile unassisted etch pit nucleation.Taken together, these results illustrate the molecular scale processes and demonstrate the first steps towards understanding of the non-linear impact of Mg 2+ on the solubility of magnesium-bearing calcites.

Supplementary Materials:
The following are available online at www.mdpi.com/2075-163X/11/4/407/s1, Figure S1: An example of the histograms obtained from all umbrellas along the CV, Figure S2: VDOS of interfacial water on top of a pure calcite (blue) and HMgF calcite (orange).For comparison and forcefield specific peak positions, in green the VDOS spectrum of pure bulk water, Figure S3: Zoom of spectrum range related to characteristics of the hydrogen bond network.Interfacial water on top of a pure calcite (blue) and HMgF calcite (orange).Pure bulk water (green), Figure S4: Variation in exchange frequency of water on top of pure calcite.The colours follow the exchange frequencies indicated in the legend, Figure S5: Variation in exchange frequency of water on top of HMgF calcite.The colours follow the exchange frequencies indicated in the legend.The green big balls represent Mg 2+ as indicated with the red arrow, Figure S6: Radial distribution function for the cation with the oxygen of carbonate in pure, LMg and HMg calcite crystals, Figure S7: Top view of the 4×8 atoms island on top of Pure 1014 surface carbonate groups in grey (C) and red (O) and Ca 2+ in small green or colour and symbol coded larger spheres to link them to the free energy profiles (Figure 6b) in orange (obtuse edge, ▬), turquoise (acute edge, ▼), pink (obtuse corner, •) and purple (acute corner, ), Table S1: Potential parameters, Table S2: Average exchange frequencies and diffusion coefficients of water coordinated to Ca 2+ in a flat calcite surface, Table S3: Surface energies of after relaxation at 300K with water on top.The Us in equation ( 1) was calculated as the energy of the configuration without water, Table S4: A summary of computed relaxed surface energies of calcite, Table S5: Coordination distances in the bulk crystal structure, Table S6: Relevant distances between surface calcium and the neighbouring atoms and calcium coordination numbers at the different Ca 2+ positions near Mg 2+ in the low Mg 2+ 101 4 calcite surface, Table S7: Comparison of relevant distances and coordination in the surface of PureI versus HMgI.Between the selected calcium (Figure S7) and the neighbouring atoms.

Figure 1 .
Figure 1.Snapshot of simulation cells after equilibration.(a) full simulation box with the calcite 101 4 slab in the bottom, carbonate groups in grey (C) and red (O) and Ca 2+ in green.The 6.5 nm water column with red (O) and white (H).(b) the calcite slab with a 16 cation-anion units island; water molecules were left out for visibility.

Figure 2 .
Figure 2. Overall water z-density profile of pure (dashed green, overlaps with dotted red line), low Mg 2+ flat surface (LMgF; dotted red) and high Mg 2+ flat surface (HMgF; yellow) with respect to the flat crystal surface (= 0.0 nm).The lines with the marker () is the integral of the water density.

Figure 3 .
Figure 3. (a) Radial distribution function for the cation with the carbon of carbonate in pure, low Mg 2+ concentration (LMg) and high Mg 2+ concentration (HMg) calcite crystals.(b) Side view snapshots of the calcite crystal showing the deformation

Table 3 . 4 ) 10 Ca 2+ with 1 Figure 4 .
Figure 4. (a) Top view of the LMgF 101 4 surface with Mg 2+ in yellow, carbonate groups in grey (C) and red (O) and Ca 2+ in small green or colour-coded (cf.Table3) larger spheres to link them to the free energy profiles in turquoise (▼), pink (•), orange (no marker) and purple ().(b) Overlay of free energy profile (continuous line, left y-axis) of Ca 2+ extracted from its equilibrium position in the crystal (0.0 nm) with the water z-density (dashed line, right y-axis) along the same axis.The colours are referring to the corresponding Ca 2+ equilibrium position in the left image and Ca 2+ in pure calcite in green (+).

Figure 5 .
Figure 5. (a) Top view of HMgF 101 4 surface with Mg 2+ in yellow, carbonate groups in grey (C) and red (O) and Ca 2+ in small green or colour-coded (cf.Table3) larger spheres to link them to the free energy profiles in purple (), turquoise (▼) and pink (•).(b) Overlay of free energy profiles (continuous line, left y-axis) of calcium ions extracted from their equilibrium position in the crystal (0.0 nm) with the water z-density (dashed line, right y-axis) along the same axis.The colours are referring to the corresponding Ca 2+ equilibrium position in the left image and Ca 2+ in pure calcite in green (+).

Figure 6 .Figure 7 .
Figure 6.(a) Top view of the 4×8 atoms island on top of HMg 101 4 surface with Mg 2+ in yellow, carbonate groups in grey (C) and red (O) and Ca 2+ in small green or colour and symbol coded larger spheres to link them to the free energy profiles in orange (obtuse edge, ▬), turquoise (acute edge, ▼), pink (obtuse corner, •) and purple (acute corner, ).(b) Overlay of free energy profile with the water density along the same axis.The dashed line represents the results for the equivalent calcium site in the pure system, i.e., the Ca 2+ are in a similar position as in (a), but no Mg 2+ is present in the surface.

Figure 8 .
Figure 8. Free energy profiles of Ca 2+ extraction from PureI (dashed lines), HMgI (continuous line) and HMgI in which the island-carbonate molecules are frozen (continuous blue line).

Table 1 .
Interfacial water dipole angle and H-bond properties.

Table 2 .
Overview of the surface energies of the simulated systems