Modelling silicate mineral interfaces for carbon dioxide sequestration: structure and stability of orthoenstatite surfaces

Abstract

The atomic structure and stability of the surfaces of orthoenstatite (Mg₂Si₂O₆), a silicate often present in natural environments where carbon dioxide is captured and stored into minerals or reduced, have been investigated for the first time using density functional theory and a classical rigid-ion force field model. Only non-polar stoichiometric surfaces have been considered, with a specific focus on the {210}, {010}, {100}, {110} and {120} surfaces that have been most frequently reported in the literature. It turns out that there are only three main surfaces ({210}, {010}, {100}), all the others being microfacets (i.e. surfaces composed of multiple crystalline facets) exhibiting combinations of the three possible surfaces. The {100} surface, which is the main one reported in the literature, is the most unstable and undergoes surface reconstruction due to major rearrangement of magnesium atoms and their coordination shell. Results suggest that this surface would either be hydroxylated or undergo other forms of reconstruction altering the Si–O bonding pattern (i.e. forming a phyllosilicate layer) depending on the chemical environment. This work provides the background knowledge and model that can be used to investigate CO₂ reduction at the enstatite-gas interfaces and to develop future models able to describe the enstatite-liquid interface in the context of CO₂ capture.

---

1 Introduction

The excessive release of CO₂ into the atmosphere is a problem plaguing our society, due to its direct link to global warming and climate change.^1,2 Accelerating the removal of CO₂ from the atmosphere through carbon dioxide capture techniques, together with limiting its emissions by transitioning to carbon pollution-free energy, play fundamental roles in tackling this challenge.

Among the several promising methodologies and technologies that are being developed (e.g. Berstad et al.,³ Goli et al.,⁴ Sreedhar et al.,⁵ Gao et al.,⁶ Zhang et al.,⁷ Todorova et al.⁸), the natural processes occurring at the interface between silicate-rich minerals and aqueous fluids are extremely interesting due to their demonstrated potential to address both the removal of carbon dioxide and the production of clean fuel. In fact, depending on the conditions and chemical environment, these mineral phases can undergo either mineral carbonation (i.e. dissolve and precipitate as carbonates), or catalyse the reduction of CO₂ to organic molecules.⁹ The former process has already found its first large-scale application,¹⁰ although challenges related to costs, safety, and overall feasibility remain at other locations.^11–14 The latter process has received more limited attention in this context, while attracting much interest in the fields of organic matter in space, hydrothermal systems, and early Earth conditions.^9,15–17

Due to the interdisciplinary and complex nature of these processes and their direct implications being predominantly within the domain of the geosciences (element cycling, geophysical processes, origin and support of life^14,17–19), the fundamental understanding of the chemical reactions involved with carbon dioxide sequestration and reduction is very limited. This includes knowledge of the atomic details of the structure, morphology, and reactivity of the mineral interfaces, which are the places where carbon dioxide binds to either form new minerals or react.

Within this context, MgFeSiO₄ olivine and Mg₂SiO₄ forsterite are the most widely studied phases. While many fundamental details of their interface with liquid water and CO₂ are still a matter of investigation, much progress has been achieved in the definition of their surface structure, morphology and adsorption properties.^20–26 In contrast, the other silicate phases that are present in the CO₂ sequestration and reduction environments have not received much attention to date.

In this paper, we investigate the surface structure and stability of orthoenstatite, a common rock forming silicate that is abundant in both the terrestrial and extra-terrestrial environments where the processes described above occur.^9,27 Orthoenstatite (Pbca) is the orthorhombic Mg end-member of the pyroxene group XYSi₂O₆ (where X and Y are divalent cations). It is a single-chain inosilicate, where the chains are made of SiO₄ tetrahedra with a Si : O ratio of 1 : 3. Two six-fold coordination cation sites are present: an irregular octahedral site that is able to accommodate intermediate- to large-size cations, and a more regular octahedral site able to accommodate intermediate- to small-size cations. The structure of orthoenstatite is well known (see, for example, Ghose et al.,²⁸ Demichelis et al.²⁹), and so are many properties of the crystalline phase (e.g. Demichelis et al.,²⁹ Angel and Hugh-Jones,³⁰ Zucker and Shim,³¹ Stangarone et al.³²). Studies on the dissolution rates of orthoenstatite along different hkl directions are also available.³³ However, to the best of our knowledge, its surface structure, morphology, and stability have so far been completely unexplored at the atomic level.

Due to the lack of experimental results and computational predictions, first principle calculations based on density functional theory have been run on a subset of orthoenstatite surfaces to provide reference data. The latter were used to validate a rigid-ion force field model that has been developed and used to study the atomic structure and stability of a larger set of non-polar stoichiometric surfaces in vacuum.

2 Computational methods

2.1 Force field development

All classical lattice dynamics calculations were run using the General Utility Lattice Program^34,35 (GULP). For the solid, geometry optimisation and phonon calculation were run using the default options and thresholds and a shrinking factor of 8. For the surfaces, geometry optimisation was run at constant volume and also using the default setup. The search for non-polar surfaces was performed using the GDIS package.³⁶

A two-region strategy is used to model the surfaces, with region 1 containing the atomic layers that exhibit a significant relaxation (the surface) and region 2 containing the atomic layers that preserve the same structure of the bulk (the bulk). A size of 1 to 4 layers (depending on the surface area) for region 1 and region 2 was sufficient to achieve convergence on the surface energy at the second decimal digit of the considered surfaces.

Here, the total energy of the system (U_system) and of the surface (U_surface) can be written as:

