# Equilibrium isotope fractionation in carbonate minerals: Role of Mg-Ca distribution and thermal effects

Paola Patron, Carolina Martinez, Alin M Elena, Marc Blanchard, Carlos Pinilla

► To cite this version:

Paola Patron, Carolina Martinez, Alin M Elena, Marc Blanchard, Carlos Pinilla. Equilibrium isotope fractionation in carbonate minerals: Role of Mg-Ca distribution and thermal effects. Geochimica et Cosmochimica Acta, 2025, 394, pp.53-69. 10.1016/j.gca.2025.01.031. hal-05017963

HAL Id: hal-05017963

https://hal.science/hal-05017963v1

Submitted on 2 Apr 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
temperature conditions. This approach neglects structural multiplicity and probability
24
distributions that may identify other relevant configurations under ambient conditions. In this
25
study, we investigate Mg and Ca isotope fractionation in the Ca₁-xMgxCO₃ solid solution
26
(calcite, dolomite, magnesite) and in aragonite. Our focus is on the impact of Mg
27
substitutions, structural multiplicity, and thermal probability on magnesium (26/24Mg) and
28
calcium (44/40Ca) isotope fractionation factors. Through a comprehensive analysis of all the
29
configurations for several compositions, we observe a significant concentration effect with
30
an increase in Mg reduced partition function ratio β associated with Mg substitution for Ca
31
when using the minimum total energy or maximum probability approaches. However, this
32
dependence is almost negligible when a more dedicated thermal probability weighted
33
average is carried out, thus highlighting the effect of structural multiplicity and thermal
34
probability when estimating equilibrium isotope fractionation factors in solid solutions. We
35
find that configurations in which the Mg atoms are grouped together exhibit reduced lattice
36
distortion, resulting in lower Mg β-factor. Furthermore, our findings highlight a direct
37
correlation between stress and disorder minimization and elevated Mg β-factors in ordered
38
dolomite structures. We demonstrate that a weighted average of fractionation factors on
39
symmetrically non-equivalent structures of Mg bearing calcite at a given composition
40
generally yields values that can differ up to ≈1 ‰. This approach proves to be reliable in
41
capturing the complexity of composition variation when computationally estimating isotope
42
fractionation within solid solution systems. It emphasises the importance of understanding
43
the structural implications of Mg incorporation on isotope fractionation, shedding light on
44
the complexities inherent in carbonate mineral formation. Overall, our research provides
45
valuable insights into the parameters influencing Mg isotope fractionation in carbonate
46
minerals, contributing to a deeper understanding of the geochemical processes governing
their formation.

49 **Keyword**

Mg and Ca isotopes, carbonate minerals, ab-initio calculations, concentration effects,
structural multiplicity, thermal effects.

52 **1. Introduction**

Carbonate minerals such as calcite and aragonite (two polymorphs of CaCO3), whether of
biological or abiotic origin, are likely to record the environmental conditions at the time of
their formation. They can therefore be used to reconstruct paleoclimates and
paleoenvironments, based on their elemental and isotopic compositions. Historically, oxygen
isotopes in carbonate minerals have been used to gain access to past temperatures. More
recently, thanks to analytical advances, the isotopic compositions of cations have opened up
new perspectives (D. Buhl et al., 2007; De La Rocha and DePaolo, 2000; Eisenhauer et al.,
2009; Farkaš et al., 2007; Heuser et al., 2005; Hippler et al., 2009; Mavromatis et al., 2019;
Saenger and Wang, 2014). Here, we are mainly interested in Ca and Mg, which, by
substituting in the calcite structure, form the partial solid solution with the magnesian pole,
i.e. magnesite (MgCO3), where the intermediate composition corresponds to dolomite
(CaMg(CO3)2). To properly interpret isotopic compositions measured in sedimentary
archives, it is necessary to understand all the parameters involved, such as temperature, pH,
speciation of elements in solution, nature and composition of the precipitated mineral phase,
precipitation rate or mineral growth kinetics (Gussone et al., 2005; Schott et al., 2016; Tang
et al., 2008).

With regard to Ca isotopes, previous experimental studies first focused on the role of
temperature on isotopic fractionation between carbonate mineral (calcite or aragonite
precipitated inorganically or from foraminiferal culture) and fluid (Gussone et al., 2003).

Other works have subsequently demonstrated that mineral growth kinetics play a major role in this isotopic fractionation (Alkhtib and Eisenhauer, 2017a, b; Lemarchand et al., 2004; Tang et al., 2008), as does aqueous calcium speciation (Moynier and Fujii, 2017). More recently, the three-isotope method has been employed to better constrain equilibrium Ca isotope fractionation between calcite/aragonite and fluid, resulting in $\Delta^{44/42}\text{Ca}_{\text{solid-fluid}} = -0.02 \pm 0.13\‰$ and $-0.80 \pm 0.10\‰$ at 25°C, for calcite and aragonite respectively (Harrison et al. 2023).

Similar experimental work was carried out to investigate the parameters controlling Mg isotope fractionation. Inorganic calcite precipitation experiments by homogeneous nucleation or via an amorphous calcium carbonate precursor phase enabled to assess the impact of growth kinetics and formation mechanism (Immenhauser et al., 2010; Li et al., 2012; Mavromatis et al., 2013; Mavromatis et al., 2017). A remarkable result is that Mg isotope fractionation decreases with increasing calcite growth rate, while conversely Ca, Sr and Ba isotope fractionation increases (Mavromatis et al., 2013; Schott et al., 2016). The distinct behavior of Mg stems from the high hydration free energy of the $Mg^{2+}$ ion relative to $Ca^{2+}$, $Ba^{2+}$ and $Sr^{2+}$. Although aqueous $Mg^{2+}$ undergoes less hydrolysis than other divalent transition metals, it nevertheless forms relatively strong complexes with inorganic and organic ligands which affects mineral-solution fractionation (Schott et al., 2016). Magnesian calcite formation pathways either by direct precipitation from solution, or via an amorphous calcium carbonate phase, lead to calcite formation with highly variable Mg contents of up to 20 mol% Mg in the latter case (Purgstaller et al., 2016). These variations in the chemical
93

composition of the mineral could impact the crystal structure and consequently equilibrium
isotope fractionation factors.

94

95

For dolomite, the study of natural samples has led to inconsistent Mg isotope fractionation
factors (Fantle and Higgins, 2014; Geske et al., 2015; Higgins and Schrag, 2010; Mavromatis
et al., 2014). To shed light on this, experimental work by Li et al. (2015) showed that the
degree of structural order in dolomite influences the fractionation factor, with $^{26}$Mg
enrichment in ordered dolomite. These dolomite synthesis experiments were carried out
under hydrothermal conditions (i.e. with a temperature range from 130 to 220°C) in order to
overcome the obstacle of abiotic dolomite synthesis at room temperature (Land, 1998).
Beyond the issue of cation order, the question of the effect of chemical composition remains,
since natural dolomites commonly display variable Mg contents (Drits et al., 2005).

96

97

98

100

101

102

103

104

Computational methods complement experimental works by determining isotopic
fractionations associated with well-constrained systems or mechanisms. Some studies have
focused on predicting equilibrium isotope fractionation between aqueous Ca²⁺ or Mg²⁺ and
carbonate minerals, employing various computational approaches (Gao et al., 2018; Pinilla
et al., 2015; Rustad et al., 2010; Schauble, 2011). Other studies, as here, focus on the possible
effect of mineral chemical composition on isotopic fractionation (Wang et al., 2017; Wang
et al., 2019) and the effect of relative cation position, i.e. structural order (Son et al., 2020).

105

106

107

108

109

110

111