U_system = U₁₁ + U₁₂ + U₂₂ (1)

U_surface = U₁₁ + (1/2)U₁₂ (2)

with U₁₁ and U₂₂ being the energies of regions 1 and 2, and U₁₂ the interaction between the two regions.

The surface energy γ is then obtained through:

where U_bulk is the energy of the bulk and A is the area of the unit cell.

The potential model developed by Kerisit et al.²⁴ for Mg₂SiO₄ forsterite was initially used, which is a modified version of ClayFF³⁷ associated with SPC/E water³⁸ and EPM2 carbon dioxide.³⁹ This model has the advantage of including all the interactions needed to model the interface with liquid water and carbon dioxide. However, due to the different nature of the silicate units, the original charge of −1.525|e| on oxygen was modified through assigning all oxygen atoms in orthoenstatite the same atomic type and an average charge of −1.37|e|. We also considered assigning the two different oxygen atoms in orthoenstatite two distinct atomic types with charges of −1.4 and −1.3|e|, but found that the results were very similar, hence we report only those obtained in the former case.

Table 1 shows that Kerisit et al.²⁴'s model is not transferable as such to orthoenstatite, as it provides differences by as much as 10% on the Mg–O distances in the irregular octahedral site.

Table 1 Experimental and calculated structural properties of orthoenstatite

	Exp.^41,42	DFT	FF	FF – Kerisit et al.^24
Calculated values (DFT and force field – FF) are reported as relative differences. V, a, b, c are the unit cell volume [Å^3] and lattice parameters [Å], Mg–O and Si–O are bond distances (maximum and minimum) [Å].
V	831.98	+0.4	+0.0	+1.3
a	18.225	+0.0	−1.3	−0.6
b	8.8129	+0.2	+0.3	+0.0
c	5.18	+0.2	+1.0	+1.9
Mg–Omin	1.992	−0.1	+0.8	+2.4
Mg–Omax	2.445	+0.1	−1.7	−10.4
Si–Omin	1.585	+0.7	−0.1	+2.2
Si–Omax	1.679	+0.4	+0.4	−0.3

---

A new model was developed through refitting the inter-atomic interactions in Kerisit et al.²⁴ using the Buckingham potentials and re-defining the charges of the two oxygen types, while maintaining the charges on the other atoms unchanged. All structural parameters obtained with this model, including bond distances, are in very good agreement with the experimental values.

The re-definition of the oxygen charges means that the parameters between the mineral oxygen atoms and the SPC/E water and EPM2 carbon dioxide models would need refitting with respect to those provided in Kerisit et al.²⁴ This is outside of the scope of this paper, which is focused on the mineral surfaces, the first step towards a longer term goal of building a broader, comprehensive model able to describe the silicates-fluid interface with a similar level of accuracy that has been achieved for other systems, e.g. the carbonate–water interface.⁴⁰ Indeed, while classical force fields are unable, by definition, to describe events that involve chemical reactivity, they can be very accurate to describe the structure and dynamics of fluids at the interface with minerals, which include physisorption, non-reactive dissolution and binding.

The rigid-ion (i.e. point charges with fixed positions) nature of both the models mentioned above does not allow to reproduce other bulk properties (e.g. mechanical, elastic and vibrational properties). While a shell potential model would provide a better estimation of bulk properties, the high computational cost would limit its usefulness in modelling the orthoenstatite-fluid interface. Therefore, we suggest to limit the use of this force field to model surface properties and interfaces, and discourage its use to investigate bulk properties.

The full set of force field parameters developed in this work is reported as ESI.†

2.2 Density functional theory (DFT) calculations

All quantum mechanical calculations were performed with CRYSTAL23.^43,44 The PBE0 ⁴⁵ hybrid functional and Gaussian-type basis set originally optimised for pyrope⁴⁶ and then used for similar calculations involving silicate minerals and their surfaces,^25,29,47,48 were adopted.

Thresholds on the accuracy of the Coulomb and Hartree–Fock (HF) exchange series and on the bipolar approximation of Coulomb bielectronic integrals were tightened with respect to default values (to TOLINTEG 8 8 8 8 16 and BIPOLA 18 18). The default DFT integration grid ensured a satisfactory accuracy of the numerically integrated electron charge density (the error is on the order of 10⁻⁴|e| on a total of 800 to 2400|e|, depending on the system). The reciprocal space was sampled according to a Monkhorst–Pack mesh⁴⁹ with shrinking factor 6, corresponding to 64 and 16 to 20 k points in the first irreducible Brillouin zone of the solid and the surfaces, respectively. The threshold on the Self Consistent Field (SCF) energy was set to 10⁻⁸ hartree, with initial direct inversion of the iterative subspace (DIIS).^50,51 To achieve SCF convergence a Fock/Kohn–Sham matrix mixing between subsequent SCF cycles of 30% was selected for most calculations, and was increased to 50% for some surfaces.

Structures were optimized by using the analytical energy gradients with respect to either atomic coordinates (for the surfaces) or atomic coordinates and lattice parameters (for the solid) within a quasi-Newton scheme, combined with the Broyden–Fletcher–Goldfarb–Shanno (BFGS)^52–55 scheme for Hessian updating. Convergence was checked on energy, gradient components and nuclear displacements using default thresholds.

A slab was used to model the surfaces. Surface energy (γ) at T = 0 K was obtained from the equation

where n is the number of the d_hkl layers; U_slab(n) and U_bulk are the energy of an n − layer slab and of the bulk, respectively; A is the area of the surface unit cell; the factor 2 accounts for the upper and lower surfaces of the slab model. U_s(n) is the energy per unit area required to form the surface, which converges to its value as more d_hkl layers are added. The thickness of the slab (i.e. n) was set so that the difference between U_s(n) and U_s(n − 1) is converged at the second decimal digit. Similarly to what observed for the classical calculations, the number of layers to achieve convergence for the considered surfaces is n = 1 to n = 3 (see Table 2). DFT calculations were run with 1 additional layer with respect to what reported in the table. Despite large size unit cells (320 atoms for n = 4), several surfaces maintain some symmetry elements, lowering the number of irreducible atoms and contributing to making these calculations feasible.

Table 2 The surfaces of enstatite examined in this study, ordered by increasing surface energy

{hkl}	d hkl		γ		A		n	Microfacets
	FF	DFT	FF	DFT	FF	DFT
A line separates surfaces that have been mentioned more often in the literature from those that have rarely or never been observed. Miller indices ({hkl}), distance between planes (dhkl [Å]), surface energy (γ [J m^−2]) and surface area (A [Å^2]) are reported for each set of surfaces, as calculated with the force field (FF) and DFT. The minimum number of layers needed to achieve convergence on the surface energy (n) are also reported. Note that n < 3 has not been considered for surfaces with A > 200 Å^2. Orientations of the microfacets are shown in the last column.
{210}	6.306	6.343	1.56	0.95	131.94	131.73	2	—
{010}	8.840	8.831	1.80	1.03	94.11	94.61	2	—
{110}	7.934	7.948	1.87	1.08	104.85	105.13	2	{210}, {010}
{120}	4.293	4.291	1.92	1.10	193.81	194.70	3	{010},{100}
{100}1	17.994	18.232	2.48	1.30 (1.37)	46.23	45.83	1	—
{100}2	17.994	18.232	2.60	1.25 (1.32)	46.23	45.83	1	—
{100}3	17.994	18.232	2.98	1.35 (1.40)	46.23	45.83	1	—
{320}	3.558	3.572	1.72	—	233.80	233.90	3	{210},{010}
{430}	2.465	2.473	1.77	—	337.50	337.90	3	{210}, {010}
{130}	2.908	2.906	1.89	—	286.08	287.52	3	{010}, {100}
{230}	2.800	2.801	1.93	—	297.08	298.28	3	{010},{210}
{310}	4.963	5.006	2.25	—	167.62	166.90	4	{100},{010}
{410}	4.010	4.050	2.34	—	207.51	206.29	4	{100},{010}

---

The counterpoise correction scheme proposed by Boys and Bernardi⁵⁶ was used to estimate the Basis Set Superposition Error (BSSE), following a similar procedure as in De La Pierre et al.⁵⁷ and in Demichelis et al.²⁵ In particular, for each hkl direction, the slab was cut from the bulk and layers of ghost atoms were added to both surfaces of the slab, with atomic functions centred at the atomic positions of the atoms in the bulk. A single point energy calculation was then run for each slab with progressive increase of ghost layers, until convergence on surface energy was reached (10⁻² J m⁻²). The error on the surface energy calculated as such was then subtracted from the surface energy of the optimised slab.

3 Results and discussion

Enstatite started attracting the attention of mineralogists nearly two hundred years ago, when it was first identified and then found in several locations, as well as in meteorites. Here, a rich literature on its possible crystal forms is available, reporting more or less elongated prisms exhibiting mostly the {100}, {010}, {210} faces, and sometimes {110} and {120} faces, in different proportions. The shallow vaulted summits of the prisms tend to include a great number of repeated faces with small inclinations, i.e. {112}, {122}, {121}, {124}, {324}, {111}, {108}, {104}, {101} and, more generally, {h0l}.^58–64 Crystal models, extracted by the Mindat database⁶⁵ team from literature that we were unable to access also exhibit the {310}, {410}, {102}, {211}, {221} and {331} cleavage.

It is important to highlight that this information is of qualitative nature, mostly relying on physical measurements of their geometry using a goniometer on natural samples containing both Mg and Fe in different proportions (most samples are actually hypersthene), often associated with clinopyroxene and containing impurities (e.g. Al and water). To the best of our knowledge, there is no additional study that elucidates the surface structure, stability, and morphology of orthoenstatite.

Eleven surfaces were identified in this study and are reported in Table 2. They correspond to all the possible non-polar surfaces from d_hkl = a to d_hkl = 2.5 Å that preserve Si–O bonding. All surfaces cutting along the z direction would break one or more SiO₄ units and require surface hydroxylation or major reconstruction. Although we cannot exclude that these could be possible, especially when in contact with fluids, we have assumed that the breaking of a strong, stable Si–O bond should be energetically disfavoured with respect to other cuts where all covalent bonds are preserved. In fact, the crystal shapes of enstatite described in the aforementioned literature suggest that most of the crystal surface is represented by {hk0} surfaces. Additionally, studies on dissolution of pyroxenes⁶⁶ show that even in extreme pH conditions and high temperatures, the Mg–O bonds are the first to break, without altering, at least initially, the structure of the Si–O framework. The dissolution would then proceed with a relatively slow detachment of silica from partially liberated tetrahedral chains. Therefore, in this initial phase of model development, we have disregarded surfaces with a Miller index l ≠ 0.