All of these studies show that the reduced partition function ratio (i.e. isotopic fractionation
factor $\beta$ of Mg increases in the order calcite < magnesite $\le$ dolomite, with the aqueous Mg$^{2+}$
ion close to or above dolomite. In detail, with somewhat different methods and limited data,
this has led to some discrepancies particularly for the effect of Mg content in calcite on Mg
isotope fractionation between calcite and solution. Indeed, data from Pinilla et al. (2015) and

112

113

114

115
Wang et al. (2019) suggest opposite trends, most likely explained by different structure optimization strategies (Wang et al., 2019; Son et al., 2020). However, all works agree that the effect of mineral chemical composition on equilibrium isotope fractionation is non-negligible.

Most of these computational works have been limited by the system size as a result of expensive DFT calculations for systems beyond a hundred of atoms and are therefore limited to the substitution of one or two Ca atoms by Mg atoms within the simulation cell.

Furthermore, in these works, the discussion has been based on results obtained from minimum-energy configurations that in principle should be the most stable at zero temperature conditions, thus neglecting structural multiplicity as well as probability distributions that can lead to other structures being the most relevant at ambient conditions.

The only exception is the Boltzmann-averaged model of Mg-aragonite in Son et al. (2020).

In the present work, we use a combination of classical force field and ab-initio results to study the Mg and Ca isotopes fractionation in the Ca₁-xMgxCO₃ solid solution (calcite, dolomite, magnesite) and aragonite. Our work has been partially inspired by a previous work by Wang et al. (2017). However, our methodology has allowed us to investigate a wider range of Mg concentrations by increasing the number of Ca atoms substituted by Mg and the simulation cell size, in order to focus not only on the Mg-O bond distance but also on Mg-Mg coordination as parameters to understand structure stability, as well as, compositional effects on the isotopic fractionation factors. Furthermore, our work presents a detailed discussion on the effect of Mg ordering relation to the stability and the Mg isotope fractionation in dolomite.

## 2. Methodology

### 2.1 Equilibrium isotope fractionation factor from harmonic vibrational modes
The equilibrium isotope fractionation factor of an element displaying two isotopic forms Y and Y* between two phases a and b is related to the ratio of the isotope concentration ratios

$$ \alpha(a, b, Y) = \frac{\left(\frac{n_{Y*}}{n_Y}\right)_a}{\left(\frac{n_{Y*}}{n_Y}\right)_b} \quad (1) $$

where $(n_Y)_a$ is the mole fraction of isotopes Y in phase a. The equilibrium fractionation factor between two phases can be related to the reduced partition function ratio $\beta(a, Y)$ of each phase by: $\ln(\alpha(a, b, Y)) = \ln(\beta(a, Y)) - \ln(\beta(b, Y))$. Isotopic reduced partition function ratios are usually presented as $10^3 \ln(\beta(a, Y))$ in permil. The reduced partition function ratio $\beta(a, Y)$ describes the isotopic fractionation properties of a given phase a normalized to the properties of an ideal gas of Y atoms, having the mean isotropic concentration (Méheut et al., 2007):

$$ \beta(a, Y) = \frac{\left(\frac{n_{Y*}}{n_Y}\right)_a}{\left(\frac{n_{Y*}}{n_Y}\right)_{\text{gas}}} \quad (2) $$

Rigorously, the $\beta$-factor must be obtained from the reduced partition function ratio of each isotopologue of the chemical specie under study and its relative abundance (Richet et al., 1977). In general, the $\beta$-factor can be calculated computationally assuming that the mean number of Y* atoms is within the diluted limit, and that the free energy of substitution at a given site does not depend on the same site at neighbouring cells. Under this condition, the $\beta$-factor can be seen as a geometric mean of the reduced partition function ratios of all the singly substituted isotopologues (Méheut et al., 2007). However, it can be shown that if the free energy associated to the substitution does not depend on the isotope nature of the surrounding isotopes and it weakly depends on the site substituted, the $\beta$-factor can be estimated as a single calculation with all sites of the atom of interest replaced at once and leading to similar results (Méheut et al., 2007). It can be shown then that within the harmonic
approximation and using the Teller-Redlich rule for vibrational frequencies (Wilson et al., 1980), the β(a, Y) factor of a solid phase can be written as (Méheut et al., 2007):

163

164

165

where $v_{q,i}$ are the phonon frequencies of the i-th branch with wave vector $q$. $v^*_{q,i}$ correspond to the frequencies of the phase containing the heavy isotope. $N_{at}$ is the number of atoms in the unit-cell and $N$ is the number of sites occupied by the Y atom in the unit-cell. $k_B$ is the Boltzmann constant, and $T$ is the temperature. The second product in equation 3 is usually performed on a grid containing $Nq$ q-vectors to ensure full convergence.

171

## 2.2 Effective bond distance and effective coordination number

A crystal structure is a solid formed by a group of atoms (unit cell) repeated periodically in the three directions of space (Dana and Hurlbut, 1960). An important structural parameter is the coordination number (CN), that is, the total number of closest neighbours of a central atom inside the crystal structure (Allan et al., 2019; Poplavko, 2019).

The effect of structural optimization and the creation of a solid solution structure can affect symmetry operations and atomic distances. When working with solid solutions at different concentrations, the crystallographic network is deformed by elongation of bond lengths, consequently the CN can change, as a result of this, the crystalline system has a certain amount of disorder, than in this case could be more conveniently studied through the use of the effective coordination number (ECN) and effective interatomic distances (d) for the Mg-O (Mg-Mg) defined as (Hoppe, 1979):

182
$$ECN = \sum_i \exp\left[1 - \left(\frac{l_i}{d(Mg-O)^*}\right)^6\right] (4)$$

where d(Mg-O)* (d(Mg-Mg)* ) is the effective interatomic distance, defined as:

$$d(Mg-O)^* = \frac{\sum_i l_i \exp[1 - (l_i/l_{min})^6]}{\sum_i \exp[1 - (l_i/l_{min})^6]} \quad (5)$$

with lᵢ the bond distance between the central atom and its i-th neighbour, and lmin is shortest

Mg-O (Mg-Mg) bond distance found. For this search we have used cut-off distances of 3 Å

for the Mg-O and 7 Å the Mg-Mg nearest neighbours.

## 2.3 Computational Details

The calculations were carried out using the General Utility Lattice Program (GULP) (Gale, 1997; Gale and Rohl, 2003). Ionic interactions were considered in the form of Coulombic attraction (long range), Pauli repulsion, and van del Waals attraction (short range). Full geometry and cell optimization at 1 bar was performed using the Broyden-Fletcher-Goldfarb-Shannon (BFGS) algorithm (Fletcher, 1987). To represent ionic interactions, we have used the classical force field parameterised by Bruneval et al. (2007) that is shown to reproduce correctly the structural and vibrational properties of Ca polymorphs and aqueous systems.

The quasi-harmonic approximation was used to calculate the thermal expansion of calcite, dolomite, and magnesite structures, in order to compare with available data and validate the force field. However, all reduced partition function ratio calculations have been performed using structures optimized at 0 K. In order to create crystalline structures that were independent by symmetry operations we have used the method that D'Arco et al. (2013); Mustapha et al. (2013) proposed initially for the analysis of symmetry in disordered systems.

This method uses a combinatory approach to determine the number of symmetry-independent classes that are subsequently used to calculate structural properties. This study
comprised the estimation of structural properties for a total of 1400 inequivalent structures for Mg concentrations in Mg-bearing calcite in a range from 1 mol % to 25 mol %. Most of the Mg-bearing structures were created based on the rhombohedral calcite supercell with Ca²⁺ ions substituted by Mg²⁺ ions. Information on supercell sizes can be found in figure 1 or the supplementary material with the total number of Ca²⁺ and Mg²⁺ cations provided as guide for each studied system.

The lowest-energy configurations at different Mg concentrations can be seen in Figure 1. Apart from using the minimum-energy configuration to estimate isotope fractionation factors, we have also considered the contribution of each inequivalent configuration weighted by its Boltzmann factor at 298.15 K. The probability $P_m$ of a given structure is thus given by:

$$P_m = \frac{\Omega_m e^{-E_m/RT}}{\sum_m \Omega_m e^{-E_m/RT}} \quad (6)$$

where $\Omega_m$ is the degeneracy or multiplicity of each one of the studied structures for a given concentration and $E_m$ the total energy for this structure obtained using either DFT or FF methods. The sum of the probability-weighted contributions corresponds to the reported average values for the β-factor.

Density functional theory (DFT) calculations were carried out for a selected set of configuration using the Vienna ab-initio simulation package (VASP) (Kresse and Furthmüller, 1996; Kresse and Hafner, 1994), which expands the electronic wavefunction on a plane wave basis set and incorporates the projected augmented wave (PAW) method for pseudopotentials (Kresse and Joubert, 1999). The energy of electronic exchange and correlation was calculated using the generalized gradient approximation by the Perdew-Burke-Ernzerhof revised for solids (PBEsol) (von Barth, 2004). DFT calculations for calcite, dolomite, aragonite, and magnesite were performed using the conventional unit cell. Results
reported in figure 3 and 6 were obtained from supercell configurations sizes as reported in the Supplementary Material. An energy cutoff of 600 eV and a sampling mesh grid of at least a 4x4x4 Monkhorst-Pack grid centred at gamma point (Monkhorst and Pack, 1976). We found the total energy to converge within 0.01 eV/formula unit for all the studied structures.

For comparison, we also include DFT estimated isotope fractionation values for some selected configurations.

## 3. Results

### 3.1 Effect of Mg and Ca concentrations on the stability and structural properties of Mg-bearing calcite

The calculated structural parameters for the systems of calcite, dolomite, magnesite and aragonite estimated from the classical force field are in good agreement with the experimental values with differences less than 1% (see Table 1). The calculated vibrational frequencies, which will ultimately affect the calculation of the isotope fractionation factors, agree also well with the experimental results, as shown in Figure 2a. These results allow us to conclude that the classical forcefield selected in this work represents a good approximation for the estimation of isotope fractionation factors, at least within the quasiharmonic approximation.

The calculated volume for all carbonate minerals has been obtained considering thermal expansion effects calculated from static condition to 298K, where we see a volume increase of 1.0% compared with the 0K static values. The temperature correction for the volume values is consistent with values reported by other authors (Wang et al., 2017) as shown in figure 2b. This figure also shows the reduction of total volume with the increase in Mg concentration due to a shorter average Mg-O bond distance ($d_{Mg-O}=2.1$ Å) in comparison to the Ca-O bond distance ($d_{Ca-O}=2.4$ Å). This behaviour seems to follow a linear trend whose equation is $V_{MgxCa(1-x)CO_3}=61.065-16.691x$.
The effect of Mg concentration on calcite-like structures has been studied by preparing supercells with 1 to 4 Ca-Mg substitutions, with Mg / (Mg + Ca) (molar ratio) ranging from 1/96 to 4/16 (Figure 1). Ca atoms in a calcite structure were replaced following the method implemented by D'Arco et al. (2013); Mustapha et al. (2013) to obtain a full set of symmetrically non-equivalent configurations. Subsequently, we performed structural optimizations for all the independent structures obtained. The profiles in Figure 3 show two examples of relation between the distribution of the atoms and the total energy. To test the force field (FF) performance on the estimation of ground-state properties, we have also performed density functional theory calculations for the same sets of structures. Results show a good agreement between the FF and DFT calculations. The FF shows a similar energy profile for the concentrations studied and identifies the same minimum-energy structures as those obtained from the DFT methods. This suggests that the used FF can correctly reproduce the energetic order and stability of the cation mixing/ordering behaviour.

The total energy is the main parameter used to determine the stability of a crystalline structure at zero temperature conditions. However, temperature can also be an important parameter when discussing structure stability and isotope fractionation. Furthermore, the configurational multiplicity or equivalent structural configurations in the case of solid solutions can be also a factor to consider. Because of this, the contribution of each inequivalent configuration was weighted by its Boltzmann probability at 298.15 K (see SI for information about the multiplicity (Ωm) and the probability for each of the studied structures at a given composition). The sum of the probability weighted contributions represents the average property of the Mg-bearing carbonate at a finite temperature. The use of this probability allows us to characterize the estimated fractionation factors as a function
of temperature and energy. In this work we found that the structural configuration with the lowest total energy does not necessarily correspond to the one with the highest probability at 298.15 K (see SI). In the case of Mg-bearing calcite, the Ca-Mg substitution leads to a significant distortion in the structure, especially as the Mg concentration increases. This distortion is evident in the distribution of values obtained for the effective average Mg-O bond distance (d(Mg-O)*), especially in cases involving more than one Mg substitution. In these cases, the d(Mg-O)* distance is consistently equal to or lower in the configuration with the minimum-energy compared to the one with the highest probability (Table 2). Additionally, when comparing the Mg-Mg effective coordination number (ECN<sub>Mg-Mg</sub>), it exhibits the opposite behaviour (Table 2). Furthermore, the effective coordination number of oxygen atoms surrounding a Mg atom (ECN<sub>Mg-O</sub>) is distributed in the range 5.9 - 6.0 for almost all structures, except the 4/16 structure that shows a ECN<sub>Mg-O</sub> near to 5.8. In pure calcite, each Ca atom is surrounded by six oxygens from six monodentate CO<sub>3</sub><sup>2-</sup> ion groups. When an Mg<sup>2+</sup> substitutes a Ca<sup>2+</sup>, there is a stronger electrostatic attraction between Mg<sup>2+</sup> and O<sup>2-</sup> leading to a shorter distance of d(Mg-O)* compared to the original d(Ca-O)*. Because of the Mg<sup>2+</sup> substitution, some of the CO<sub>3</sub><sup>2-</sup> underwent a slight rotation and displacement, resulting in much shorter Mg-O distances and longer Ca-O distances depending of Mg distribution (see Figure S1).

## 3.2. Effect of Mg and Ca concentrations on the equilibrium fractionation factors

To understand how Mg concentration affects the equilibrium isotope fractionation factors of <sup>44</sup>Ca/<sup>40</sup>Ca and <sup>26</sup>Mg/<sup>24</sup>Mg, we have calculated the β-factor as a function of Mg concentration for the configuration with the lowest total energy, the configuration with the highest probability, as well as, the value obtained from a weighted average using the Boltzmann
probability (see SI for individual β-factors). Values for these three approaches are reported in Fig. 4. First, the Ca β-factors for calcite and aragonite are different and are related to the difference in structure and the Ca-O coordination number. In the case of calcite, there are six oxygens surrounding a Ca atom while there are nine in the case of pure aragonite (Table 1). Consequently, the Ca-O bond distance is larger in the case of aragonite than in calcite, leading to a lower <sup>44</sup>Ca/<sup>40</sup>Ca β-factor for the former mineral compared to the latter. Secondly, the β-factor for <sup>44</sup>Ca/<sup>40</sup>Ca calculated for dolomite is lower than the values reported for calcite (see figures 5a and 8a for value comparison). This suggests that in a similar way to stability and total energy, the fractionation factor is strongly correlated with the way Ca and Mg atoms are distributed within the structure, as well as the length of their atomic bonds.

For Mg-bearing calcite, the method used, i.e. methods based either on the lowest-energy configuration, the highest-probability configuration or on all configurations averaged according to their probability, must be considered when discussing the results. Differences are observed between the results of these methods for cases where more than one simultaneous substitution is performed, especially for concentration values beyond 2/24. Ca and Mg isotope fractionation vary differently as a function of the Mg content. Overall, we see on Figure 4a that Ca fractionation factor remains almost constant within 1.0 ‰ over the whole range of Mg content investigated here (from 0 to 25 mol%). In more detail, one can note a downward trend with increasing Mg content when considering the lowest-energy configurations, while the average β values at 298K show an upward trend.

In the case of Mg isotope fractionation, depicted in Figure 4b, the points are more scattered. The lowest-energy configurations show a general and significant upward trend with increasing Mg content (with an increase of about 5‰ over the whole compositional range).
while the average $\beta$ values at 298K show a much more moderate trend with variations within ~1‰. The differences relate to how Mg atoms are distributed in the simulation cell leading to different isotope $\beta$-factor for a given concentration. In order to further understand the effect caused by the Mg ordering on the fractionation factor we have calculated the $^{44}$Ca/$^{40}$Ca and $^{26}$Mg/$^{24}$Mg $\beta$-factor for configurations with the same Mg concentration (2/16, 3/24 and 2/24, 3/36) but with different configurational distribution and cell size. It is clear from Figure 4b that the $\beta$-factors for these concentrations are different and possibly closely related to the way atoms are distributed within the structure. Indeed, a look into the distribution of Mg atoms and its relation to the minimum configurational energy as displayed in figure 5, shows that the Mg $\beta$-factors depend on the Mg-Mg ECN (see panel 5g and 5h) and the Mg-O effective distance (panel 5e and 5f) as seen for two selected Mg concentrations. The Mg-Mg ECN shows that the structure with the largest ECN i.e. the one that contains the Mg atoms arranged close to each other as second and third neighbours as well as the lowest Mg-O distance lead to be the energetically most stable configuration possible due to the localization of any stress and structural distortion within the lattice. Thus, $\beta$-factor shows a clear correlation with the total energy (Fig. 5a and 5b) whereas the $\beta$-factor shows an anti-correlation (Fig. 5c and 5d).

Previous works by Pinilla et al. (2015) and Wang et al. (2019) have provided conflicting results. Pinilla et al. (2015) concluded that Mg-bearing calcite is enriched in the heavier Mg isotope at concentrations of 3.12 ‰mol relative to a system with a higher concentration of 6.25 mol%. A more recent computational work by Wang et al. (2019) able to explore larger concentrations, has found the opposite behaviour and attributed the results to effects due to calculations performed at constant volume that induce spurious pressure effects. Our present
study performed at constant pressure suggests that results depend on the way the atoms are distributed within the simulation cell and the energy surface explored for a given Mg composition. Our results indicate an upward trend of the Mg β-factors with increasing Mg content with a similar slope as in Wang et al (2017, 2019) when considering the lowest-energy configurations. However, the method based on the probability-weighted average at 298K suggests that the concentration effect on Mg isotope fractionation is more moderate. Given that the probability to find a determined configuration in a solid solution depends on the total energy and temperature, it suggests that for the case of solid solution a Boltzmann weighted average that takes into account this information may represent a better approach than using only the minimum-energy configurations when estimating isotope fractionation factors for these systems.

### 3.3. Effect of Ca and Mg distribution on the equilibrium isotope fractionation factor in dolomite

The CaCO₃-MgCO₃ has been one of the most studied solid solutions, with the main characteristic being an immiscibility gap separated by a small stability field in dolomite (x=0.5) with a negative enthalpy value for ordered dolomite compared to several disordered structures with the same composition (Navrotsky and Capobianco, 1987). Dolomite is a complex mineral with multiple origins: primary precipitate, diagenetic replacement or metamorphic or hydrothermal phase (Warren, 2000). Ancient sedimentary dolomites usually show highly-ordered structures with a Mg/Ca compositional ratio nearing 1, whilst modern dolomites are mostly Ca-rich with a lower degree of cationic order (Ning et al., 2019; Warren, 2000). Furthermore, dolomite cannot be synthesized in the laboratory under ambient
conditions whatever the saturation state of the solution and without bacterial mediation (Land, 1998). Current models of dolomite formation propose that, although seawater is supersaturated with respect to dolomite, the high energy barrier that controls Mg²⁺ hydration prevents its precipitation. In order to better understand these observations, experimental works employing various approaches mainly based on high-temperature mineral synthesis, have explored the evolution of the mineral composition and the cation structural order during the calcite-dolomite replacement reaction (Kaczmarek and Sibley, 2011; Navrotsky and Capobianco, 1987; Perrin et al., 2016). There is consensus on the presence of a local order in magnesian calcites, even if the nature of this order is still debated (intercalation of Ca-rich layers and Mg-rich layers, or Mg clustering in certain basal layers).

We have used the methodology described above to generate the structures of CaMg[CO₃]₂ in order to understand the ordering effect on the isotopic fractionation factor. Figure 6 shows the structures studied, where structure number 9 corresponds to ordered Dolomite. Figure 7a shows the total energy difference between the lowest-energy configuration (structure # 9) and the rest. To evaluate the force field's capacity for generating structures in the correct energy order, we present outcomes derived from both FF and DFT methods. In both approaches, the energy ranking of structures aligns closely, wherein structure #9 representing ordered dolomite exhibits the lowest formation energy. A further inspection into the effective bond distance shows that the Ca-O and Mg-O distances for ordered dolomite are the largest and smallest of all given structures, respectively (Fig. 7c). Furthermore, the effective coordination number (Fig. 7b) indicates that most of the structures show a certain number of Mg atoms surrounding a Ca atom, but this number is never bigger than in the case of ordered dolomite (six atoms). The fact that the ordered configuration number 9 corresponds to the one with the lowest energy is then due to an interplay between the Mg and Ca positions that
391 guarantees less bond strain due to CO₃ rotations, leading to a more homogeneous distribution
392 of the Mg-O and Ca-O bonds within the structure. This same result has also been pointed out
393 recently by Allan et al. (2019) whose work revealed how low enthalpy configurations are
394 associated with reduction of strain due to the rotation of the CO₃ groups, possible only when
395 carbonate lies between layers formed solely by Mg or Ca atoms.

396 The calculated β-factors for 44/40Ca and 26/24Mg isotopic ratios can be seen in Figure 8. Like
397 for Mg-bearing calcite, the β-factor trend closely resembles that of the total energy shown in
398 Figure 7. Ca and Mg β-factors are respectively correlated and anticorrelated with the
399 configurational total energy. The way in which Mg and Ca atoms are accommodated within
400 the structure affects the bonds and consequently the isotope fractionation factors. Shorter
401 Mg-O bonds lead to larger Mg β-factors and larger Ca-O bonds lead to smaller Ca β-factors.

402 The highest and lowest fractionation factors for Mg and Ca are achieved for the ordered
403 dolomite structure.

404

405 **4. Discussion**

406
407 **4.1 The effect of Mg distribution on the Mg and Ca equilibrium isotope fractionation factors in carbonates**

408 Stable isotope compositions of carbonates provide important records on many low-
409 temperature geochemical processes. Our results indicate that the fractionation of Mg and Ca
410 isotopes is sensitive to the Mg-O (Ca-O) bond distances, linked to the difference in the ionic
411 radii of Ca and Mg, as well as to the way Mg atoms are distributed within the structure. For
412 a given composition, our results show that even slight variations in the Mg arrangement lead
413 to fluctuations in the 26/24Mg β-factor of several permil at 298 K (Fig. 5). Based on DFT
calculations, Wang et al. (2019) showed that the Mg-O bond decreased sharply when the Mg concentration increased up to 16 mol% (4/24) and that the bond distance decreases smoothly beyond this value. Our work agrees with these findings. We see a decrease (with fluctuations in our results) of the effective Mg-O distance as a function of Mg concentration estimated from structures using the lowest energy and the maximum probability (Fig. 9b). In this last method the Mg-O distance may have plateau towards a distance of ~ 2.13 Å beyond 15 mol%. Likewise, the Mg-Mg effective coordination number, Fig. 9d, shows that the substitution of Mg atoms has a major distorting effect on the calcite lattice, which benefits the decrease in the Mg-O bond distance (Fig. 9b). Over the same compositional range, Ca-O bonds exhibit smaller variations but are complementary to Mg-O bonds (anti-correlation, Fig. 9b and 9c). Furthermore, the enrichment of Mg heavy isotopes in Mg-bearing calcite (Fig. 4b) is enhanced by the packing Mg atoms as second or third neighbours within the structure (Fig. 9d).

These observations at the scale of each configuration make it possible to identify the physicochemical parameters that determine the value of the local β-factor. At the solid solution scale, the approach based on the probability-weighted average shows that the effect of concentration on isotopic fractionation is much weaker than previously suggested (e.g., Wang et al. 2017, 2019). The methodology used in this work has allowed us to explore in a consistent way a larger set of configurations for each Mg concentration, which in principle is a cumbersome task to carry out using DFT methods. Because of this, it is probable that some of the structures found here do not coincide with those found as local energy minimum in the work of Wang et al. (2019). In any case, Figure 4 shows an average effect for Mg fractionation in calcite lower than the permil level at 298K. This result agrees with the
experimental study by Li et al. (2012) who concluded that there was no effect of Mg content in their experiment of inorganic precipitation of calcite containing up to 14.9 mol% MgCO₃.

In the context of the formation of Mg-calcite via the precipitation of an amorphous calcium carbonate, Mavromatis et al. (2017) observed a decrease from ~ -3.0‰ to -3.6‰ in the ²⁶Mg/²⁴Mg fractionation between calcite and reactive fluid. It was proposed that this trend was due to the simultaneous increase in the Mg content of the calcite (up to 20 mol %) and/or to the evolution towards thermodynamically more stable solid phases. Although this isotopic variation of ~ -0.6‰ at 25°C cannot be totally excluded when considering the uncertainty of the theoretical results (i.e., dispersion of data around general trends), our results suggest almost no variation in the equilibrium isotopic composition of calcite with increasing Mg content. The second hypothesis proposed by Mavromatis et al (2017), i.e. the thermodynamic evolution of the solid phases formed, therefore seems the most likely. It is important to note, however, that precipitation kinetics can influence cationic disorder and therefore vary this compositional effect.

## 4.2 Structural and isotopic properties of ordered and disordered dolomite

Dolomite shows contrasted structural and isotopic results between ordered and disordered structures. The ordered structure shows a balance between the cation distribution, the distortions and stresses created to the structure. The rotation of CO₃ groups due to an unevenly distributed series of Mg-O and Ca-O bonds leads to structures that in principle are larger in volume than the ordered dolomite. As already shown by previous studies (e.g., Liu et al., 2016; Wang et al. 2017), we find that ordered dolomite has the smallest Mg-O distance and the largest Ca-O distance from the calcite-magnesite solid solution. This order-disorder also affects the isotopic properties with ordered dolomite favouring the incorporation of the heavy ²⁶Mg isotope in comparison to other structures with similar composition, while the
opposite effect is seen for the case of <sup>44</sup>Ca (Fig. 8). The isotope fractionation disparities between ordered and disordered dolomite yield a 5‰ enrichment in the heavy Mg isotope at 25 °C. Our findings exhibit a larger but consistent trend with the experimentally observed 0.25‰ variation at 160-220°C by Li et al. (2015). This consistency emerges when comparing the results obtained from the Mg isotope fractionation of aqueous solutions with those from both poorly-ordered and ordered dolomite systems.

### 4.3. Mg isotope fractionation between calcite and dolomite

Figure 9 shows the equilibrium Mg isotope fractionation factor between calcite and dolomite and its correlation with the changes in the Mg-O bond distance in Mg-bearing calcite. When considering the lowest-energy configurations, this isotope fractionation factor follows a trend with a negative dependence on the Mg-O distances. However, it should be noted that the calcite-dolomite fractionation is almost constant for the composition range below 5.0 mol% MgCO<sub>3</sub>. Our results show that for these smallest concentrations such as 1/24 and 1/48, the variation of the 10<sup>3</sup>lnα<sub>Calcite-Dolomite</sub> fractionation factor is of the order of ~ 0.2‰ and is related to dilute substitutions that generate small disturbance within the crystal lattice.

When considering the results based on the probability-weighted average, the calcite-dolomite fractionation shows a more moderate trend with small fluctuating values, differing by less than 1.0‰ and remaining relatively constant for most of the studied concentrations. To investigate whether this behavior is temperature-dependent, we have also generated a figure similar to Figure 4 but using data obtained at a higher temperature of 433 K (Fig. SI2). We observe that, although there is a reduction in the fractionation factor as expected at higher
483 temperature, the behaviour of the β-factor as a function of Mg concentration resembles that
484 at lower temperatures for the studied methods.

485 In order to further understand the effect of volume, chemical exchange and ion position when
486 generating solid solutions Ca₁₋ₓMgₓCO₃ as already presented for calcite and dolomite (Fig. 4
487 and Fig. 8) we follow the work by Liu et al. (2016) and use the results from figure 2b to
488 estimate the deviation of the calculated volumes for the Ca₁₋ₓMgₓCO₃ solid solutions (See
489 Supplementary Material for volume values) from a linear behaviour of an ideal mixing where
490 the volume can theoretically be expressed as:

$$V_{sol}(Ca_{1-x}Mg_xCO_3) = (1-x)V(CaCO_3) + xV(CaCO_3) \quad (6)$$

491 Figure 10a displays the percentage of volume deviation from linearity for the studied
492 systems. We see that deviations for solid solutions are in general negative and lower than
493 1%. Deviations for ordered and disordered dolomite tend to be more negative stressing the
494 importance of Mg distribution in the structure. Furthermore, the total formation energy per
495 formula unit of any of the solid solution can be written using the two end members CaCO₃
496 and MgCO₃ as:

$$\Delta E_f = E(Ca_{1-x}Mg_xCO_3) - (1-x)E(CaCO_3) - xE(MgCO_3) \quad (7)$$

497 where E(Ca₁₋ₓMgₓCO₃), E(CaCO₃) and E(MgCO₃) are the total energies per formula unit
498 obtained for the solid solution, calcite and magnesite (see supplementary material for values).

500 This formation energy can be used together with the derived equations of state (EOS) for
501 the end members CaCO₃ and MgCO₃ (see Figure S3 and Table S1 in Supplementary Material
502 for fitted EOS and further information) to split the formation energy ΔEₜ in a volume
503 expansion (ΔEₜₚ) and ion exchange-relaxation (ΔEₜ₋ₚ) contributions in such a way that
504 ΔEₜ = ΔEₜₚ + ΔEₜ₋ₚ (Liu et al. 2016). ΔEₜₚ measures the effect on the total formation
505
energy resulting from volume change in the solution when compared to contributions from
the end members. $\Delta E_{xc-rlx}$ measures the effect in total formation energy of the Mg and Ca
exchange, solution cell shape and atomic position relaxation within the lattice. Figure 10b
displays $\Delta E_f$ and $\Delta E_{vc}$ for all the studied solutions, including the ordered and
disordered dolomites. We see that in the case of Ca₁ₓMgₓCO₃ solid solutions the $\Delta E_{vc}$ and
$\Delta E_{xc-rlx}$ are almost identical but of opposite sign leading therefore to small deviations in
formation energy (less than 40 meV/f.u.) for the solid solutions when compared to the end
members CaCO₃ and MgCO₃. However, in the case of ordered and disordered dolomite, it is
noticeable that all structures have similar $\Delta E_{vc}$ contributions, but the $\Delta E_{xc-rlx}$ is different,
with the ordered dolomite having a more negative contribution than the disordered ones. In
general, ion exchange produces a volume mismatch that is compensated by the CO₃²⁻ group
relaxation and distortion within the structure. However, ordered dolomite is the only case
where CO₃²⁻ groups are not distorted, thus not causing further strain on the structure (Chan
and Zunger 2009), which explains why energetically this structure is more stable than
disordered versions. This result also agrees with previous work by Liu et al. (2016) and Chan
and Zunger (2009). This final result shows that the ion exchange effects and volume
mismatch can compensate in the case of ordered dolomite leading to an isotope fractionation
factor values larger than those where these two contributions do not compensate evenly as
seen in figure 8.

**4.4 Methodological aspects for the estimation of equilibrium isotope fractionation factors in solid solutions using computer simulations**
During recent years, several studies have been carried out to understand the effect of cation concentration on the equilibrium isotope fractionation of several species. These works have been based on DFT methods calculations containing tens of atoms. Even for single-atom substitution models, discrepancies have been found in the literature (Pinilla et al., 2015; Wang et al. 2019). Son et al. (2020) performed a thorough analysis of the structure optimization scheme (fixed or relaxed volume, constrained or unconstrained site symmetry) to explain these discrepancies. Overall, the usual approach has been to explore the energy surface space for possible structures and use those with the minimum total energy to obtain isotope fractionation factors and to compare these with experiments. Our research demonstrates that this method is valid when the element of interest is statistically isolated in the crystal structure, i.e. for small concentrations close to the endmembers of the solid solution. We can roughly estimate this validity range at 5 or 10 mol% for Mg in calcite. However, the most reliable approach involves utilizing a statistical weighted average that considers all potential non-equivalent structures generated for a given concentration. This approach mimics what happens in nature when Mg and Ca are incorporated within the structure and that thermal effects can produce a distribution of configurations with the same concentration but differing atomic order whose effect is changing the fractionation factor in the permil order of magnitude in the case of Mg in calcite at room temperature. In the case of Mg-bearing calcite, using this approach gives a more moderate concentration effect compared to the approach based on lowest-energy configurations. This approach holds also for higher temperatures where we also see a similar trend for isotope fractionation values.

549

## 5. Conclusions
Equilibrium Mg and Ca isotope fractionation has been calculated among the carbonates
551

calcite, dolomite, magnesite and aragonite using classical force field and DFT methods. From
552

a methodological point of view, this study shows that the usual approach based on the lowest
553

energy configurations is valid for studying the isotopic properties of trace or minor elements,
554

i.e. when these elements can be considered isolated in the structure. However, as the
555

concentration increases, it is essential to take account of configurational disorder, with a
556

weighted average of the contributions of each symmetrically non-equivalent configuration.
557

In this way, the calculated isotopic properties more accurately reflect the cationic distribution
558

and its effect on crystal structure.
559

Our work has shown that substitution of Mg into carbonate lattice leads to an increase in the
560

Mg fractionation factor, when either the minimum energy or maximum probability structures
561

are used for the calculation. These results agree with the trends reported by Wang et al.
562

(2017); Wang et al. (2019). However, we show that considering thermal probability as well
563

as structure multiplicity leads to an almost negligible dependence of the Mg fractionation
564

factor on Mg concentration. In all approaches, it is clear that the presence of Mg causes a
565

strain in the structure resulting in larger Ca-O bonds and shorter Mg-O bonds. Those
566

configurations where the Mg incorporation is localized to a region of the structure (cluster of
567

Mg atoms) produce smaller perturbations in the crystal structure, leading to lower Mg
568

fractionation factors. In this study, we have also investigated the fractionation of Mg isotopes
569

in ordered and disordered dolomite structures. We have shown that minimization of stress
570

and disorder in relation to the distribution of Mg and Ca cations in the dolomite lattice leads
571

to the ordered structure and the highest Mg equilibrium isotope fractionation factor.
572

573

574 ACKNOWLEDGEMENTS
This work was performed using Granado-HPC from the Universidad del Norte and SCARF
from the STFC of the UK. C.P. acknowledges funding from MINCIENCIAS (Contract No.
112721-296-2022) and the Royal Society of the UK (Grant ICA-R1-201174) and thanks the
International Centre for Theoretical Physics (ICTP) for an Associate Fellowship. Finally, we
thank the associate editor and three anonymous reviewers for their comments, which have
significantly improved this manuscript.

581

## 582 APPENDIX A. SUPPLEMENTARY MATERIAL

Supplementary materials contain additional data for this work such as further information on
interatomic distances of Mg-O and Ca-O for a selected set of Mg-substituted calcite (Ca₁ₓMgₓCO₃) solid solutions, the reduced partition function ratio (β-factor) for 44/40Ca and
²⁶/²⁴Mg as a function of Mg concentration at 458 K, the equations of state for calcite and
magnesite used to estimate the volume and exchange-relaxation contributions to the
formation energy of the studied solid solutions and, lists the vibrational frequencies obtained
with the classical force field classified according to different vibration modes.

590

591

## 592 APPENDIX B. SUPPLEMENTARY MATERIAL

It contains information concerning estimated total energies, multiplicities, and probabilities
for all the studied structures. It also contains data related to the pure phases and point
calculations used to obtain the equations of state presented in this work.

596

## 597 AUTHOR CONTRIBUTIONS STATEMENT
P.P. and C.M. carried out the simulations and analysed data. A.M.E software and visualization. C.P. and M. B. conceptualization, results analysis, and supervision. All authors contributed to the writing- reviewing and editing of the manuscript.

## 6. DATA AVAILABILITY

Data used for figures and tables are available through Mendeley Data at doi: 10.17632/9ytj4s59gd.1

Alkhtib, M. and Eisenhauer, A. (2017a) Calcium and strontium isotope fractionation during precipitation from aqueous solutions as a function of temperature and reaction rate; II Aragonite. Geochim. Cosmochim. Acta 209, 320-342.

Alkhtib, M. and Eisenhauer, A. (2017b) Calcium and strontium isotope fractionation in aqueous solutions as a function of temperature and reaction rate; I Calcite. Geochim. Cosmochim. Acta 209, 296-319.

Allan, N.L., Thomas, L., Hart, J.N., Freeman, C.L. and Mohn, C.E. (2019) Calcite-magnesite solid solutions: using genetic algorithms to understand non-ideality. Phys. Chem. Miner 46, 193-202.

Barth, K., & Gomm, B. (2004) Gruppentest zur Früherkennung von Lese- und Rechtschreibschwierigkeiten: Phonologische Bewusstheit bei Kindergartenkindern und Schulanfängern (PB-LRS). Reinhardt

Blöchl, P.E. (1994). Projector augmented-wave method. Phys. Rev. B, 50(24), 17953-17979.

Bruneval, F., Donadio, D. and Parrinello, M. (2007) Molecular dynamics study of the solvation of calcium carbonate in water. J. Phys. Chem. B 111(42) 12219-12227.

D. Buhl, A., Immenhauser, G., Smeulders, L., Kabiri, D.K. and Richter (2007) Time series δ26Mg analysis in speleothem calcite: Kinetic versus equilibrium fractionation, comparison with other proxies and implications for palaeoclimate research. Chem. Geol 244, 715-729.

Chan, J. A., and Zunger, A. (2009) II-VI oxides phase separate whereas the corresponding carbonates order: The stabilizing role of anionic groups. Phys. Rev. B, 80(16), 165201.
629 D'Arco, P., Mustapha, S., Ferrabone, M., Noël, Y., De La Pierre, M. and Dovesi, R. (2013)

630 Symmetry and random sampling of symmetry independent configurations for the

631 simulation of disordered solids. J. Phys.: Condens. Matter 25(35), 355401.

632 Dana, J.D. and Hurlbut, C.S. (1960) Manual de Mineralogía. Rev. S.A 2a Edición.

633 De La Rocha, C.L. and DePaolo, D.J. (2000) Isotopic evidence for variations in the marine

634 calcium cycle over the Cenozoic. Science 289 (5482), 1176-1178.

635 Drits, V.A., Mccarty, D.K., Sakharov, B. and Milliken, K.L. (2005) New insight into structural

636 and compositional variability in some ancient excess-ca dolomite. Can. Mineral 43(4),

637 1255-1290.

638 De Villiers, J. P. R. (1971) Crystal structures of aragonite, strontianite, and witherite. Am.

639 Mineral. 56(5-6), 758-767.

640

641 Effenberger, H., Mereiter, K., & Zemann, J. (1981) Crystal structure refinements of

642 magnesite, MgCO₃, and calcite, CaCO₃, with discussion of some aspects of the

643 stereochemistry of calcite-type carbonates. Z. Kristallogr. 156(1-4), 233-243.

644

645 Eisenhauer, A., Kisakürek, B. and Böhm, F. (2009) Marine Calcification: An Alkali Earth

646 Metal Isotope Perspective. Elements 5(6), 365-368.

647 Fantle, M.S. and Higgins, J. (2014) The effects of diagenesis and dolomitization on Ca and

648 Mg isotopes in marine platform carbonates: Implications for the geochemical cycles of Ca

649 and Mg. Geochim. Cosmochim. Acta 142, 458-481.

650 Farkaš, J., Böhm, F., Wallmann, K., Blenkinsop, J., Eisenhauer, A., Van Geldern, R. and

651 Veizer, J. (2007) Calcium isotope record of Phanerozoic oceans: Implications for chemical

652 evolution of seawater and its causative mechanisms. Geochim. Cosmochim. Acta 71(21),

653 5117-5134.

654 Fletcher, R. (1987) Practical methods of optimization. John Wiley & Sons.

655 Gale, J.D. (1997) GULP: A computer program for the symmetry adapted simulation of

656 solids. J. Chem. Soc., Faraday Trans 93(4), 629.

657 Gale, J.D. and Rohl, A.L. (2003) The general utility lattice program (GULP). Mol. Simul

658 29(5), 291-341.

659 Gao, C., Cao, X., Liu, Q., Yang, Y., Zhang, S., He, Y., Tang, M. and Liu, Y. (2018) Theoretical

660 calculation of equilibrium Mg isotope fractionations between minerals and aqueous

661 solutions. Chem. Geol 488, 62-75.
662 Geske, A., Goldstein, R.H., Mavromatis, V., Richter, D.K., Buhl, D., Kluge, T., John, C.M. and 663 Immenhauser, A. (2015) The magnesium isotope (δ26Mg) signature of dolomites. 664 Geochim. Cosmochim. Acta 149, 131-151.

665 Grzechnik, A., et al. (1999) Raman spectroscopy of magnesite and dolomite: A 666 comparative study. Am. Mineral. 84(7-8), 1074-1083

667

668 Gussone, N., Eisenhauer, A., Heuser, A., Dietzel, M., Bock, B., Böhm, F. and Nägler, T.F. 669 (2003) Model for kinetic effects on calcium isotope fractionation (δ44Ca) in inorganic 670 aragonite and cultured planktonic foraminifera. Geochim. Cosmochim. Acta 67(7), 1375- 671 1382.

672 Gussone, N., Böhm, F., Eisenhauer, A., Dietzel, M., Heuser, A., Teichert, B.M.A., Reitner, J., 673 Wörheide, G. and Dullo, W.C. (2005) Calcium isotope fractionation in calcite and 674 aragonite. Geochim. Cosmochim. Acta 69(18), 4485-4494.

675 Harrison, A. L., Heuser, A., Liebetrau, V., Eisenhauer, A., Schott, J., and Mavromatis, V. 676 (2023) Equilibrium Ca isotope fractionation and the rates of isotope exchange in the 677 calcite-fluid and aragonite-fluid systems at 25° C. Earth. Planet. Sci. lett, 603, 1179-603, 117985.

678

679 Heuser, A., Eisenhauer, A., Böhm, F., Wallmann, K., Gussone, N., Pearson, P.N., Nägler, T.F. 680 and Dullo, W.C. (2005) Calcium isotope (δ44/40Ca) variations of Neogene planktonic 681 foraminifera. Paleooceanography 20(2).

682 Higgins, J.A. and Schrag, D.P. (2010) Constraining magnesium cycling in marine sediments 683 using magnesium isotopes. Geochim. Cosmochim. Acta 74(17), 5039-5053.

683

684 Hippler, D., Buhl, D., Witbaard, R., Richter, D.K. and Immenhauser, A. (2009) Towards a 685 better understanding of magnesium-isotope ratios from marine skeletal carbonates. 686 Geochim. Cosmochim. Acta 73(20), 6134-6146.

686

687 Hoppe, R. (1979) Effective coordination numbers (ECoN) and mean fictive ionic radii 688 (MEFIR). Z. Kristallogr. Cryst. Mater 150(1-4), 23-52.

689 Immenhauser, A., Buhl, D., Richter, D., Niedermayr, A., Riechelmann, D., Dietzel, M. and 690 Schulte, U. (2010) Magnesium-isotope fractionation during low-Mg calcite precipitation in 691 a limestone cave – Field study and experiments. Geochim. Cosmochim. Acta 74(15), 4346- 692 4364.

692

693 Kaczmarek, S.E. and Sibley, D.F. (2011) On the evolution of dolomite stoichiometry and 694 cation order during high-temperature synthesis experiments: An alternative model for the 695 geochemical evolution of natural dolomites. Sediment. Geol 240(1-2), 30-40.

696 Kresse, G., & Hafner, J. (1994). Ab initio molecular-dynamics simulation of the liquid- 697 metal-amorphous-semiconductor transition in germanium. Phys. Rev. B, 49(20), 14251.

697
Kresse, G. and Furthmüller, J. (1996) Efficient iterative schemes for ab initio total-energy calculations using a plane-wave basis set. Phys. Rev. B, 54(16), 11169-11186.

Kresse, G. and Joubert, D. (1999) From ultrasoft pseudopotentials to the projector augmented-wave method. Phys. Rev. B, 59(3), 1758-1775.

Land, L.S. (1998) Failure to Precipitate Dolomite at 25 °C from Dilute Solution Despite 1000-Fold Oversaturation after 32 Years. Aquat. Geochem 4, 361-368.

Lemarchand, D., Wasserburg, G and Papanastassiou, D. (2004) Rate-controlled calcium isotope fractionation in synthetic calcite. Geochim. Cosmochim. Acta 68(22), 4665-4678.

Li, W., Chakraborty, S., Beard, B.L., Romanek, C.S. and Johnson, C.M. (2012) Magnesium isotope fractionation during precipitation of inorganic calcite under laboratory conditions. Earth Planet. Sci. Lett 333, 304-316.

Li, W., Beard, B.L., Li, C., Xu, H. and Johnson, C.M. (2015) Experimental calibration of Mg isotope fractionation between dolomite and aqueous solution and its geological implications. Geochim. Cosmochim. Acta 157, 164-181.

Liu, Z. T. Y., Burton, B. P., Khare, S. V., and Sarin, P. (2016) First-principles phase diagram calculations for the carbonate quasibinary systems CaCO3-ZnCO3, CdCO3-ZnCO3, CaCO3-CdCO3 and MgCO3-ZnCO3. Chem. Geol 443, 137-145.

Markgraf, S. A., & Reeder, R. J. (1985) High-temperature structure refinements of calcite and magnesite. Am. Mineral. 70(5-6), 590-600.

Markgraf, S. A., & Reeder, R. J. (1986) High-temperature crystal chemistry of dolomite. Am. Mineral. 71(5-6), 795-804.

Mavromatis, V., Meister, P. and Oelkers, E.H. (2014) Using stable Mg isotopes to distinguish dolomite formation mechanisms: A case study from the Peru Margin. Chem. Geol 385, 84-91.

Mavromatis, V., Gautier, Q., Bosc, O. and Schott, J. (2013) Kinetics of Mg partition and Mg stable isotope fractionation during its incorporation in calcite. Geochim. Cosmochim. Acta 114, 188-203.

Mavromatis, V., González, A.G., Dietzel, M. and Schott, J. (2019) Zinc isotope fractionation during the inorganic precipitation of calcite – Towards a new pH proxy. Geochim. Cosmochim. Acta 244, 99-112.
732 Mavromatis, V., Immenhauser, A., D. Buhl, A., Purgstaller, B., Baldermann, A. and Dietzel, M (2017) Effect of organic ligands on Mg partitioning and Mg isotope fractionation during low-temperature precipitation of calcite in the absence of growth rate effects. Geochim. Cosmochim. Acta 207, 139-153.

736 Méheut, M., Lazzeri, M., Balan, E. and Mauri, F. (2007) Equilibrium isotopic fractionation in the kaolinite, quartz, water system: Prediction from first-principles density-functional theory. Geochim. Cosmochim. Acta 71, 3170-3181.

739 Monkhorst, H. J. and Pack, J.D. (1976) Special points for Brillouin-zone integrations. Phys. Rev. B, 13(12), 5188-5192.

741 Moynier, F. and Fujii, T. (2017) Calcium isotope fractionation between aqueous compounds relevant to low-temperature geochemistry, biology and medicine. Sci. Rep 7(1), 44255.

744 Mustapha, S., D'Arco, P., De La Pierre, M., Nöel, Y., Ferrabone, M. and Dovesi, R. (2013) On the use of symmetry in configurational analysis for the simulation of disordered solids. J. Phys.: Condens. Matter 25(10), 105401.

747 Navrotsky, A. and Capobianco, C. (1987) Enthalpies of formation of dolomite and of magnesian calcites. Am. Mineral 72(7-8), 782-787.

749 Ning, M., Huang, K., Lang, X., Ma, H., Yuan, H., Peng, Y. and Shen, B. (2019) Can crystal morphology indicate different generations of dolomites? Evidence from magnesium isotopes. Chem. Geol 516, 1-17.

752 Perdew, J.P., Ruzsinszky, A., Csonka, G.I., Vydrov, O.A., Ireta, J. and Scuseria, G.E. (2008) Restoring the density-gradient expansion for exchange in solids and surfaces. Phys. Rev. Lett. 100(13), 136406.

755 Perrin, J., Vielzeuf, D., Laporte, D., Ricolleau, A., Rossman, G.R. and Floquet, N. (2016) Raman characterization of synthetic magnesian calcites. Am. Mineral 101(11), 2525-2538.

757 Pilati, T. A., et al. (1998). High-temperature vibrational spectra of dolomite, calcite, and magnesite: A comparative study. Am. Mineral. 83(5-6), 689-696.

760 Pinilla, C., Blanchard, M., Balan, E., Natarajan, S.K., Vuilleumier, R. and Mauri, F. (2015) Equilibrium magnesium isotope fractionation between aqueous Mg2+ and carbonate minerals: Insights from path integral molecular dynamics. Geochim. Cosmochim. Acta 163, 126-139.

764 Poplavko, Y.M. (2019) Structure of electronic materials. Electro. Mater, 1-69.
765 Purgstaller, B., Mavromatis, V., Immenhauser, A. and Dietzel, M. (2016) Transformation of Mg-bearing amorphous calcium carbonate to Mg-calcite – In situ monitoring. *Geochim. Cosmochim. Acta* 174, 180-195.

768 Reeder, R. J., & Markgraf, S. A. (1986) High-temperature crystal chemistry of dolomite. *Am. Mineral.* 71(5-6), 795-804.

771 Richet P., Bottinga Y. and Javoy M. (1977) A Review of Hydrogen, Carbon, Nitrogen, Oxygen, Sulphur, and Chlorine Stable Isotope Fractionation Among Gaseous Molecules. *Ann. Rev. Earth Planet. Sci.* 5, 65-110.

774 Rustad, J.R., Casey, W.H., Yin, Q.Z., Bylaska, E.J., Felmy, A.R., Bogatko, S.A., Jackson, V.E. and Dixon, D.A. (2010) Isotopic fractionation of Mg2+(aq), Ca2+(aq), and Fe2+(aq) with carbonate minerals. *Geochim. Cosmochim. Acta* 74(22), 6301-6323.

777 Saenger, C. and Wang, Z. (2014) Magnesium isotope fractionation in biogenic and abiogenic carbonates: Implications for paleoenvironmental proxies. *Quat. Sci. Rev* 90, 1-21.

780 Schauble, E.A. (2011) First-principles estimates of equilibrium magnesium isotope fractionation in silicate, oxide, carbonate and hexaaquamagnesium(2+) crystals. *Geochim. Cosmochim. Acta* 75(3), 844-869.

783 Schott, J., Mavromatis, V., Fujii, T., Pearce, C.R. and Oelkers, E.H. (2016) The control of carbonate mineral Mg isotope composition by aqueous speciation: Theoretical and experimental modeling. *Chem. Geol* 445, 120-134.

786 Son, S., Li, W., Lee, J.Y. and Kwon, K.D. (2020) On the coordination of Mg2+ in aragonite: Ab-initio absorption spectroscopy and isotope fractionation study. *Geochim. Cosmochim. Acta* 286, 324-335.

789 Tang, J., Dietzel, M., Böhm, F., Köhler, S.J. and Eisenhauer, A. (2008) Sr2+/Ca2+ and 44Ca/40Ca fractionation during inorganic calcite formation: II Ca isotopes. *Geochim. Cosmochim. Acta* 72(15), 3733-3745.

792 Wang, W., Qin, T., Zhou, C., Huang, S., Wu, Z. and Huang, F. (2017) Concentration effect on equilibrium fractionation of Mg-Ca isotopes in carbonate minerals: Insights from first-principles calculations. *Geochim. Cosmochim. Acta* 208 185-197.

795 Wang W., Zhou C., Liu Y., Wu Z. and Huang F. (2019) Equilibrium Mg isotope fractionation among aqueous Mg 2+, carbonates, brucite and lizardite: Insights from first-principles molecular dynamics simulations, *Geochim. Cosmochim. Acta*. Elsevier Ltd, pp. 117-129.
799 Warren, J. (2000) Dolomite: occurrence, evolution and economically important associations. Earth. Sci. Re 52(1-3), 1-81.

801 Weir, F. W., & Lippincott, E. R. (1961). Vibrational frequencies of calcite, aragonite, and other minerals in the carbonate group. J. Chem. Phys. 35(3), 890-897.

803

804 Wilson, E.B., Decius, J.C. and Cross, P.C. (1980) Molecular vibrations: the theory of infrared and Raman vibrational spectra. Courier Corporation.

806 Ye, Y., Smyth, J. R., & Boni, P. (2012) Crystal structure and thermal expansion of aragonite-group carbonates by single-crystal X-ray diffraction. Am. Mineral. 97(4), 707-712.

808
<figure><img src="image_2.png" /><figcaption>**Figure 1.** Lowest-energy configurations of calcite, aragonite, magnesite, dolomite and Mg-bearing calcite. Ca atoms are represented in blue, Mg atoms in orange, O atoms in red and C atoms in dark brown. All supercells have been obtained using the VESTA code.</figcaption></figure>
<figure><img src="image_3.png" /><figcaption>**Figure 2.** a) Calculated vibrational frequencies using the classical force field for calcite, dolomite, magnesite and aragonite in comparison with experimental values reported by several works (Weir and Lippincott 1961, Grzechnik et al 1999), b) calculated volume using the classical force field as a function of Mg concentration for the studied carbonates minerals (orange). We also include values from Effenberger et al 1981 (blue), Markgraf, S. A., & Reeder, R. J. 1986 (darkgreen), Pilati et al. 1998 (Magenta), Bruneval et al. 2007 (violet), Wang et al. 2017 (red) for comparison. The straight line corresponds to a linear fitting of our calculated values whose equation is $V_{MgXCa1-xCO3}=61.065-16.691x$.</figcaption></figure>
<figure><img src="image_4.png" /></figure>