In order to preserve the Si–O bonding structure, many {hk0} surfaces turn out to be rough and microfaceted. In particular, there are three main surfaces: {210}, {100}, and {010}, which are also the ones most often reported in the literature. All the other surfaces are microfacets exhibiting two of these three orientations. The first part of the table reports the {hk0} surfaces that we have found in the literature.^{58–61,63,64} The second part of the table reports surfaces that are in principle also possible, but have never been observed.

Given that many of these surfaces have large unit cell areas, DFT calculations were performed only on a subset of surfaces, which includes the three main surfaces and the two most commonly reported microfacets. The stability order of the surfaces shown in Table 2 is the same with both the force field and DFT, which confirms that our potential model is suitable for the simulation of orthoenstatite surfaces. However, the relative differences between the surface energy calculated with the force field and DFT (including the BSSE correction) is high. Past work on a similar system, forsterite,^22,25 shows differences in the order of 20 to 40%, though in this case the comparison is with a polarisable force field. If we calculate the values of forsterite surface energy with a rigid ion model, e.g. Kerisit et al.,²⁴ relative differences with DFT become >100%, while still maintaining the correct stability order. This is a direct consequence of rigid-ion models inaccuracy to model mechanical and elastic properties. As discussed in the previous section, polarisable models come with high computational cost that would limit the applicability of this model to simulate the orthoenstatite-fluid interface.

The most stable surface is {210}, which also undergoes very minor surface rearrangement during geometry optimisation, as shown in Fig. 1 and in Table 3. A similar comment holds for {010}, {110}, and {120}, which are slightly less stable.

Fig. 1 View along x of {210} (top left), {010} (top right), {110} (bottom left), and {120} (bottom left) surfaces. The b lattice parameter is shown as a measure of the size of the unit cell; the microfacets are highlighted with a dotted line and their orientation is reported. SiO₄ tetrahedra are coloured in brown, oxygen atoms in red and magnesium atoms in green. The difference between the optimised and unoptimised surfaces is minimal and not visually significant, hence, for each {hkl} set only one representation is provided.

Table 3 Relative (%) differences, obtained with DFT, between selected interatomic distances of the optimised surfaces and those of the surfaces obtained from the bulk: Δ_max and Δ_min are the maximum and minimum relative differences, σ is the standard deviation, and is the average of the absolute relative differences

	{210}	{010}	{110}	{120}	{100}1	{100}2	{100}3
The Si atoms of the SiO4 tetrahedra on the top surfaces have been chosen as references. Si–O is the bond distance between each Si atoms and the O atoms of its tetrahedron (a measure of the rigidity of the SiO4 units during the optimization); Si⋯Si is the distance between each Si atom with its two closest Si atoms (a measure of the relative displacement of the SiO4 centres on the surface); Si⋯Mg is the distance between each Si atom with its four closest Mg atoms (a measure of the relative displacement of Mg atoms); Si⋯O is the distance between each Si atom and its five closest non-bonded O atoms (a measure of the rotation of SiO4 units on the surface). Total is the statistics run on the whole set of the aforementioned distances. A similar atomic rearrangement is observed with the force field.
Si–O
Δ max	3.4	3.7	5.0	4.8	2.2	3.2	1.8
Δ min	−1.3	−1.0	−2.8	−2.6	−2.0	−2.4	−2.0
σ	1.2	1.4	1.6	1.6	1.6	1.5	1.4
	0.9	1.0	1.1	1.2	1.3	1.2	1.1
Si⋯Si
Δ max	2.1	2.2	2.7	2.5	2.6	2.2	1.0
Δ min	0.2	−0.6	−1.1	−1.0	−2.9	−3.2	−2.3
σ	0.7	1.1	1.2	1.1	3.9	2.1	2.4
	1.1	1.4	1.3	1.4	2.8	1.6	1.7
Si⋯Mg
Δ max	5.6	3.4	3.4	3.6	2.6	3.5	7.0
Δ min	−7.7	−5.0	−9.6	−9.1	−17.0	−20.5	−12.3
σ	3.2	2.3	3.0	2.7	7.8	7.0	5.6
	2.5	2.0	2.4	2.3	6.0	4.9	3.7
Si⋯O
Δ max	7.8	6.9	7.4	8.1	11.7	6.3	16.2
Δ min	−3.2	−1.3	−3.5	−4.1	−14.7	−4.2	−12.3
σ	2.7	2.0	2.5	2.3	6.7	2.7	5.2
	2.5	1.8	2.2	1.8	4.5	2.2	3.3
Total
Δ max	7.8	6.9	7.4	8.1	11.7	6.3	16.2
Δ min	−7.7	−5.0	−9.6	−9.1	−17.0	−20.5	−12.3
σ	2.6	2.0	2.6	2.4	6.2	4.4	5.2
	1.8	1.6	1.9	1.8	3.9	2.6	3.3

---

Of the three main surfaces, and of all the surfaces that have been observed experimentally, {100} is the least stable. A higher surface energy is also associated with all surfaces that have {100} as a main microfacet (i.e. {410} and {310}, see Fig. 1 in the ESI†). The {120} surface is not affected because the {100} facet is very small compared to the {010} facet, though the main structural rearrangement observed during geometry optimisation is in proximity of the {100} facet.

Fig. 2 shows that the {100} surface terminates with Mg ions with low coordination. In order to achieve dipole neutrality without changing the SiO₄ bonding, there are four distinct ways of arranging 2 Mg atoms on the surface. One corresponds to a higher surface energy with both force fields and DFT, leading to some major surface reconstruction. Three undergo some reconstruction and provide reasonable structures and energies, though they have a different stability order with FF and DFT. Therefore, we have reported the data for all these arrangements. In all these cases, the {100} surface undergoes surface reconstruction due to rebuilding the coordination environment of the exposed Mg atoms, resulting in rotations of SiO₄ tetrahedra and Mg displacement (as also seen in Table 3).

Fig. 2 View along x of {100}₁ (top), {100}₂ (middle) and {100}₃ (bottom) as cut from the bulk (left) and optimised with DFT (right). Force field optimization leads to similar structures. Colours and labels as in Fig. 1.

With the adopted basis set, a BSSE of about 25% on the surface energy is estimated for forsterite surfaces²⁵ and for the surfaces of orthoenstatite that do not undergo any main reconstruction. Such contribution increases to 30–40% for the 100 surface. Here, there is a limitation in the way the BSSE is calculated, which is that its accuracy decreases if the surface undergoes major reconstruction (as it is estimated from the unoptimised surface cut). For this reason, the value of the surface energy in the hypothesis of a 25% BSSE correction is reported into parenthesis for the {100} surfaces in Table 2.

The fact that {100} has often been observed experimentally and not just as “one of” the surfaces, but rather as “the” main surface exhibited by orthoenstatite, is in contrast with its high surface energy and unlikely atomic arrangement. Here, a reasonable hypothesis is that when in contact with fluids or other materials, the upper Mg atoms may leave the surface and the surface may undergo hydroxylation (an hypothesis that would be in line with what was observed in Oelkers and Schott⁶⁶ and Zakaznova-Herzog et al.⁶⁷), or other forms of reconstruction altering the Si–O bonding (e.g. the formation of a phyllosilicate layer, as observed in some experiments^68–71).

In order to test these hypotheses, two models have been designed. In one, magnesium atoms have been completely removed from the top and bottom of the {100} slab (2 per side) and replaced with four hydrogen atoms per side bound to the most external non-bridging oxygen atoms (i.e. bound to one silicon atom only, see Fig. 3). In the other, magnesium atoms and the two apical oxygen atoms in the unit cell have been removed. In both cases, constant volume geometry optimization was run using DFT. In the first case, the hydroxyl groups form a network of hydrogen bonds that connect the chains of tetrahedra on the surface. In the second case, a distorted phyllosilicate layer forms from joining the SiO₄ chains through rotation and rearrangement of the tetrahedra to compensate for the oxygen vacancies. Both these arrangements are represented in Fig. 3 and do not affect the structure of the slab beyond the superficial layer.

Fig. 3 View along z of the most external silicate layer of the optimised {100} surface: stoichiometric (left), hydroxylated (middle, H₄-S₁₀₀-H₄ in Table 4), and with removal of apical oxygen (right, S₁₀₀ in Table 4). Hydrogen bonds are shown in green dashed lines, H atoms are shown in white colour. The other colours are as in Fig. 1.

While it is important to highlight that these are just two preliminary models and that further investigation is required to achieve a conclusive understanding of the {100} surface structure and stability, the results of the optimisations (run in vacuum) are in good agreement with what is observed experimentally during weathering. Through titration and electrokinetic measurements, Oelkers et al.⁷² observe loss of divalent cations coupled with protonation as a common feature for silicate surfaces, which is also in line with the outcomes of the XPS experiments by Zakaznova-Herzog et al.⁶⁷ and Schott et al.⁷³ Numerous other studies based on XPS, AES, TEM, HRTEM, XRD SEM-EDS, SAED investigations show that the chemical components of the {100} surface of orthopyroxenes are different from their inner compositions,⁶⁸ and that polymerisation of silicate chains on the surface⁶⁹ leads to the formation of a phyllosilicate layer.^70,71 Most of these studies highlight that the surface structure and chemistry of weathering depend on external factors such as temperature, pH and chemical composition of the solution.

The results of the simulations discussed in the previous paragraphs show that the way the {100} surface re-arranges in solution is also likely to occur when orthoenstatite is in vacuum or in contact with gases. This is particularly relevant in the context of organic matter in space, as orthoenstatite is often found in meteorites and is present in extra-terrestrial setups where the production of hydrogen and methane takes place. As these surfaces are non-stoichiometric (i.e. the formula unit no longer corresponds to that of the solid due to loss or addition of atoms), there is no obvious way of computing their surface energies and of comparing their relative stabilities. Table 4 shows that a number of reactions leading to the formation of both the non-stoichiometric surfaces are highly favoured, both in gas phase and in water. These data should be considered semi-quantitative: the zero point energy and the thermal contributions have not been calculated due to the size of the system; no BSSE correction has been included; the non-stoichiometric surface models might correspond to local minima (we did not test all the possible hypothesis); and DFT has known limitations in reproducing relative stabilities of solids that are too different in their structure and density.^74,75 Additionally, the estimation of the reactions in water have been obtained through including the experimental hydration free energies for Mg²⁺ (−1830 kJ mol⁻¹ (ref. 76)), H₃O⁺ and H₂O (−461.30 kJ mol⁻¹ and −26.46 kJ mol⁻¹ (ref. 77)), while the hydration of the two orthoenstatite surfaces is neglected as it is expected to cancel nearly exactly. However, while these sources of inaccuracy would affect the absolute number of the reaction energy, it is unlikely that they would change the sign.

Table 4 Energy [kJ mol⁻¹] associated to possible reaction pathways leading to the {100} surface reconstruction via hydroxylation (top) and formation of a phyllosilicate-like layer (bottom) in gas phase and in water