<figure><img src="image_5.png" /><figcaption>**Figure 3.** a) and b) Structural total energy for the Mg-bearing calcite with 0.187 (3/16) and 0.125 (2/16) mol% respectively. DFT (green-triangle-up) and FF (blue-circle) methods show that the FF can predict correctly the stability surface energy. The total energy is given relative to the lowest-energy configuration.</figcaption></figure>
<figure><img src="image_6.png" /></figure>

<figure><img src="image_7.png" /></figure>

**Figure 4.** Reduced partition function ratio (β-factor) for <sup>44/40</sup>Ca (top) and <sup>26/24</sup>Mg (bottom) as a function of Mg concentration at 298K. Blue circles correspond to β-factor estimated using the lowest energy configuration while green tri-fine are the β-factor obtained from the highest probability structure obtained using a Boltzmann distribution. Purple triangle-right symbols correspond to the β-factor obtained from a probability-weighted average. Light green star corresponds to β-factor estimated in this work using DFT. We also include values from Pinilla et al. 2015 (square), Wang et al. 2017 (diamond), Wang et al. 2019 (triangle-down) and Song et al. 2020 (triangle-up) for comparison.
<figure><img src="image_8.png" /><figcaption>**Figure 5.** Panels a) and b) correspond to the total configurational energy and the Ca β-factor estimated for all the symmetrically inequivalent structures at 0.187 mol% (3/16) and 0.125 mol% (3/24) Mg concentration. Panels c) and d) correspond to the same plot for Mgβ-factor. Panels e) and f) show effective Mg-O bond distance for the same concentrations. Panels g) and h) show the corresponding Mg-Mg effective coordination number for these systems.</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>**Figure 6:** Ordered (9) and disordered (1-8, 10) structures with the stoichiometric composition of dolomite CaMg[CO₃]₂ with different Mg and Ca atoms distributions. Corresponding space groups are specified in brackets. The blue and orange colours correspond to Ca and Mg octahedra respectively. Brown atoms correspond to carbon and the red atoms oxygen. The black lines show the rhombohedral unit cell. The hexagonal c-axis direction (stacking direction) corresponds to the long diagonal of the rhombohedral cell.</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>**Figure 7.** a) Energy difference for all dolomite configuration (Fig. 6) with respect to the ordered configuration (number 9) as found from the classical FF and DFT calculations. b) correspond to Ca-Mg effective coordination number and c) are the Ca-O and Mg-O effective bond distances.</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>**Figure 8.** a) Reduced partition function ratio (β-factor) for <sup>44/40</sup>Ca and b) <sup>26/24</sup>Mg as a function of studied dolomite configurations as found from the classical FF at 298K. Ordered dolomite corresponds to configuration 9 in these plots. The star symbols represent the β-factor values estimated for ordered dolomite using DFT. We also include values from Pinilla et al. 2015 (square), Wang et al. 2017 (pentagon), Wang et al. 2019 (triangle-down), Rustad et al. 2010 (diamond), Schauble et al. 2011 (cross).</figcaption></figure>
<figure><img src="image_12.png" /><figcaption>**Figure 9.** a) the equilibrium Mg isotope fractionation (10³ln(α)) between dolomite and Mg-substituted calcite (Ca₁-xMgxCO₃) at 298 K. Blue circles correspond to values estimated using the lowest energy configuration. Green circles are those obtained from the highest probability structure obtained using a Boltzmann distribution. Purple triangle-right symbols correspond to the β-factor obtained from a probability-weighted average. We also include values from Wang et al. 2017 (diamond) and Wang et al. 2019 (triangle-down) for comparison. Panels b) and c) correspond to effective Mg-O and Ca-O distances and panel d) is the effective Mg-Mg coordination number as a function of concentration for the solid solution of calcite type carbonates.</figcaption></figure>
<figure><img src="image_13.png" /></figure>

<figure><img src="image_14.png" /><figcaption>**Figure 10.** a) Volume deviation from linearity for all the studied Ca<sub>1-x</sub>Mg<sub>x</sub>CO<sub>3</sub> solid solutions (red tringles). Green circles and orange filled square correspond to data for disordered and ordered dolomite respectively. b) volume change energies (triangles) and ion exchange-relaxation energies (crosses) for Ca<sub>1-x</sub>Mg<sub>x</sub>CO<sub>3</sub> and disordered dolomite. Circles correspond to the total formation energy per formula unit as described by equation (7). The values of the ordered dolomite are given by filled squares.</figcaption></figure>