Reaction	ΔUel
Mg2-S100-Mg2 and Mg2O2-S100-O2Mg2 represent the unaltered {100} slabs; H4-S100-H4 is the slab where 2Mg^2+ atoms per side have been replaced with 4H+ atoms; S100 is the slab where 2MgO per side have been removed. Energies do not include thermal contributions, the zero point energy, the hydration free energy of the solids, and the BSSE correction.
Mg2-S100-Mg2 + 4H2O(g) → H4-S100-H4 + 4MgO	−1035
Mg2-S100-Mg2 + 8H2O(g) → H4-S100-H4 + 4Mg(OH)2	−1437
Mg2-S100-Mg2 + 4Mg(OH)2 → H4-S100-H4 + 8MgO	−633
Mg2-S100-Mg2 + 8H2O(aq) → H4-S100-H4 + 4Mg(OH)2	−1226
Mg2-S100-Mg2 + 8H3O(aq)^+ → H4-S100-H4 + 4Mg(aq)^2+ + 8H2O(aq)	−475
Mg2O2-S100-O2Mg2 → S100 + 4MgO	−201
Mg2O2-S100-O2Mg2 + 4H2O(g) → S100 + 4Mg(OH)2	−603
Mg2O2-S100-O2Mg2 + 4H2O(aq) → S100 + 4Mg(OH)2	−497
Mg2O2-S100-O2Mg2 + 4H3O(aq)+ → S100 + 2Mg(OH)2 + 4H2O(aq) + 2Mg(aq)^2+	−122

---

The surfaces in the second part of Table 2 are all shown in Fig. 1 of the ESI.† Given that they are all microfacets of the three main surfaces, we can assume that what was discussed above for {210}, {010} and {100} will apply to this set of surfaces. In particular, most surfaces do not undergo any major reconstruction with the force field, except for {410} and {310} due to {100} being their main facet. This is why they also require more layers to converge, due to their considerably larger area with respect to that of the {100} unit cell. These facets would most likely be hydroxylated or reconstructed in a simular way as described above for the {100} surface.

4 Conclusions

Modelling is becoming an increasingly powerful tool able to provide accurate predictions and atomic insights into experimental findings. In order to model complex phenomena, such as the adsorption and reactivity of CO₂ at the mineral-fluid interface, detailed knowledge of the mineral surfaces is required. Within this context, this work:

• provides a detailed analysis of the atomic structure and stability of the surfaces of orthoenstatite, which represents the background knowledge necessary to start modelling the interface with gas-phase and dissolved CO₂;

• shows that only three stoichiometric non-polar surfaces exist, the others being microfacets. The restricted number of possible main surfaces and large number of possible microfacets, together with their similar surface energy, could justify why such a variety of faces are observed in natural samples;

• shows that one of the aforementioned surfaces, {100}, is the least stable cut despite being that most frequently observed in natural samples, meaning it will likely exists as a non-stoichiometric surface (i.e. either hydroxylated or with Si–O bonding structure altered);

• provides a force field that models structural parameters and surface stability order with similar accuracy as DFT and that can be used, together with the surface models discussed above, as a starting point to develop a model of the orthoenstatite-CO₂–water interface.

These fundamental details on orthoenstatite surface structure, stability and morphology add to the knowledge we have on other better characterised silicate phases (e.g. olivines) and can be used together in future works to characterise non-stoichiometric surfaces, to identify binding sites for water and CO₂ (also in presence of cationic substitutions and other surface defects), to investigate physical and chemical adsorption mechanisms of these and other molecules, and to shed light on the CO₂ reduction mechanism in both terrestrial and extra-terrestrial environments. The models and tools developed in this work can also be expanded in future works to characterise the fundamental interactions leading to the dissolution/reprecipitation processes involved in CO₂ geosequestration.

Author contributions

J. A. and R. D. developed the force field and wrote the paper. I. M. and R. D. run the DFT calculations.

Conflicts of interest

There are no conflicts to declare.

Data availability

Force field description has been included as ESI.† The information that is needed to make the input files for GULP and CRYSTAL are reported in the text. The full set of input files necessary to reproduce the calculations and a selection of output files are stored on the data storage facilities available through the Pawsey Supercomputing Centre and Curtin University, and can be made available upon request.

The codes used for this work can be found at: 1. Lattice dynamics: https://gulp.curtin.edu.au/index.html. 2. DFT, CRYSTAL: https://www.crystal.unito.it/. 3. Visualisation: https://jmol.sourceforge.net/ and https://github.com/arohl/gdis.

Acknowledgements

RD thanks Curtin University and the Australian Research Council (FT180100385) for supporting this research. The Pawsey Supercomputing Centre and the Australian National Computational Infrastructure are also acknowledged for the provision of computing time through the NCMAS and Pawsey Partners merit allocation schemes. Curtin Library and Annalen der Physik teams are also acknowledged for facilitating access to historical bibliographic references.

References

1. K. O. Yoro and M. O. Daramola, in Chapter 1 - CO₂ emission sources, greenhouse gases, and the global warming effect, ed. M. R. Rahimpour, M. Farsi and M. A. Makarem, Woodhead Publishing, 2020, pp. 3–28.
2. S. I. Seneviratne, M. G. Donat, A. J. Pitman, R. Knutti and R. L. Wilby, Nature, 2016, 529, 477–483.
3. D. Berstad, R. Anantharaman and P. Nekså, Int. J. Refrig., 2013, 36, 1403–1416.
4. A. Goli, A. Shamiri, A. Talaiekhozani, N. Eshtiaghi, N. Aghamohammadi and M. K. Aroua, J. Environ. Manage., 2016, 183, 41–58.
5. I. Sreedhar, R. Vaidhiswaran, B. M. Kamani and A. Venugopal, Renewable Sustainable Energy Rev., 2017, 68, 659–684.
6. W. Gao, T. Zhou, Y. Gao, B. Louis, D. O'Hare and Q. Wang, J. Energy Chem., 2017, 26, 830–838.
7. X. Zhang, S.-X. Guo, K. A. Gandionco, A. M. Bond and J. Zhang, Adv. Mater., 2020, 7, 100074.
8. T. K. Todorova, M. W. Schreiber and M. Fontecave, ACS Catal., 2020, 10, 1754–1768.
9. O. O. Müntener, Geology, 2010, 38, 959–960.
10. Climeworks, Orca: the first large-scale plant, https://climeworks.com/plant-orca – accessed on 14 Jan 2025.
11. P. Kelemen, S. M. Benson, H. Pilorgé, P. Psarras and J. Wilcox, Front. Clim., 2019, 1,  DOI:10.3389/fclim.2019.00009.
12. P. B. Kelemen and J. Matter, Proc. Natl. Acad. Sci. U. S. A., 2008, 105, 17295–17300.
13. I. M. Power, S. A. Wilson and G. M. Dipple, Elements, 2013, 9, 115–121.
14. O. Plümper and J. Matter, Elements, 2023, 19, 165–172.
15. A. Neubeck, N. T. Duc, D. Bastviken, P. Crill and N. G. Holm, Geochem. Trans., 2011, 12, 6.
16. T. M. McCollom and J. S. Seewald, Elements, 2013, 9, 129–134.
17. T. M. McCollom, F. Klein and M. Ramba, Icarus, 2022, 372, 114754.
18. T. M. McCollom, Rev. Mineral. Geochem., 2013, 75, 467–494.
19. T. Tomkinson, M. R. Lee, D. F. Mark and C. L. Smith, Nat. Commun., 2013, 4, 2662.
20. J. Olsson, N. Bovet, E. Makovicky, K. Bechgaard, Z. Balogh and S. Stipp, Geochim. Cosmochim. Acta, 2012, 77, 86–97.
21. T. Liu, S. Gautam, H. Wang, L. M. Anovitz, E. Mamontov, L. F. Allardd and D. R. Cole, Phys. Chem. Chem. Phys., 2018, 20, 27822–27829.
22. N. H. de Leeuw, S. C. Parker, C. R. A. Catlow and G. D. Price, Phys. Chem. Miner., 2000, 27, 332–341.
23. G. Watson, P. M. Oliver and S. C. Parker, Phys. Chem. Miner., 1997, 25, 70–78.
24. S. Kerisit, J. H. Weare and A. R. Felmy, Geochim. Cosmochim. Acta, 2012, 84, 137–151.
25. R. Demichelis, M. Bruno, F. R. Massaro, M. Prencipe, M. De La Pierre and F. Nestola, J. Comput. Chem., 2015, 36, 1439–1445.
26. J. Navarro-Ruiz, M. Sodupe, P. Ugliengo and A. Rimola, Phys. Chem. Chem. Phys., 2014, 16, 17447–17457.
27. M. Cameron and J. J. Papike, Am. Mineral., 1981, 66, 1–50.
28. S. Ghose, V. Schomaker and R. K. McMullan, Z. Kristallogr. – New Cryst. Struct., 1986, 176, 159–175.
29. R. Demichelis, H. Suto, Y. Noël, H. Sogawa, T. Naoi, C. Koike, H. Chihara, N. Shimobayashi, M. Ferrabone and R. Dovesi, Mon. Not. R. Astron. Soc., 2012, 420, 147–154.
30. R. J. Angel and D. A. Hugh-Jones, J. Geophys. Res.: Solid Earth, 1994, 99, 19777–19783.
31. R. Zucker and S. Shim, Am. Mineral., 2009, 94, 1638–1646.
32. C. Stangarone, M. Tribaudino, M. Prencipe and P. P. Lottici, J. Raman Spectrosc., 2016, 47, 1247–1258.
33. A. Bouissonnie, D. Daval, F. Guyot and P. Ackerer, J. Phys. Chem. C, 2020, 124, 3122–3140.
34. J. D. Gale and A. L. Rohl, Mol. Simul., 2003, 29, 291–341.
35. J. D. Gale, Z. Kristallogr., 2005, 220, 552–554.
36. S. Fleming and A. L. Rohl, Z. Kristallogr. – Cryst. Mater., 2005, 220, 580–584.
37. R. T. Cygan, J. J. Liang and A. G. Kalinichev, J. Phys. Chem. B, 2004, 108, 1255–1266.
38. H. J. C. Berendsen, J. R. Griger and T. P. Straatsma, J. Phys. Chem., 1987, 91, 6269.
39. J. G. Harris and K. H. Yung, J. Phys. Chem., 1995, 99, 12021–12024.
40. B. I. Armstrong, A. Silvestri, R. Demichelis, P. Raiteri and J. D. Gale, Philos. Trans. R. Soc., A, 2023, 381, 202220250.
41. Y. Ohashi, Phys. Chem. Miner., 1984, 10, 217–229.
42. J. M. Jackson, S. V. Sinogeikin and J. D. Bass, Phys. Earth Planet. Inter., 2007, 161, 1–12.
43. R. Dovesi, F. Pascale, B. Civalleri, K. Doll, N. M. Harrison, I. Bush, P. D'Arco, Y. Nöel, M. R′erat, P. Carbonni'ere, M. Caus'a, S. Salustro, V. Lacivita, B. Kirtman, A. M. Ferrari, F. S. Gentile, J. Baima, M. Ferrero, R. Demichelis and M. De La Pierre, J. Chem. Phys., 2020, 152, 204111.
44. A. Erba, J. K. Desmarais, S. Casassa, B. Civalleri, L. Donà, I. J. Bush, B. Searle, L. Maschio, L. Edith-Daga, A. Cossard, C. Ribaldone, E. Ascrizzi, N. L. Marana, J.-P. Flament and B. Kirtman, J. Chem. Theory Comput., 2023, 19, 6891–6932.
45. C. Adamo and V. Barone, J. Chem. Phys., 1999, 110, 6158–6170.
46. F. Pascale, C. M. Zicovich-Wilson, R. Orlando, C. Roetti, P. Ugliengo and R. Dovesi, J. Phys. Chem. B, 2005, 109, 6146–6152.
47. M. Bruno, F. R. Massaro, M. Prencipe, R. Demichelis, M. De La Pierre and F. Nestola, J. Phys. Chem. C, 2014, 118, 2498–2506.
48. R. Demichelis, B. Civalleri, M. Ferrabone and R. Dovesi, Int. J. Quantum Chem., 2010, 110, 406–415.
49. H. J. Monkhorst and J. D. Pack, Phys. Rev. B, 1976, 13, 5188–5192.
50. P. Pulay, Chem. Phys. Lett., 1980, 73, 393–398.
51. P. Pulay, J. Comput. Chem., 1982, 3, 556–560.
52. C. G. Broyden, J. Inst. Math. Appl., 1970, 6, 76–90.
53. R. Fletcher, Comput. J., 1970, 13, 317–322.
54. D. Goldfarb, Math. Comput., 1970, 24, 23–26.
55. D. F. Shanno, Math. Comput., 1970, 24, 647–656.
56. S. F. Boys and F. Bernardi, Mol. Phys., 1970, 19, 553–566.
57. M. De La Pierre, M. Bruno, C. Manfredotti, F. Nestola, M. Prencipe and C. Manfredotti, Mol. Phys., 2014, 112, 1030–1039.
58. G. vom Rath, Ann. Phys., 1869, 214, 515–550.
59. W. C. Brögger and G. vom Rath, London, Edinburgh Dublin Philos. Mag. J. Sci., 1876, 2, 379–387.
60. V. Goldschmidt, Atlas der Kristallformen, Universitätsverlag Winters Heidelberg, 1913–1923, vol. 7.
61. W. E. Troger, H. U. Bambauer, F. Taborszky and H. D. Trochim, Optical determination of rock-forming minerals. part 1, Determinative tables, Schweizerbartsche Verlagsbuchhandlung, Stuttgart, West Germany, 1979, pp. 72–76.
62. W. A. Deer, R. A. Howie and J. Zussman, An introduction to the rock-forming minerals, Longman Scientific & Technical, Harlow, England, 2nd edn, 1992, p. 155.
63. J. W. Anthony, R. A. Bideaux, K. W. Bladh and M. C. Nichols, Handbook of Mineralogy, Mineralogical Society of America, Chantilly, VA 20151-1110, USA, 2001, https://www.handbookofmineralogy.org/pdfs/enstatite.pdf.
64. F. Zambonini, Mineralogia Vesuviana, Reale Accademia delle Scienze Fisiche e Matematiche di Napoli, 1910, pp. 140–143.
65. J. Ralph, D. Von Bargen, P. Martynov, J. Zhang, X. Que, A. Prabhu, S. M. Morrison, W. Li, W. Chen and X. Ma, Am. Mineral., 2024, 110, 833–844.
66. E. H. Oelkers and J. Schott, Geochim. Cosmochim. Acta, 2001, 65, 1219–1231.
67. V. P. Zakaznova-Herzog, H. W. Nesbitt, G. M. Bancroft and J. S. Tse, Geochim. Cosmochim. Acta, 2008, 72, 69–86.
68. G. L. J. Xue, Y. Pan and H. Lu, Chin. Sci. Bull., 2003, 48, 931–934.
69. Y. Nakajima and P. H. Ribbe, Contrib. Mineral. Petrol., 1981, 78, 230–239.
70. R. A. Eggleton, Clays Clay Miner., 1982, 30, 11–20.
71. I. Ohnishi and K. Tomeoka, Meteorit. Planet. Sci., 2007, 42, 49–61.
72. E. H. Oelkers, S. V. Golubev, C. Chairat, O. S. Pokrovsky and J. Schott, Geochim. Cosmochim. Acta, 2009, 73, 4617–4634.
73. J. Schott, R. A. Berner and E. L. Sjöberg, Geochim. Cosmochim. Acta, 1981, 45, 2123–2135.
74. R. Demichelis, P. Raiteri, J. D. Gale and R. Dovesi, J. Phys. Chem. C, 2013, 117, 17814–17823.
75. R. Demichelis, B. Civalleri, P. D'Arco and R. Dovesi, Int. J. Quantum Chem., 2010, 110, 2260–2273.
76. Y. Marcus, Biophys. Chem., 1994, 51, 111–127.
77. J. R. Pliego Jr and J. M. Riveros, Phys. Chem. Chem. Phys., 2002, 4, 1622–1627.

---

Footnote

† Electronic supplementary information (ESI) available. See DOI: https://doi.org/10.1039/d5nr01501d

---