Structure, thermal expansion and incompressibility of MgSO₄·9H₂O, its relationship to meridianiite (MgSO₄·11H₂O) and possible natural occurrences

Since being discovered initially in mixed-cation systems, a method of forming end-member MgSO₄·9H₂O has been found. We have obtained powder diffraction data from protonated analogues (using X-rays) and deuterated analogues (using neutrons) of this compound over a range of temperatures and pressures. From these data we have determined the crystal structure, including all hydrogen positions, the thermal expansion over the range 9–260 K at ambient pressure, the incompressibility over the range 0–1.1 GPa at 240 K and studied the transitions to other stable and metastable phases. MgSO₄·9D₂O is monoclinic, space group P2₁/c, Z = 4, with unit-cell parameters at 9 K, a = 6.72764 (6), b = 11.91154 (9), c = 14.6424 (1) Å, β = 95.2046 (7)° and V = 1168.55 (1) ų. The structure consists of two symmetry-inequivalent Mg(D₂O)₆ octahedra on sites of symmetry. These are directly joined by a water–water hydrogen bond to form chains of octahedra parallel with the b axis at a = 0. Three interstitial water molecules bridge the Mg(D₂O)₆ octahedra to the SO₄²⁻ tetrahedral oxyanion. These tetrahedra sit at a ≃ 0.5 and are linked by two of the three interstitial water molecules in a pentagonal motif to form ribbons parallel with b. The temperature dependences of the lattice parameters from 9 to 260 K have been fitted with a modified Einstein oscillator model, which was used to obtain the coefficients of the thermal expansion tensor. The volume thermal expansion coefficient, α_V, is substantially larger than that of either MgSO₄·7D₂O (epsomite) or MgSO₄·11D₂O (meridianiite), being ∼ 110 × 10⁻⁶ K⁻¹ at 240 K. Fitting to a Murnaghan integrated linear equation of state gave a zero-pressure bulk modulus for MgSO₄·9D₂O at 240 K, K₀ = 19.5 (3) GPa, with the first pressure derivative of the bulk modulus, K′ = 3.8 (4). The bulk modulus is virtually identical to meridianiite and only ∼ 14% smaller than that of epsomite. Above ∼ 1 GPa at 240 K the bulk modulus begins to decrease with pressure; this elastic softening may indicate a phase transition at a pressure above ∼ 2 GPa. Synthesis of MgSO₄·9H₂O from cation-pure aqueous solutions requires quench-freezing of small droplets, a situation that may be relevant to spraying of MgSO₄-rich cryomagmas into the surface environments of icy satellites in the outer solar system. However, serendipitously, we obtained a mixture of MgSO₄·9H₂O, mirabilite (Na₂SO₄·10H₂O) and ice by simply leaving a bottle of mid-winter brine from Spotted Lake (Mg/Na ratio = 3), British Columbia, in a domestic freezer for a few hours. This suggests that MgSO₄·9H₂O can occur naturally – albeit on a transient basis – in certain terrestrial and extraterrestrial environments.

1. Introduction

1.1. Background

Crystalline hydrates of the general form MgSO₄·nH₂O with n = 1, 1.25, 2, 2.5, 3, 4, 5, 6, 7 and 11 have been known for over a century and their ranges of stability and metastability are reasonably well established (D'Ans, 1933; Hodenberg & Kühn, 1967). The lack of any analogues with n = 11 in other divalent metal sulfates, and the curious gap in hydration states between n = 7 and n = 11 led us to investigate methods of synthesizing hitherto unknown compounds and new hydration states. Our attention was originally drawn to this gap between n = 7 and n = 11 by high-pressure studies of MgSO₄·11H₂O where we observed an apparent decomposition at ∼ 0.9 GPa, 240 K to a mixture of ice VI and what was – at the time – an unknown hydrate (Fortes et al., 2017). The unidentified phase did not correspond to any of the high-pressure phases of MgSO₄·7D₂O we had encountered in other high-pressure experiments (Gromnitskaya et al., 2013) and one possible explanation was the occurrence of an intermediate hydration state that became stable only under pressure.

In our experiments we initially adopted chemical substitution as a crude proxy for hydrostatic stress, the expectation being that this would create analogue structures of the phase observed at high pressure. The work paid dividends in terms of characterizing the extent to which various divalent metal cations can substitute for Mg²⁺ in the meridianiite structure and their effect on the unit-cell parameters (Fortes, Browning & Wood, 2012a,b), and also revealed crystals isotypic with meridianiite formed from MgCrO₄ and MgSeO₄ (Fortes & Wood, 2012: Fortes, 2015). Furthermore, by substituting Mg²⁺ with Zn²⁺, Ni²⁺ or Fe²⁺, we determined that novel hydration states could indeed be formed in mixed-cation compounds [i.e. (Mg_1 − xNi_x)SO₄ for 0.3 < x < 0.8] by rapid quenching of aqueous solutions in liquid nitrogen. These included a monoclinic enneahydrate for all of the aforementioned cations and a triclinic octahydrate for Ni²⁺. The discovery of MgSO₄·9H₂O also resolved a long-standing question over the origin of weak parasitic peaks observed in our first neutron powder diffraction measurement on MgSO₄·11D₂O, which did not appear to match any other known hydrate or plausible contaminant (Fortes et al., 2008, and Fig. S1 of the supporting information).

Following its identification, the new enneahydrate became our best candidate for the phase that forms when meridianiite decomposes at high pressure, implying that MgSO₄·9H₂O might occur naturally at depth inside the largest icy satellites of Jupiter, where MgSO₄/Na₂SO₄-brine oceans are believed to be present (Khurana et al., 1998; Kivelson et al., 2002; Vance & Brown, 2013; Vance et al., 2014) and where MgSO₄ and Na₂SO₄ hydrates have been reported in surface deposits from near-IR spectroscopy (Orlando et al., 2005; Dalton et al., 2005; Shirley et al., 2010). Since it may be a `rock-forming' mineral in the outer solar system, and since we have also found it in slowly frozen natural brine (see §4), this compound and its relationship to other cryohydrates acquires mineralogical significance.

For brevity, the composition of these phases will be written in the form MSn, where M = Mg, S = SO₄ and n identifies the number of water molecules per formula unit. Hence, MS9 would mean MgSO₄·9H₂O, MS11 = MgSO₄·11H₂O and so on. The addition of `(D)' indicates a deuterated analogue, such that MS9(D) means MgSO₄·9D₂O.

1.2. Experimental objectives

The objectives of this work are: (i) to determine the complete structure of MS9, including all hydrogen positions for comparison with the recently determined structure of MgSeO₄·9H₂O (Fortes et al., 2015); (ii) to determine the unit-cell parameters of MS9 as a function of temperature and pressure in order to characterize the thermal expansion and compressibility of the structure; (iii) to obtain data suitable for testing our hypothesis that MS11 decomposes to MS9 + ice at high pressure. We have used a combination of techniques to achieve these goals, including X-ray and neutron powder diffraction and density functional theory (DFT) calculations. A necessary consequence of using neutron scattering is the requirement to employ a deuterated analogue, MS9(D), in order to eliminate the undesirable incoherent neutron scattering from ¹H. Although there are precedents where materials exhibit substantially different phase behaviour on deuteration, our experience of working with many different forms of ice and salt hydrates is that deuteration has a negligible effect on the molar volume and shifts by just a few percent the values of thermoelastic parameters and locations of phase boundaries (e.g. Pistorius, 1968; Fortes, Wood, Tucker & Marshall, 2012).

This work forms the middle part of a series of three papers. The first deals with the high-pressure behaviour of meridianiite (MS11) and its decomposition to the enneahydrate: the third paper, which will be presented elsewhere, deals with the Ni-bearing species, NiSO₄·9H₂O and its solid solution series with MS9, and with the structure of NiSO₄·8H₂O.

2. Experimental method

2.1. Sample synthesis

Aqueous solutions were prepared from MgSO₄ (Sigma-Aldrich M7506, anhydrous β-phase, ReagentPlus^®, ≥ 99.5%) in either H₂O (Alfa-Aesar, ACS Reagent Grade, 36645) or D₂O (Aldrich 151882 99.9 atom % D). We endeavoured to prepare the most concentrated solutions possible, since this reduces the amount of accessory ice in the final sample. However, highly supersaturated solutions of either MgSO₄ tend to be quite labile, crystallizing very quickly, so the typical concentrations used here are in the region of 30 wt % dissolved solids, yielding mixtures containing ∼ 70 wt % MS9 and 30 wt % water ice. Once it became clear how deleterious the presence of ice in the sample can be, efforts were made to synthesize and flash-freeze a starting solution with the same stoichiometry as the 9-hydrate itself, i.e. 42.6 wt % MgSO₄. Such a solution may be prepared by very slow evaporation, without agitation, of an initially undersaturated solution (e.g. 25 wt %) at a temperature of 373–383 K. Ordinarily, this process will lead to the crystallization of metastable MS3, but rapid quenching of the solution prior to the appearance of trihydrate crystals can produce a virtually phase pure specimen of the enneahydrate.

Droplets were pipetted directly into liquid nitrogen held in a steel cryomortar, producing glassy solid spherules in the region of 2–4 mm in diameter. The amorphous spherules are pulverized with a nitrogen-chilled steel pestle and then transferred to an appropriate sample holder. This work was carried in the UCL Earth Sciences Cold Rooms at an air temperature of around 261 K so as to minimize the incorporation of ice formed from atmospheric water vapour into the sample.

2.2. X-ray powder diffraction

X-ray powder diffraction data were obtained in Bragg–Brentano parafocussing reflection geometry using a PAN­analytical X'Pert Pro multipurpose powder diffractometer with germanium-monochromated Co Kα₁ radiation (λ = 1.788996 Å, and an X'Celerator multi-strip detector). Data were collected with variable divergence and receiving slits, converted to fixed-slit geometry with the proprietary X'Pert Pro HighScore Plus software package, and exported in an appropriate format for analysis in the GSAS/Expgui package (Larsen & Von Dreele, 2000; Toby, 2001).

Stable low-temperature measurements were achieved using a thermoelectrically cooled cold stage (Wood et al., 2012). This portable cold stage was held in a plastic box filled with dry-ice pellets whilst the powder specimen was prepared and loaded. Powder samples were transferred to the cold stage and packed down using a liquid nitrogen-cooled spatula to form a top-loaded pressed-powder specimen. A cover with an X-ray transparent window was screwed into place, which preserves a cold and dry atmosphere over the specimen during mounting on the diffractometer and subsequent measurement. The presence of a substantial open space (∼ 52 cm³) over the sample and beneath the cover, which is not actively cooled, means that there is scope for sublimation of ice and transport of vapour over the top of the specimen. In early studies of ice-rich enneahydrate samples, we found that the MS9 would transform irreversibly to the more stable undecahydrate, MS11, on time scales of hours at 250 K, presumably by reaction with water vapour. This motivated our efforts to synthesize ice-free samples, simply in order to allow for measurements over a large angular range with little background noise without the specimen transforming to another phase.

2.3. Neutron powder diffraction (ISIS)

Neutron powder diffraction data were collected using the high-resolution powder diffractometer (HRPD; Ibberson et al., 1992; Ibberson, 2009) at the STFC ISIS spallation neutron source, Rutherford Appleton Laboratory, UK. For ambient-pressure analysis, the powdered samples were transferred directly into a pre-cooled sample holder mounted on a cryostat centre stick. The sample holder consisted of an aluminium holder with a slab-geometry space of dimensions (relative to the incident beam) of width = 18 mm, height = 23 mm and depth = 15 mm. Steel-framed vanadium windows were screwed to the front and back of the holder, the exposed Al and steel components on the front face being masked by Gd foil.

Samples were mounted in the cryostat at 100–150 K, where they were observed to be amorphous. Slow incremental warming to 230 K resulted in crystallization, after which the temperature was reduced to 10 K to collect powder diffraction data suitable for structure refinement.

2.3.1. Cation-pure MgSO₄ samples at ambient pressure

An attempt was made initially to synthesize phase-pure MS9(D), which resulted in a degree of disproportionation and formation of a sample containing substantial MS11(D) and some water ice, the latter exhibiting rather broad Bragg peaks. Despite the sub-optimal quality of the specimen, data were collected at base temperature and in 5 K increments on warming back up to 230 K using the long time window (100–200 ms time-of-flight) on HRPD; in backscattering this provides high-resolution data covering d-spacings from 2.2 to 4 Å. This dataset is referred to as Series 1.

Subsequently, a second sample was prepared with a more water-rich stoichiometry. As before, this was loaded at 150 K and crystallized after warming to 220 K. Inevitably, this sample contained more water ice, but the co-existing material was phase-pure 9-hydrate and all Bragg peaks were sharper than those in Series 1. Powder diffraction data were obtained on warming in 5 K increments from 220 K to 265 K using the 100–200 ms time window. These data (Series 2) recorded the onset of MS9(D) transforming to MS11(D) at 250 K, a process that was complete upon reaching 265 K on the ca 1 h timescale of each measurement.

A third specimen consisting, like the second, of ∼ 70 wt % MS9(D) and 30 wt % deuterated ice Ih, was prepared in the same fashion and cooled to 9 K where diffraction data were acquired in both the short and long time windows, 30–130 ms and 100–200 ms, for ∼ 17 and 13 h, respectively. Data were then collected in the 100–200 ms window on warming to 220 K in 10 K increments. In order to evaluate the transformation of the 9-hydrate to the 11-hydrate, the sample was warmed first from 220 to 250 K and data were collected in 20 min snapshots for ∼ 2 h, and then up to 260 K with 20 min shots being acquired sequentially for ∼ 5.5 h. This dataset comprises Series 3.

2.3.2. Cation-pure MgSO₄ samples at non-ambient pressure

A fourth aliquot of flash-frozen material was loaded into a TiZr gas-pressure vessel embedded in dry-ice snow. Crystallization occurred as the sample was warmed from 220 to 240 K under 10 MPa of He gas to a mixture dominated by MS9(D) with ∼ 25 wt % MS7(D). Data were collected in the 30–130 ms time window in HRPD's 90° detector banks, which have a larger solid angle coverage than the backscattering detectors but a lower resolution (Δd/d ≃ 2 × 10⁻³). Using He gas as the pressure-transmitting medium, diffraction patterns were acquired in roughly 100 MPa increments along the 240 K isotherm up to ∼ 540 MPa with counting times of ∼ 12 h each. The unmasked pressure vessel contributes a number of parasitic features to the diffraction pattern, as well as a substantial Q-dependent background, which was corrected by subtraction of an empty-cell measurement collected at the end of the experiment.

2.4. Ab initio calculations

In order to aid in the interpretation of some of our experimental data, and to predict certain quantities we are unable to measure, we carried out a series of first-principles calculations using DFT and the plane-wave pseudopotential method (Hohenberg & Kohn, 1964; Kohn & Sham, 1965). The calculations were carried out using CASTEP (Payne et al., 1992; Segall et al., 2002; Clark et al., 2005) in conjunction with the analysis tools in the Materials Studio software package (https://accelrys.com). A basis-set cut-off of 1200 eV and a 4 × 2 × 2 -point grid (∼ 0.04 Å⁻¹ reciprocal lattice spacing) were required to achieve convergence of better than 1 × 10⁻² GPa in the stress and better than 1 × 10⁻⁴ eV per atom in total energy. We also tested two different gradient-corrected functionals (GGAs), specifically those referred to as `PBE' (Perdew et al., 1996, 1997) and `Wu–Cohen' or `WC' (Wu & Cohen, 2006). In the past, we have successfully used the PW91 or closely related PBE GGAs for these types of materials (Fortes et al., 2014, 2015). However, there is reason to believe that the Wu–Cohen GGA functional is more accurate than PBE for specific types of calculations (e.g. Tran et al., 2007; Haas et al., 2009, 2011), correcting the under-binding commonly encountered in PW91 and PBE simulations.

Structural relaxations under zero-pressure athermal conditions were carried out, starting from the experimental crystal structure obtained at 9 K for MS9 (see §3.1), using the BFGS method (Pfrommer et al., 1997). The relaxations were considered to have converged when the forces on each atom were less than 5 × 10⁻³ eV Å⁻¹ and each component of the stress tensor was smaller than 0.005 GPa. Details of how other specific properties were calculated are given in subsequent sections; the exception to this is the calculation of the elastic constants, which appears in the supporting information.

3. Results

3.1. Structure determination and completion

From prior work (Fortes, Browning & Wood, 2012a,b) it was known that MS9 is monoclinic, with unit-cell parameters at 250 K a = 6.75 Å, b = 11.95 Å, c = 14.65 Å, β = 95.1°, V = 1177 ų with Z = 4. Systematic absences indicated that the most probable space group was P2₁/c; the space group Pc remained a possibility due to the limited number of well separated k-odd 0k0 reflections observable by eye in the low-angle portion of the pattern and since Pawley refinements in both space groups only weakly favoured P2₁/c. We also observed that in X-ray powder data cation-pure MS9 produces 100 peaks of measurable intensity, whilst the 011 peak is of zero (or negligible) intensity. However, in Ni-doped samples of MS9 the pattern is reversed, 011 is present and 100 has no measureable intensity (Fig. S2). We interpreted this as being due to the preferential occupancy of particular Mg²⁺ sites in the structure by Ni²⁺, which requires there to be at least two symmetry-independent octahedral sites. This is possible in the space group Pc with the octahedrally coordinated cations on general positions; should the space group be P2₁/c, then this obliges the cations to sit on special positions and these sites (2a and 2c) occupy separate 011 planes, thus permitting site ordering.

In an effort to avoid unconscious bias in the structure solution process concerning either the space group or the siting of Mg(H₂O)₆ octahedra on general or special positions, the structure was solved initially in the space group Pc. Since the earlier identification of this phase as a 9-hydrate was based purely on unit-cell volume considerations, we also chose to `test' the stoichiometry by solving the structure without any interstitial water and then evaluating difference Fourier maps phased on any partial solution. X-ray powder data obtained by quenching a highly concentrated solution, which contained largely MS9 and only a very small amount of MS11 (< 5 wt %), were used as the basis for the structure solution, which was achieved using the parallel tempering algorithm in FOX, Version 1.9.7.1 (Favre-Nicolin & Černý, 2002, 2004). For the 9-hydrate, FOX was used to construct ideal MgO₆ octahedra with Mg—O distances of 2.08 Å, and ideal SO₄ tetrahedra with S—O distances of 1.48 Å; these were treated as rigid bodies throughout the solution process. In runs of half a million trials each, the crystal structure was optimized against the powder diffraction data, consistently producing very similar structures with chemically sensible arrangements of the ionic polyhedra even after only a few 10 s of thousands of trials. Difference Fourier maps phased on these structures revealed three peaks per formula unit that correspond to the additional water molecules necessary to form a 9-hydrate (Fig. S3). Ultimately the structure with the lowest overall cost function was exported as a CIF file to form the basis for further analysis.

Examination of the most favourable solution revealed that the cations occupied centres of symmetry; consequently, the true space group must indeed be P2₁/c, as expected, with Mg²⁺ ions occupying the 2a and 2c special positions ( site symmetry) and so the structure was transformed into the higher-symmetry setting.

The trial heavy-atom structure was refined against the X-ray powder dataset using GSAS/Expgui with quite stiff bond distance and bond angle restraints, Mg—O bond lengths being restrained to be 2.08 (4) Å, S—O bond lengths being restrained to 1.480 (5) Å, and the SO₄²⁻ oxyanion being forced to adopt ideal tetrahedral symmetry (O—S—O angles = 109.5°). Study of the first coordination shell of the water O atoms at the end of this refinement revealed a sensible pattern of O⋯O vectors of length 2.7–3.0 Å, consistent with hydrogen-bonded contacts, and so pairs of H atoms were sited at distances of 0.98 Å from each O atom along these vectors.

The complete structure was then refined against the 9 K neutron powder diffraction data for MS9(D) (Series 3). Since these data, particularly at short d-spacings, are dominated by Bragg peaks from water ice, it was necessary to apply similar bond distance and angle restraints to those mentioned earlier and to constrain isotropic displacement parameters for light/heavy atoms to shift together. In the later stages of the refinements, restraint weightings were removed.

A graphical representation of the fit is given in Fig. 1, structural parameters are given in the supporting CIF data and selected bond distances and angles are reported in Tables 1 and 2 for comparison with the DFT-derived values.

Table 1
Unit-cell dimensions, bond lengths, angles and other geometric properties of the ionic polyhedra in MS9 compared with the zero-pressure DFT geometry optimizations

The octahedral distortion index, quadratic elongation () and bond-angle variance () are defined by Robinson et al. (1971).

	MgSO4·9D2O (9 K exp)	MgSO4·9H2O (WC GGA)	MgSO4·9H2O (PBE GGA)
a (Å)	6.7273 (2)	6.63048	6.78924
b (Å)	11.9040 (2)	11.94593	12.19521
c (Å)	14.6405 (3)	14.63552	14.78586
β (°)	95.199 (1)	94.0222	94.3563
V (Å^3)	1167.61 (2)	1156.3832	1220.6764
ΔV (%)	–	−0.96	+4.54
S—O1 (Å)	1.479 (4)	1.4911	1.4944
S—O2† (Å)	1.475 (4)	1.4719	1.4763
S—O3 (Å)	1.491 (4)	1.4948	1.4981
S—O4 (Å)	1.479 (4)	1.5063	1.5078
O1—S—O2 (°)	109.8 (3)	110.32	110.12
O1—S—O3 (°)	108.7 (3)	109.38	109.35
O1—S—O4 (°)	110.0 (3)	108.93	109.09
O2—S—O3 (°)	109.6 (3)	110.20	110.15
O2—S—O4 (°)	110.0 (3)	109.35	109.43
O3—S—O4 (°)	108.8 (3)	108.64	108.68
SO4 volume (Å^3)	1.667	1.7009	1.7117
ΔV (%)	–	+2.0	+2.7
Mg1—Ow1 (Å)	2.028 (5)	2.0546	2.0751
Mg1—Ow2 (Å)	2.075 (5)	2.0946	2.1150
Mg1—Ow3 (Å)	2.056 (5)	2.1320	2.1417
Ow1—Mg1—Ow2 (°)	90.3 (2)	90.24	90.26
Ow1—Mg1—Ow3 (°)	92.3 (2)	90.69	90.62
Ow2—Mg1—Ow3 (°)	91.4 (2)	93.80	93.35
Mg1 volume (Å^3)	11.521	12.2054	12.5101
ΔV (%)	–	+5.9	+8.6
Dist. index	0.0080	0.0125	0.0112
λoct	1.0009	1.0020	1.0015
σ^2oct (°^2)	2.6586	5.4403	4.2356
Mg2—Ow4 (Å)	2.034 (4)	2.0696	2.0877
Mg2—Ow5‡ (Å)	2.108 (5)	2.1307	2.1448
Mg2—Ow6 (Å)	2.051 (4)	2.0858	2.1019
Ow4—Mg2—Ow5 (°)	90.1 (2)	90.02	90.60
Ow4—Mg2—Ow6 (°)	90.2 (2)	91.28	90.94
Ow5—Mg2—Ow6 (°)	90.6 (2)	90.89	91.08
Mg2 volume (Å^3)	11.725	12.2585	12.5439
ΔV (%)	–	+4.0	+7.0
Dist. index	0.0142	0.0113	0.0105
λoct	1.0005	1.0005	1.0005
σ^2oct (°^2)	0.1396	0.8798	0.8770

†Accepts two hydrogen bonds (all others accept three).
‡Hydrogen-bond acceptor in tetrahedral coordination.

Table 2
Hydrogen-bond lengths and angles in MS9 at 9 K

	D—O—D (°)	O—D (Å)	D⋯O (Å)	O⋯O (Å)	O—D⋯O (°)
Ow1—D1a⋯O2	106.5 (6)	0.968 (4)	1.761 (6)	2.699 (6)	162.3 (5)
Ow1—D1b⋯Ow9		0.970 (4)	1.811 (6)	2.775 (6)	171.9 (5)
Ow2—D2a⋯Ow5	103.8 (6)	0.971 (4)	2.051 (7)	3.014 (7)	171.0 (5)
Ow2—D2b⋯Ow9		0.974 (4)	1.834 (7)	2.786 (7)	164.9 (5)
Ow3—D3a⋯Ow8	102.2 (6)	0.978 (4)	1.880 (7)	2.855 (6)	175.8 (5)
Ow3—D3b⋯O1		0.964 (4)	1.933 (6)	2.894 (6)	174.4 (5)
Ow4—D4a⋯Ow7	106.5 (6)	0.969 (4)	1.784 (6)	2.734 (6)	165.6 (5)
Ow4—D4b⋯O1		0.976 (4)	1.822 (6)	2.797 (6)	176.7 (5)
Ow5—D5a⋯Ow7	107.1 (6)	0.980 (4)	1.695 (7)	2.671 (7)	173.6 (6)
Ow5—D5b⋯O4		0.971 (4)	1.753 (6)	2.723 (6)	176.4 (5)
Ow6—D6a⋯Ow8	106.1 (6)	0.982 (4)	1.771 (7)	2.715 (6)	160.3 (5)
Ow6—D6b⋯O3		0.957 (4)	1.774 (6)	2.722 (6)	170.5 (5)
Ow7—D7a⋯O3	103.8 (6)	0.989 (4)	1.793 (6)	2.859 (6)	164.6 (5)
Ow7—D7b⋯O2		0.988 (4)	1.750 (6)	2.730 (6)	171.7 (5)
Ow8—D8a⋯O4	103.7 (6)	0.992 (4)	1.882 (6)	2.847 (6)	163.5 (5)
Ow8—D8b⋯O1		0.985 (4)	1.942 (6)	2.925 (6)	175.1 (5)
Ow9—D9a⋯O3	102.1 (6)	0.972 (4)	1.859 (6)	2.813 (6)	166.4 (5)
Ow9—D9b⋯O4		0.979 (4)	1.792 (6)	2.766 (6)	172.6 (5)

Figure 1 Rietveld refinement of the MgSO4·9D2O structure against neutron powder diffraction data obtained at 9 K in HRPD's 30–130 ms window and 100–200 ms window. The graphical output for each histogram's refinement has been stitched together at d = 2.4 Å. Filled red circles report the observed data, solid green lines the calculated fit and the purple line underneath is obs − calc. Note that, whilst MS9 represents almost 50% by weight of the sample, its larger unit-cell and lower symmetry mean that the intensity of the strongest MS9 peaks is many times less than water ice; consequently, the vertical scale has been greatly expanded, with the effect that the difference curve – particularly in the vicinity of the ice peaks – appears noisy. Tick marks for water ice are uppermost and those for MgSO4·9D2O are below them.

The structure (Fig. 2) consists of isolated Mg(H₂O)₆ octahedra on sites of symmetry; the Mg1 octahedron, which occupies the 2a position, donates a hydrogen bond (via water molecule Ow2) to its neighbouring octahedron, situated on the 2c site, via water molecule Ow5. This effectively creates hydrogen-bonded chains of Mg(H₂O)₆ octahedra running parallel to the b-axis at a = 0. Sulfate tetrahedra sit at a ≃ 0.5 linked by the interstitial water molecules Ow7 and Ow9 in a pentagonal motif to form ribbons, which also run parallel to the b-axis (Fig. 3). The third interstitial water molecule, Ow8, cross-links adjacent chains of MgO₆ octahedra along the c-axis, accepting hydrogen bonds from Ow3 (Mg1-coordinated) and Ow6 (Mg2-coordinated).

Figure 2 Asymmetric unit of MS9 with some symmetry-related atoms included for clarity. Symmetry codes: (i) x-1, y, z; (ii) -x, -y, -z; (iii) -x, 1-y, -z. Drawn using Diamond (Putz & Brandenburg, 2006) and rendered with PovRay (POV, 2004).

Figure 3 Packing and hydrogen-bond connectivity of the polyhedra in MS9 in the bc plane. Symmetry codes: (i) ; (ii) -x, -y, -z; (iii) ; (iv) ; (v) -x, 1-y, -z. Drawn using Diamond (Putz & Brandenburg, 2006) and rendered with PovRay (POV, 2004).

The sulfate O atoms O1, O3 and O4 each accept three hydrogen bonds, whilst O2 accepts only two. The effect of this is typically manifested in the S—O bond lengths, those with the least number of hydrogen-bond connections tending to have the shortest S—O length. Whilst this is evident from the DFT calculations (Table 1), the differences in S—O bond lengths obtained experimentally show no such systematic variation.

Both Mg(H₂O)₆ octahedra are elongated along the axis of the inter-octahedral hydrogen-bond donor–acceptor (Mg1—Ow2 and Mg2—Ow5). The occurrence of hydrogen-bond donation to cation-coordinated water is fairly common in divalent metal sulfate hydrates, although it does not occur in the non-isotypic crystal MgSeO₄·9H₂O (Fortes et al., 2015). The 2c octahedral sites are larger and more regular (in terms of bond-angle variance) than the 2a sites in both Mg and Ni endmembers.

In the majority of inorganic hydrates the cation is either fully saturated (i.e. its first coordination shell is entirely filled by water O atoms) or is undersaturated, such that some of the cation-coordinated O atoms are O²⁻ ions rather than neutral H₂O; this tends to result in corner- and/or edge-sharing polyhedra. In `cation-saturated' hydrates, the water molecules generally donate hydrogen bonds to the anions (or oxyanions) and less frequently donate hydrogen bonds to other water molecules, thus the role of water is to act as a bridge between the ions rather than associating with any neighbouring neutral water.

In hydrates where there is excess water above that required to saturate the cation, the question then arises as to whether the additional water just serves to extend these bridges or whether there is an opportunity for association with neighbouring water to form some kind of polymeric unit (e.g. clusters, ribbons or sheets). The structures of both MgSeO₄·9H₂O and MS11 contain quite large clusters of hydrogen-bonded water; the selenate 9-hydrate contains a centrosymmetric dodecamer and the sulfate 11-hydrate contains a centrosymmetric hexadecamer. The structure of MgSO₄·9H₂O, on the other hand, has no such extended polymeric water structure. What we see instead are a trimer and a linear hexamer (Fig. S4). The trimers are very common in highly hydrated materials such as these, and identical polymeric units are found in MS11, MgSeO₄·9H₂O (Fortes et al., 2015) and Fe³⁺(ClO₄)₃·9H₂O (Hennings et al., 2014a), for example.

Since we expect the local structure of the aqueous solution from which these compounds are obtained to be rich in five- and six-sided rings, the occurrence of closed ring-like structures in MS11 and MgSeO₄·9H₂O is unsurprising. However, it appears that these are substantially less common than finite branched and unbranched chains. Branched or decorated (H₂O)₉ chains occur in both trivalent and divalent metal bromide and iodide enneahydrates (Hennings et al., 2013; Schmidt et al., 2014), whereas branched (H₂O)₈ chains appear in CaBr₂.9H₂O (Hennings et al., 2014b). Of course, the principal difference with our title compound is the 1:2 or 1:3 cation-to-anion ratio in these other materials. We might anticipate that in the 1:1 cation-to-anion compounds with increasing numbers of water molecules per formula unit that the interstitial water would be free to adopt a more liquid-like arrangement, comprising closed rings. That it does not do so in MS9, forming instead a linear six-membered chain, may be related to its metastability.

3.2. DFT structural relaxations

The structure of MS9 was optimized using both the PBE and WC functionals starting from the 9 K experimental structure. The unit-cell parameters and atomic coordinates output from these simulations are provided in the supporting CIFs. We found that PBE led to an athermal unit-cell volume ∼ 4.5% larger than the experimental value measured at 9 K, whereas the WC functional gave a zero-pressure athermal cell just 1% smaller (Table 1). Both S—O and Mg—O bond lengths are in closer agreement with experiment in the WC GGA calculations, but are still over-estimated. In the case of both functionals, the correlation between observed and calculated O—H⋯O angles and H⋯O distances is very good, which is to say that the overall intermolecular connectivity is well reproduced. PBE and WC yield closely similar hydrogen-bond angles but differing hydrogen-bond lengths; PBE produces hydrogen bonds that are too weak and WC produces hydrogen bonds that are slightly too strong. The difference in the calculated molar volumes between the two GGA functionals therefore arises principally from an average decrease of 4 ± 2% in the length of H⋯O hydrogen bonds.

The calculations capture differences in S—O bond lengths due to differences in the number of hydrogen bonds accepted by the apical O atoms (notably, S—O2), a feature that is not seen clearly in the experimental powder refinement. This is indicative of a level of `noise' in the structural parameters arising from the multi-phase powder refinement, since these subtle differences are typically seen very clearly in single-crystal data from other similar materials (e.g. Fortes et al., 2015). The shape of the M2 octahedron is well reproduced, being elongated along the Mg2—Ow5 axis as a result of the hydrogen bond accepted by Ow5. However, both the PBE and WC calculations reveal a similar distortion of the M1 octahedron, this being elongated along the Mg1—Ow3 axis, for which there is no obvious rationale since Ow3 is quite unambiguously in trigonal coordination, i.e. it does not accept a hydrogen bond. In this instance the neutron powder refinement shows the `expected' pattern of bond lengths for the M1 octahedron.

Since the overall structural parameters produced by the WC GGA calculations are in such excellent agreement with the observations a choice was made to proceed with further work using only the WC functional.

Geometry optimizations were carried out on MS9 at a range of fixed pressures from −2 to + 5 GPa. Fig. 4 shows a plot of the total energy of the crystal as a function of unit-cell volume in this pressure range; zero-pressure corresponds to the minimum in the curve, where the gradient is zero, since P = −dE/dV.

Figure 4 Calculated variation of total electronic energy against unit-cell volume for the experimentally observed structure of MgSO4·9H2O (α-MS9) and a hypothetical polymorph with the same structure as MgSeO4·9H2O, dubbed β-MS9. Solid and dashed lines are integrated third-order Birch–Murnaghan equations of state fitted to the calculation results (open and filled symbols). For α-MS9 the fit parameters are V0 = 1156.4 (5) Å^3, K0 = 24.2 (6) GPa and K′0 = 3.7 (5). For β-MS9 the equivalent fit parameters are V0 = 1177.1 (4) Å^3, K0 = 25.6 (5) GPa and K′0 = 4.4 (3).

We also evaluated the stability of a possible sulfate analogue to MgSeO₄·9H₂O. The selenate is stable in aqueous solutions at low temperature, allowing growth of large single crystals, but it is not isotypic with the sulfate described in this paper. However, it is trivial to carry out geometry optimizations on this structure with sulfur substituted for selenium. As shown in Fig. 4, this hypothetical polymorph (termed β–MS9) is marginally less dense than the observed polymorph (α-MS9), having a molar volume at zero pressure (the minimum in the dashed curve) that is ∼ 1.8% larger, and it has a lower total electronic energy at volumes greater than 1020 ų per unit cell. More specifically what this means is that the α-MS9 structure is not the thermodynamically stable form of MgSO₄·9H₂O at ambient pressure; there is an alternative structural arrangement with the same composition (β-MS9) having a lower enthalpy at 0 K. The enthalpy difference between the two phases at zero pressure amounts to 17 (1) kJ mol⁻¹ (in favour of β-MS9) and the point at which the enthalpy difference, ΔH = 0, occurs is at 4.3 (4) GPa. In §3.4 we will also show that α-MS9 is unstable with respect to both MgSO₄·11H₂O and MgSO₄·7H₂O, depending on the temperature and availability of water vapour, and we must also assume that there is a sufficient kinetic barrier even at 250 K to prevent the transformation of α-MS9 into the hypothetical β-phase. A CIF containing the zero-pressure structural parameters of β-MS9 is provided in the supporting information.

3.3. Thermal expansion

Lattice parameters were refined at 45 temperatures from 9 to 230 K, in 5 K increments (Series 1), nine temperatures from 220 to 260 K in 5 K increments (Series 2), and 24 temperatures from 9 to 260 K in 10 K increments (Series 3). These values are given in Table S1 and plotted in Fig. 5.

Figure 5 Unit-cell parameters of MS9 as a function of temperature; solid lines represent unweighted fits of a modified Einstein expression [equations (1), (2) and (3)] to the Series 1 data, done using Origin Pro (OriginLab, Northampton, MA).

It is quite obvious that not all of the cell-parameter data are continuous. This is most clear in the b-axis data, although both a and c exhibit more subtle signatures, where a relaxation process occurs on warming through ∼ 220 K. In Series 1, this effect appears to be absent, with all parameters forming smooth curves that are, importantly, contiguous with the higher-temperature portions of Series 2 and 3. We infer that the Series 1 cell parameters represent the most completely annealed strain-free state of the crystal at all temperatures. Series 3 has a significantly different b-axis length to the other observations, which `turns over' at 220 K to join Series 1 and 2. Series 2 exhibits a similar small `kink' at the lowest temperature; taken together these seem to show a structural relaxation along the unique axis of the crystal into a more stable state. The temperature at which this structural relaxation occurs is in the same region where the progenitor glass devitrifies. The origin of the phenomenon is not clear, since we have no high-precision structural datasets in this region and there are no orientationally disordered water molecules in the structure. One possible explanation is that Bjerrum defects (such as two protons occupying the Ow2⋯Ow5 bond) could lead to expansion along b, although we have no data with which to test this hypothesis.

For the purpose of making a simple density calculation (e.g. for planetary modelling) we have fitted a second-order polynomial to the calculated density of MS9 between 80 and 260 K, ρ = AT² + BT + ρ₀ assuming that the molar volumes of the protonated and deuterated species are identical. The coefficients obtained are A = −3.25 (5) × 10⁻⁷ g cm⁻³ K⁻², B = −2.4 (2) × 10⁻⁵ g cm⁻³ K⁻¹ and ρ₀ = 1.6095 (1) g cm⁻³ for protonated MS9, valid only between 80 and 260 K.

In order to interpret the subtleties of the thermal expansion tensor, we analyse the data in terms of an Einstein oscillator model similar to that used previously for analysis of the thermal expansion of MS11 (Fortes et al., 2008) where all atoms in the solid vibrate at the same characteristic angular frequency, ω_E, which is expressed in terms of an Einstein temperature, θ_E = - hω_E/k_B. The advantage of this with respect to the more complex double-Debye model required to fit these data accurately is that fewer adjustable parameters are needed, as shown below.

The temperature dependence of the molar volume, V(T) is described as

where E = 3Rγθ_E/K_T (γ is the Grüneisen ratio, K_T is the isothermal bulk modulus and V₀ is the zero-temperature molar volume). Note that this is only dimensionally correct when V is a volume in m³ mol⁻¹.

A similar expression may be used to describe the temperature dependence of the lattice parameters

where X₀ is the value of the fitted parameter at 0 K. For the purpose of deriving the coefficients of the thermal expansion tensor a sufficiently good fit is only obtained when the parameter E is allowed to vary linearly as a function of temperature

The parameters obtained from fitting equation (1) [and (3)] to the unit-cell volume of MS9, and fitting equation (2) [and (3)] to the individual cell parameters, are listed in Table 3 and the fits are depicted in Fig. 5 as solid lines.

Table 3
Einstein fit parameters [equations (4) and (6) for the volume, (5) and (6) for cell parameters]

	a-axis	b-axis	c-axis	β	V
X0 (Å, Å^3)	6.7274 (1)	11.9039 (1)	14.6404 (1)	95.199 (1)	1167.63 (1)
θ (K)	189 ± 15	153 ± 7	439 ± 14	59 ± 15	209 ± 5
E0	1.8 (5) × 10^−2	2.9 (3) × 10^−2	2.5 (3) × 10^−1	−1.4 (4) × 10^−2	12.5 (8)
E1	1.10 (4) × 10^−4	1.25 (2) × 10^−4	−2.7 (7) × 10^−4	–	4.21 (7) × 10^−2

This formalism is not as accurate as the more rigorous Debye model, which would be the preferred choice when heat capacity data are available to aid in constraining the Debye temperatures. Nevertheless, we present the fitting of two different Debye models in the supporting information.

Since for a crystal of monoclinic symmetry, two of the three principal directions of the thermal expansion tensor are free, whilst remaining orthogonal to one another, to rotate freely about the crystal's unique axis, an accurate picture of the structure's behaviour is not readily obtained by only fitting expressions to the individual axial lengths as a function of temperature. We therefore use the Einstein model fits to the cell parameters in order to derive the thermal expansion tensor coefficients. For a monoclinic crystal these are represented by a symmetrical second rank tensor of the form

the components of which are unit strains corresponding to thermal expansion coefficients. The eigenvalues and eigenvectors of this matrix, obtained by matrix decomposition methods, are the magnitudes and orientations of the three principal axes of the thermal expansion tensor (i.e. directional expansivities, , and ) with respect to an orthogonal basis. We have applied the commonly used Institute of Radio Engineers relationship between the orthogonal basis and the unit-cell of a monoclinic crystal such that X || a*, Y || b and Z || c (Boisen & Gibbs, 1990).

The fitted parameters in Table 3 were used to calculate smoothly varying unit-cell parameters as a function of pressure, from which the thermal expansion coefficients were obtained using the methods described by Schlenker et al. (1978) and Hazen & Finger (1982); the magnitudes of the principal expansivities, , and , and their spatial orientation with respect to the Cartesian basis were then obtained by standard eigenvalue decomposition methods.

The principal and volumetric thermal expansivities are shown in Fig. 6; symbols with error bars depict point-by-point derivatives of the `raw' unit-cell parameters calculated (in order to reduce the uncertainties) across moving 15 K wide windows.

Figure 6 Principal and volume thermal expansion coefficients of MS9 as a function of temperature. Solid lines are derived from smoothed values calculated using the Einstein-model fit parameters reported in Table 3. Symbols are `point-by-point' derivatives of the raw cell parameter values calculated across a moving window 15 K wide. The dashed red line in (d) depicts the volume thermal expansion of MS11 for comparison with that of MS9: the dotted blue line represents the unconstrained double-Debye model fit to the data (see text in the supporting information)

Representation surfaces (Reynolds glyphs – see Hashash et al., 2003) of the unit-strain tensor – or more specifically in this case, the thermal expansion tensor – at a range of temperatures are shown in Fig. 7, where the distance from the origin to the edge of the representation surface indicates the thermal expansivity in any given direction.

Figure 7 Representation surface of the MS9 thermal expansion tensor at selected temperatures. Coloured arrows identify the eigenvectors of the thermal expansion tensor and the black arrows show the crystallographic axes. Green shading indicates positive values and red shading indicates negative thermal expansion values. Tensor surfaces were generated in WinTensor (Kaminsky, 2004 2007) and postprocessed with Meshconv (courtesy of Patrick Min, Princeton) and Meshlab (Sourceforge).

The absolute values of the volume thermal expansion coefficient are large when compared with other similar cryohydrates. The temperature dependence of the volume expansivity of MS11(D) is drawn in Fig. 6(d); for comparison the volume expansivity of MS11(D) at 250 K is ∼ 72 × 10⁻⁶ K⁻¹, that of MS7(D) is ∼ 89 × 10⁻⁶ K⁻¹ (Fortes et al., 2006) and that of Na₂SO₄·10D₂O is ∼ 90 × 10⁻⁶ K⁻¹ (Brand et al., 2009). However, the volume thermal expansivity of MS9 is virtually identical to the non-isostructural MgSeO₄·9D₂O, these being 113 and 110 × 10⁻⁶ K⁻¹, respectively, at 250 K and the two materials exhibit large degrees of anisotropy (Fortes, 2016). For MS9(D), at 250 K, the eigenvalues of the thermal expansion tensor are fairly isotropic in the a–b plane (/ = 0.88) whilst the eigenvalue most closely aligned with c is substantially smaller (/ = 0.35). In MgSeO₄·9D₂O at 250 K the direction of least thermal expansion is aligned with the twofold axis, which we define as ; hence, /= 0.09 and / = 0.63. It is also noteworthy that both sulfate and selenate 9-hydrates lack the region of negative volume expansion observed in MgSO₄·11D₂O and Na₂SO₄·10D₂O.

In MS9 the angle between and c is negligibly small at high temperatures and only increases to > 5° below 100 K. A fairly dramatic tilt of the representation figure below 50 K corresponds to turning negative and the development of a cone of pure shear between the positive and negative lobes of the strain figure.

3.4. Phase transitions at constant T and on warming

In the vast majority of MS9 syntheses, water ice is present in substantial quantities. We observed, however, that the intensity of ice peaks would diminish during the acquisition of X-ray powder diffraction data at ca. 250 K at the same time as which additional Bragg peaks from MS11 would appear. For measurements of many hours duration, this would continue until virtually all detectable ice in the specimen had been consumed. Our interpretation of this was that MS9 reacts with water vapour in pore spaces and at the specimen's large air interface, to form the thermodynamically stable phase, MS11. Since the equilibrium vapour pressure of ice is quite substantial, approaching 1 mbar at the temperature of the X-ray measurements (92.4 Pa at 252 K: IAPWS, 2011; Wagner et al., 2011), the vapour consumed by reaction with MS9 must be replaced by further sublimation of ice. Consequently, the rate at which this reaction proceeds, being a function of the vapour pressure of ice, will be strongly dependent on temperature. Moreover, at a given temperature, it should be slower for deuterated specimens than protonated samples by virtue of the low vapour pressure over D₂O ice (P_vap[D₂O] = 67 Pa at 252 K: Matsuo et al., 1964).

A series of X-ray powder diffraction patterns were obtained from a sample containing MS9 + ice with negligible initial quantities of MS11 over a period of several hours, measuring in the 2θ range 5–50° every 20 min (Fig. S5); clearly, peaks from MS11 grow over time at the expense of both MS9 and water ice.

These XRD data were used to refine the phase fractions of each constituent as a function of time (Fig. 8b). The well known Kolmogorov–Johnson–Mehl–Avrami (KJMA) expression (Kolmogorov, 1937; Johnson & Mehl, 1939; Avrami, 1939, 1940) was fitted to these data

where X is the MS11 phase fraction, A is the saturation value of the phase fraction, k is a rate constant (units time⁻ⁿ), and n is the `Avrami exponent', which contains information on the time dependence of nucleation, the geometry of the growing crystallites and the nature of the reaction process (interfacial versus diffusional). The Avrami exponent n is the sum of two quantities, α + β. The number density of nucleation centres (N) has a model time dependence of the form N ∝ t^α, such that for α = 0 all nucleation sites are present at t = 0 (the growth medium is nuclei site saturated), and for α = 1 the nucleation rate is a constant: for α < 1 the nucleation rate slows with time, and for α > 1, the nucleation rate increases with time. The term β expresses the dimensionality of the growth geometry (i.e. one-dimensional = needle-like or acicular growth, two-dimensional = platy or tabular growth, three-dimensional = blocky or globular growth), and has the values 1, 2 or 3 for one-dimensional, two-dimensional or three-dimensional growth, respectively, when the reaction process is interfacial, and has the values 0.5, 1 or 1.5 when the reaction process is diffusional. For our purposes, we can infer to a certain extent the nature of the reaction process, and also use the known habit of MS11 to infer that the growing crystallites will be platy.

Figure 8 (a) Refined phase fractions of MS9, MS11 and ice as a function of time from X-ray powder diffraction data measured on a Peltier-cooled cold stage with ca 52 cm^3 head-space over the sample. The solid line through the symbols for MS11 reports the fit of equation (5) to the data. The two dashed lines show the expected phase fractions of ice and MS9 concomitant with the observed abundance of MS11 as a function of time. (b) Refined phase fractions of deuterated MS9, MS11 and ice as a function of time from neutron powder diffraction data measured in an indium-sealed container. Lines and symbols have the same meaning as in part (a).

The data in Fig. 8(a) are fitted with equation (5) to obtain A = 0.68 (1), k = 2.9 (4) × 10⁻⁷ s⁻ⁿ and n = 1.73 (8). Thus, from the X-ray powder data collected at 252 K, we find that the peak growth rate is 1.3 × 10⁻⁴ s⁻¹, which occurs after 1.0 h. It is clear from Fig. 8(a) that the consumption of MS9 is smaller, and that of ice larger, than is commensurate with the observed abundance of MS11. This suggests that, in the thin surface layer of the sample being probed by X-rays, water ice is being preferentially depleted by another process – most probably sublimation: in other words, the observed phase mixture is not determined only by the reaction MgSO₄·9H₂O + 2H₂O → MgSO₄·11H₂O.

In fact we have the possibility to test this hypothesis by observing the same process in a confined bulk sample using a neutron diffraction probe (Series 3). Fig. S6 shows the growth of MS11 Bragg peaks at the expense of MS9 and ice at 260 K over a period of 5.5 h in an indium-sealed can. As before, the refined phase fraction of MS11 has been fitted with equation (5) to obtain A = 0.41 (1), k = 6.0 (2) × 10⁻⁴ s⁻ⁿ and n = 0.83 (4). As shown in Fig. 8(b), we now find that the expected abundances of both MS9 and ice are in closer agreement with the observations, reflecting the limited pore space and head space in the sample can that is available for saturation with vapour from sublimated water ice.

Evidently, the peak growth rate has passed at the onset of the 260 K neutron measurements, and the rate we observe declines as a function of time. This suggests that growth began (and peaked) during the phase of warming from 250 to 260 K and/or the subsequent short period of thermal equilibration. In terms of the Avrami exponent, the most logical interpretation is that growth of MS11 is two-dimensional (platy) and diffusional, giving β = 1. The difference between the two observations then is that nucleation occurs slowly over the course of the X-ray measurement at 252 K (i.e. α ≃ 1), whereas the slightly higher temperature in the neutron measurement [and associated increase in vapour pressure, P_vap(D₂O) = 147 Pa] results in rapid saturation of nucleation centres (i.e. α ≃ zero). This produces n ≃ 2 for the X-ray measurements and n ≃ 1 for the neutron measurements, close to what we obtain and consistent with both the expected morphology and crystal growth process.

In light of the transformation of MS9 to MS11 in the presence of ice, or more accurately water vapour in equilibrium with ice, we made considerable efforts to synthesize ice-free samples. As described in §2.1, we were able to achieve this by careful evaporation of MgSO₄ solutions into a supersaturated state at ∼ 383 K followed by flash freezing in liquid nitrogen. Fig. S7 shows X-ray powder diffraction patterns obtained from just such an ice-free specimen both before and after a 12 h data collection. The objective was to form an initial liquid containing 42.6 wt % MgSO₄, which would freeze to form phase-pure MS9. The XRD data reveal a small amount of MS11 is present with a refined phase fraction of 4.8 (1) wt %; this corresponds to a bulk composition for the sample of 42.4 wt % MgSO₄ – very close indeed to our target.

After 12 h, the accessory MS11 has disappeared and we see instead Bragg peaks from MS7, epsomite, with a refined phase fraction of 15.7 (1) wt %. The bulk composition of the sample is now 43.6 wt % MgSO₄, indicating a small net loss of water. The observation requires all of the initial MS11 and ∼ 11% of the MS9 to have dehydrated to MS7 in the cold, dry conditions of the sample holder with the loss (depending on initial packing density) of between 2 and 4 mg of water. This is sufficient, if confined solely to the sample pore space, to produce a partial pressure of water vapour of some tens of Pascals. It would be useful to characterize this process in more detail, both as a function of temperature and relative humidity so as to place this process in context with what we know about the transformation of other MgSO₄ hydrates under similar conditions (Chou & Seal, 2004; Chipera & Vaniman, 2007; Wang et al., 2009, 2011). Such information is essential to determining the longevity of MS9 as a mineral in, for example, the martian regolith.

3.5. Isothermal equation of state

Data Series 4, measured in a TiZr pressure vessel under He gas, yielded unit-cell parameters for MS9(D) at six state points between 10 and 543 MPa at 240 K. One additional datum at 1.07 (2) GPa was obtained by refinement of a neutron powder pattern from a mixture of MS9(D) and deuterated ice VI formed by the dissociation of MS11(D) (see Fortes et al., 2017). Cell parameters as a function of pressure are given in Table S2.

A convenient description of the material's stiffness may be obtained by parameterizing the pressure-dependence of the unit-cell volume using a Murnaghan integrated linear equation of state, MILEOS (Murnaghan, 1944)

where K₀ is the zero-pressure isothermal bulk modulus and K′ is the first pressure derivative of the bulk modulus, (∂K₀/∂P). An unweighted fit of equation (6), done using OriginPro, to the experimental V(P) obtained solely from the gas cell study on HRPD yields V₀ = 1184.9 (2) ų, K₀ = 19.3 (6) GPa and K′ = 4 (2): inclusion of the data point from the higher-pressure study on PEARL/HiPr improves the precision most noticeably on the value of K′, with V₀ = 1184.8 (2) ų, K₀ = 19.5 (3) GPa and K′ = 3.8 (4). Hence the incompressibility of MS9(D) is virtually identical to synthetic meridianiite (K₀ = 19.9 GPa; Fortes et al., 2017) and synthetic mirabilite (K₀ = 19.1 GPa: Brand et al., 2010) and MS9(D) is somewhat more compressible than synthetic epsomite (K₀ = 22.6 GPa: Fortes et al., 2006).

As we have done for the thermal expansion, the elastic strain due to a hydrostatic stress is described by a symmetrical second rank tensor of the form

the components of which are unit strains corresponding to compressibility coefficients. The eigenvalues and eigenvectors of this matrix, are the magnitudes and orientations of the three principal axes of the elastic strain figure (i.e. directional compressibilities, , and ) with respect to the same orthogonal basis as was used in the derivation of the thermal expansion tensor.

Each of the unit-cell parameters was fitted with a Murnaghan integrated linear equation of state of the same form as equation (6): the parameters resulting from these fits are listed in Table 4. A plot of the high-pressure unit-cell parameters is given in Fig. 9 with solid and dashed lines reporting various MILEOS fits to each.

Table 4
Murnaghan integrated linear equation of state [equation (6)] fit parameters

	a-axis	b-axis	c-axis	β	V
MILEOS parameters including PEARL/HiPr datum (240 K)
X0	6.7701 (7)	11.9705 (9)	14.6799 (4)	95.106 (6)	1184.8 (2)
K0	51 (2)	39.9 (6)	139 (2)	−132 (8)	19.5 (3)
K′	40 (3)	10 (1)	−60 (3)	−292 (15)	3.8 (4)
MILEOS parameters excluding PEARL/HiPr datum (240 K)
X0	6.7694 (6)	11.9708 (12)	14.6801 (5)	95.112 (3)	1184.9 (2)
K0	55 (3)	39 (1)	135 (5)	−154 (7)	19.3 (6)
K′	24 (10)	12 (5)	−44 (18)	−195 (27)	4 (2)
MILEOS parameters (athermal WC GGA calculations)
X0	6.625 (1)	11.963 (4)	14.633 (2)	94.12 (2)	1156.9 (4)
K0	61 (1)	58 (1)	140 (2)	−287 (17)	24.2 (2)
K′	17 (1)	7 (1)	15 (2)	−49 (16)	4.1 (2)

Figure 9 Unit-cell parameters of MS9 as a function of pressure at 240 K. Open symbols correspond to data measured in a TiZr gas cell on HRPD (Series 4), whereas the filled symbol was measured in a Paris–Edinburgh press on PEARL/HiPr (Fortes et al., 2017). Solid lines are least-squares fits of equation (6) to all of the data, whereas the dashed red lines are fits only to the gas cell data. The fitted parameters are listed in Table 4. Panel (f) summarizes the relative variation of each unit-cell parameter with pressure. (g) Principal linear compressibilities and (h) volume compressibility of MS9 at 240 K. Solid lines are derived from fitting equation (6) to both the HRPD and PEARL/HiPr data, whereas the dashed lines are derived from fitting only to the HRPD gas-cell data. Symbols with error bars show the moving-average derivatives calculated from the observed unit-cell parameters. Panel (i) shows the angle between the stiffest direction in the crystal, and the c-axis.

The fitted parameters in Table 4 were used to calculate smoothly varying unit-cell parameters as a function of pressure, from which the Eulerian infinitesimal unit-strain coefficients were obtained using the method outlined in §3.3; the magnitudes of the principal compressibilities, , and , and their spatial orientation with respect to the orthogonal basis were then obtained by standard eigenvalue decomposition methods. Fig. 9(g) shows the linear compressibilities derived from the raw experimental data and from the MILEOS fits where, as in the rest of Fig. 9, dashed lines correspond to fits only to the gas cell data and solid lines include the PEARL/HiPr point at 1.07 (2) GPa.

The principal and volumetric compressibilities and incompressibilities (e.g. ) were fitted with linear functions in P (in the low-pressure limit) to obtain the zero-pressure values and their first derivatives (Table 5). It is worth observing that this method recovers with reasonable accuracy the values of K₀ and K′ found earlier by fitting the Murnaghan integrated linear equation of state. It is apparent from Fig. 9(g) that the stiffest principal direction, , softens under compression – indeed this can be seen by visual inspection of the c-axis as a function of pressure (Fig. 9c). What Fig. 9(g) shows is that the softening of is reduced above ∼ 1 GPa whilst the slope of (which had hitherto been stiffening at a high rate) turns positive. The overall effect of this is to cause the bulk compressibility of the crystal to increase above ∼ 1 GPa (Fig. 9h). This change in the eigenvalues of the strain tensor is matched by a large change in the eigenvectors; Fig. 9(i) shows that the angle between and the c-axis increases by ∼ 50° from 0 to 1.2 GPa.

Table 5
Derived linear and volume compressibilities (B) and incompressibilities (K) found by fitting linear expressions to the pressure dependences of the elastic strain tensor's eigenvalues (HRPD + PEARL/HiPr)

Values in italics are those found from DFT calculations (athermal, WC GGA).

	Compressibility		Incompressibility
	B0 (GPa^−1)	B′	K0 (GPa)	K′
	2.24 × 10^−2	−1.92 × 10^−2	44.6	40.9
	1.75 × 10^−2	−4.6 × 10^−3	57.3	14.9
	2.50 × 10^−2	−6.0 × 10^−3	40.0	9.8
	1.72 × 10^−2	−1.9 × 10^−3	58.2	6.5
	5.2 × 10^−3	7.1 × 10^−3	190.8	−234.0
	6.5 × 10^−3	−8.0 × 10^−4	153.4	18.8
βV	5.27 × 10^−2	−1.81 × 10^−2	19.0	6.7
	4.12 × 10^−2	−7.3 × 10^−3	24.3	4.3

Representation surfaces of the elastic strain figures at zero pressure, 500 and 900 MPa are shown in Fig. 10 where the distance from the origin to the edge of the representation surface indicates the compressibility in any given direction under hydrostatic stress.

Figure 10 Representation surface of the MS9 unit-strain figure at selected pressures. Coloured arrows identify the eigenvectors of the compressibility tensor and the black arrows show the crystallographic axes. Tensor surfaces generated in WinTensor (Kaminsky, 2004, 2007) and post-processed with Meshconv (courtesy of Patrick Min, Princeton) and Meshlab (Sourceforge).

There are some noteworthy similarities between the strain response of the crystal to both thermal and hydrostatic stress. The crystal in both instances exhibits the largest fairly isotropic strains in a plane that is roughly perpendicular to c whilst exhibiting the smallest strain roughly parallel with c.

3.6. DFT equation of state

The energy–volume curves (Fig. 4) were each fitted with an integrated form of the third-order Birch–Murnaghan equation of state in order to determine the zero-pressure athermal molar volume, the zero-P/T bulk modulus and the first pressure derivative of the bulk modulus. The parameters from these fits are given in the caption to Fig. 4. Since our calculations also provide us with relaxed cell-parameters at each pressure point (Figs. 11a–d) we are able to reproduce the derivation given in §3.5 of the elastic strain tensor coefficients. The pressure dependence of the principal axial compressibilities and the volume compressibility are shown in Figs. 11(g)–(i). Parameters evaluated at zero P and T are reported in italics in Table 5 for comparison with the experimental values at 240 K.

Figure 11 Unit-cell parameters of MS9 as a function of pressure from DFT calculations (athermal, WC GGA). Solid lines are least-squares fits of equation (6), the fitted parameters being listed in Table 4. Panel (f) summarizes the relative variation of each unit-cell parameter with pressure. (g) Principal linear compressibilities and (h) volume compressibility of MS9 from DFT calculations (athermal, WC GGA). Panel (i) shows the angle between the stiffest direction in the crystal, , and the c-axis.

The pressure dependences of the unit-cell parameters are in very good agreement with the experimental observations – with one important exception; they do not reproduce the increasing compressibility of the c-axis and the associated increase in bulk compressibility. This difference may be due entirely to the influence of temperature. However, obtaining further insight through finite-temperature calculations represents a formidable computational expense. Our experimental observation of the pressure-induced softening seems robust, being present in data both with and without the additional point obtained in the Paris–Edinburgh press. An improved high-pressure experimental study that straddles the two existing datasets and also samples a wide temperature space is essential before an ab initio molecular dynamics study can be justified.

4. Possible natural occurrences of MgSO₄·9H₂O

Magnesium sulfate is a quite common constituent of water on Earth, although the number of locations where it is sufficiently concentrated to produce MgSO₄ hydrate (or cryohydrate) minerals is not large, the stability of the higher hydrates being dependent on low ambient temperatures and high humidity. The highest hydrate, MgSO₄·11H₂O, occurs naturally as the mineral meridianiite in a variety of glacial and periglacial environments (Sakurai et al., 2009; Genceli et al., 2009) and in a limited number of MgSO₄-rich hypersaline lakes during the winter months; example localities where the mineral has been identified include the Basque Lakes, Clinton Lake and kłlil'x^w (aka Spotted Lake), all in British Columbia, Canada (e.g. Peterson et al., 2007; Cannon, 2012). Whilst MgSO₄-rich saline waters are comparatively rare on Earth due to the influence of continental weathering, such liquids are expected to be common on other rocky planets where the weathering of basaltic materials dominates (King et al., 2004). On Mars, abundant magnesium(II) and iron(III) sulfates are known to occur, including minerals such as kieserite and jarosite (e.g. Clark et al., 1976; Toulmin et al., 1977; Wänke et al., 2001; Foley et al., 2003; McSween, 2004; Chipera & Vaniman, 2007) and it is hypothesized that meridianiite may occur in a permafrost-like deposit, forming a substantial reservoir of bound water in the near-surface regolith (Feldman, Mellon et al., 2004; Feldman, Prettyman et al., 2004; Peterson & Wang, 2006). Similarly, water–rock interactions during the accretion and differentiation of icy planetary bodies in the outer solar system may have resulted in large brine reservoirs crystallizing substantial quantities of MgSO₄ and Na₂SO₄ hydrates (Kargel, 1991). These are apparent in near-IR spectra of their surfaces (Orlando et al., 2005; Dalton, 2007; Dalton et al., 2005; Shirley et al., 2010), although it remains unclear the extent to which some of the hydrated salts on Europa's surface are due to endogenic versus exogenic processes, such as radiolysis of MgCl₂ combined with sulfur implantation from neighbouring Io (Brown & Hand, 2013).

Although MS9 is evidently metastable in the presence of liquid water or water vapour at ambient pressures we have determined that it forms above ∼ 800 MPa at 240 K by the decomposition of MS11 into MS9 + ice VI (Fortes et al., 2017) in much the same fashion as MgSiO₃ perovskite in the Earth's mantle forms when spinel-structured Mg₂SiO₄ (ringwoodite) experiences a pressure-induced exsolution of one formula unit of its component oxide. Judging from the inferred high-pressure melting curve of MS11 derived in our companion paper, MS9 may well occur at pressures of 600–800 MPa in the region from 240 to 270 K. The implication of this is that MS9 could be present in substantial quantities as a `rock-forming' mineral in icy satellites that are speculated to have MS11-bearing outer layers of sufficient depth (e.g. Jupiter's moon Ganymede).

Furthermore, the manner in which we synthesize MS9 suggests another way in which this otherwise metastable hydrate could form and persist in nature. Salt hydrates on the surfaces of icy satellites may have been emplaced by cryovolcanic eruptions, perhaps involving production of droplets or spatter and perhaps with sufficient exsolution of volatiles to create a foam or fine spray. Under such circumstances it is easy to envisage the rapid chilling of fine droplets in the cold and very low-pressure environment (ca 100 K, few μPa) of an icy satellite's surface leading to solidification of glassy beads, similar to the orange and green glassy volcanic spherules found in the lunar regolith (e.g. Heiken et al., 1974). Annealing processes acting over gigayear timescales, or shorter periods in contact with `warm' cryovolcanic deposits could conceivably lead to crystallization of MS9 and/or other as-yet unknown metastable cryohydrates from these glassy beads.

Lastly, it is possible that the presence of other solutes could ameliorate the growth and persistence of MS9, either by providing a nucleation template or by shifting the thermodynamics of the system in its favour. The possibility that this might occur was identified serendipitously: in broader work on MgSO₄-rich hypersaline analogue environments, we had cause to collect brine samples on two occasions from Spotted Lake, near Osoyoos, British Columbia, Canada (Fig. 12).

Figure 12 View of Spotted Lake looking north-west acquired by ADF in September 2013. The lake is ∼ 700 m long and 240 m across at its widest point, having an area of 13 hectares. The bulk of this endorheic basin is filled with coarse gypsum crystals and organic mud. Concentrated Mg^2+–Na^+–SO4^2− brine occupies some 670 discrete pools from < 1 m to ∼ 40 m across. The sample described in the text was obtained at 49° 4′ 36.98″ N, 119° 33′ 55.62″ W.

The first sample was collected in late September (air temperature = 292 K), when the brine pools were precipitating large circumferential `reefs' of water-clear mirabilite, thereby depleting the remaining liquid of sodium. The second sample was collected in late December of the same year (air temperature = 271 K), when the pools were roofed over with a layer of quite soft, wet ice ∼ 3 cm thick, beneath which was liquid. Brine samples were dried at 673 K for 24 h to determine the total dissolved solid content and the residue was then analysed for major elements using a Jeol JXA8100 microprobe. In both samples, the brine composition was dominated by Mg²⁺, Na⁺ and SO₄²⁻ with very small quantities of K⁺ and Cl⁻ (∼ 0.5 and 0.2 mol L⁻¹, respectively, in both samples). The autumn brine (concentration 26.6 wt %) had a Mg/Na ratio of 0.9, whereas in the winter brine (concentration 15.3 wt %) this had increased significantly to Mg/Na = 3.0, reflecting the precipitation of mirabilite throughout the summer and autumn. The remainder of the winter brine (∼ 30 g) was decanted into a plastic bottle and left in a chest freezer for several hours. At the end of this period the contents had solidified entirely. These were extracted, ground to a powder, and examined by X-ray powder diffraction on our Peltier cold stage.

Fig. 13 shows the diffraction pattern obtained from this specimen, which is clearly an excellent match for a mixture of MS9, mirabilite and water ice. Since our treatment of this brine is no different from what may occur during a cold snap at any of the MgSO₄-rich lakes in British Columbia, we venture to suggest that MS9 forms in these environments, albeit on a transient basis, as a natural material. It also suggests that there may be a route to preparing single crystals by equilibrium growth from an appropriately doped aqueous solution. Testing the hypothesis that MS9 can occur naturally in cold hypersaline lakes requires collection, transport and analysis of frozen specimens ideally at liquid nitrogen temperatures or else in situ analysis with a portable X-ray diffractometer as long as the specimen is kept below 250 K during the measurement. In situ Raman spectroscopy would be the easiest method of identifying MS9 and work is in progress to characterize diagnostic differences between the vibrational spectra of MS7, MS9 and MS11. These efforts will pay dividends in future efforts to unambiguously identify cryohydrates by Raman spectroscopy done in situ on planetary surfaces and at very large stand-off distances, i.e. from orbit (Angel et al., 2012).

Figure 13 Top: X-ray powder diffraction data (λ = 1.788996 Å) from a slowly frozen specimen of winter sub-ice brine from Spotted Lake, British Columbia. Bottom: Calculated X-ray powder diffraction patterns of MS9 (blue) and mirabilite (red). Strong peaks from water ice are shaded grey. The phase fractions obtained by Rietveld refinement are 23.0 (5) wt % MgSO4·9H2O, 7.3 (5) wt % Na2SO4·10H2O, 69.7 (5) wt % ice.

It is worthwhile commenting on a recently published paper by Vu et al. (2016) in which various brines thought to be representative of the ocean composition of Jupiter's icy satellite Europa were frozen and examined by Raman microscopic methods. These authors observed a prevalence of mirabilite over MgSO₄ hydrates in all brine compositions, which ranged from Mg/Na ratios of 0.47 to 1.64. In the most Mg-rich frozen brines, Vu et al. report (their Fig. 5) a ν₁ sulfate band at 990 cm⁻¹, which they attribute solely to mirabilite. However, this band is quite asymmetric on the low-frequency side, as one would expect if there were a substantial unresolved shoulder near 980–985 cm⁻¹ from a MgSO₄ cryohydrate. Secondly, these authors could not be aware (reasonably) of either the existence or the Raman signature of a possible MgSO₄·9H₂O. Finally, their Mg/Na ratios are substantially smaller than the Spotted Lake sample from which we obtained MS9 + mirabilite, and are also smaller than the ratio predicted for a Europan ocean by the FREZCHEM Pitzer potential calculations of Marion et al. (2005). If nothing else, this work emphasizes (i) the value of establishing definitive vibrational signatures of these cryohydrates, since the internal sulfate modes do not differ substantially amongst the higher hydrates (Wang et al., 2006) and so the less commonly measured low–frequency lattice modes or THz spectra are more likely to have interpretive value, and (ii) the importance of doing simultaneous diffraction and vibrational spectroscopy to provide certainty about the phases or phase mixtures being studied.

5. Concluding remarks

We have determined the crystal structure of MgSO₄·9H₂O, including all hydrogen positions using a combination of X-ray and neutron powder diffraction. Additionally, we have determined the thermal expansion and incompressibility tensors of MgSO₄·9H₂O over the range 10–260 K at ambient pressure and 0–1 GPa at 240 K. Our observations have been instrumental in understanding the decomposition of meridianiite at high pressure and provide useful information on structural systematics and the transformation of metastable cryohydrates at low temperature. Further work is required to verify and further characterize an apparent elastic softening observed under compression at 240 K, but which is not reproduced in athermal ab initio calculations.

In future work we intend to characterize both the Zn²⁺ and Fe²⁺ analogues of MS9; since both of these have larger ionic radii than Mg²⁺ and different electronic structures, we expect to observe a difference in site substitution preferences and in the chemically induced strain of the unit cell. More importantly, we aim to provide the mineralogical and planetary science community with the spectroscopic data required to identify MS9 in the field – both on Earth and elsewhere.

6. Related literature

References cited in the supporting information include: Alexandrov et al. (1963), Arbeck et al. (2010), Cook-Hallett et al. (2015), Fortes, Wood et al. (2003), Geruo et al. (2014), Gowen et al. (2011), Hill (1952), Nimmo (2004), Nimmo & Matsuyama (2007), Nimmo & Schenk (2006), Reuss (1929), Tsuji & Teanby (2016), Voigt (1910), Wahr et al. (2009).

Supporting information

Computing details

Data collection: HRPD control software for 9_K_MS9. Cell refinement: GSAS/Expgui (Larsen & Von Dreele, 2000: Toby, 2001) for 9_K_MS9. Data reduction: Mantid (Arnold et al., 2014: Mantid, 2013) for 9_K_MS9. Program(s) used to solve structure: n/a for 9_K_MS9. Program(s) used to refine structure: GSAS/Expgui (Larsen & Von Dreele, 2000: Toby, 2001) for 9_K_MS9. Molecular graphics: Diamond (Putz & Brandenburg, 2006) for 9_K_MS9. Software used to prepare material for publication: publCIF (Westrip, 2010) for 9_K_MS9.

(9_K_MS9) hexaaquamagnesium(II) sulfate trihydrate

Crystal data

[Mg(D2O)6](SO4)(D2O)3	V = 1168.55 (1) Å^3
Mr = 300.60	Z = 4
Monoclinic, P21/c	Neutron radiation
Hall symbol: -P 2ybc	T = 9 K
a = 6.72764 (6) Å	Particle morphology: fine powder
b = 11.91154 (9) Å	colourless
c = 14.64240 (14) Å	Specimen preparation: Prepared at 80 K and 100 kPa
β = 95.2046 (7)°

Data collection

HRPD, High Resolution Powder diffractometer	Distance from source to specimen: 95000 mm
Radiation source: ISIS neutron spallation source	Distance from specimen to detector: 965 mm
Scan method: time of flight

Fractional atomic coordinates and isotropic or equivalent isotropic displacement parameters (Ų)

	x	y	z	Uiso*/Ueq
D9b	0.4930 (7)	0.2880 (4)	0.3579 (3)	0.0179*
D9a	0.4169 (7)	0.4054 (4)	0.3347 (3)	0.0179*
D8b	0.0921 (7)	0.2353 (4)	0.2630 (3)	0.0179*
D8a	−0.1387 (7)	0.2444 (4)	0.2620 (3)	0.0179*
D7b	0.5152 (7)	0.2628 (4)	0.1110 (3)	0.0179*
D7a	0.5267 (6)	0.3897 (4)	0.1362 (3)	0.0179*
D6b	0.1959 (7)	0.4696 (4)	0.1671 (3)	0.0179*
D6a	0.0189 (7)	0.3862 (4)	0.1658 (4)	0.0179*
D5b	0.2345 (7)	0.6577 (4)	0.0854 (3)	0.0179*
D5a	0.2690 (7)	0.6540 (4)	−0.0192 (3)	0.0179*
D4b	0.2676 (7)	0.4077 (4)	−0.1043 (3)	0.0179*
D4a	0.3518 (7)	0.3827 (4)	−0.0037 (3)	0.0179*
D3b	0.1176 (7)	0.0372 (4)	0.1852 (3)	0.0179*
D3a	0.0295 (7)	−0.0792 (4)	0.1738 (4)	0.0179*
D2b	0.2306 (7)	0.1692 (4)	−0.0397 (4)	0.0179*
D2a	0.0398 (7)	0.2246 (4)	−0.0121 (3)	0.0179*
D1b	0.3723 (7)	−0.0907 (4)	0.0349 (3)	0.0179*
D1a	0.3298 (7)	−0.0800 (4)	−0.0708 (3)	0.0179*
Ow9	0.4315 (7)	0.3527 (4)	0.3855 (3)	0.00106*
Ow8	−0.0178 (7)	0.2901 (3)	0.2550 (3)	0.00106*
Ow7	0.5354 (7)	0.3375 (4)	0.0842 (3)	0.00106*
Ow6	0.0740 (6)	0.4488 (4)	0.1325 (3)	0.00106*
Ow5	0.1670 (7)	0.6485 (4)	0.0243 (3)	0.00106*
Ow4	0.2460 (7)	0.4217 (4)	−0.0403 (3)	0.00106*
Mg2	0.0	0.5	0.0	0.00106*
O4	0.6340 (6)	0.1806 (3)	0.3073 (3)	0.00106*
O3	0.5643 (6)	−0.0146 (3)	0.2732 (3)	0.00106*
O2	0.5130 (5)	0.1254 (3)	0.1533 (2)	0.00106*
O1	0.2926 (5)	0.1176 (4)	0.2746 (3)	0.00106*
S1	0.5008 (7)	0.1032 (4)	0.2517 (3)	0.00106*
Ow3	0.0402 (7)	−0.0085 (4)	0.1407 (3)	0.00106*
Ow2	0.1244 (7)	0.1594 (4)	0.0006 (3)	0.00106*
Ow1	0.2713 (7)	−0.0691 (4)	−0.0135 (3)	0.00106*
Mg1	0.0	0.0	0.0	0.00106*

Geometric parameters (Å, º)

D9b—Ow9	0.979 (4)	D3a—Ow8^v	1.880 (6)
D9b—O4	1.792 (6)	D3a—Ow3	0.978 (4)
D9a—Ow9	0.972 (4)	D2b—Ow9^iv	1.834 (6)
D9a—O3^i	1.859 (6)	D2b—Ow2	0.974 (4)
D8b—Ow8	0.985 (4)	D2a—Ow5^vi	2.051 (6)
D8b—O1	1.942 (6)	D2a—Ow2	0.971 (4)
D8a—Ow8	0.992 (4)	D1b—Ow9^vii	1.811 (6)
D8a—O4^ii	1.882 (6)	D1b—Ow1	0.970 (4)
D7b—Ow7	0.988 (4)	D1a—O2^viii	1.761 (6)
D7b—O2	1.750 (6)	D1a—Ow1	0.968 (4)
D7a—Ow7	0.989 (4)	S1—O4	1.478 (4)
D7a—O3^i	1.893 (6)	S1—O3	1.491 (4)
D6b—Ow6	0.957 (4)	S1—O2	1.475 (4)
D6b—O3^i	1.774 (5)	S1—O1	1.479 (4)
D6a—Ow8	1.771 (6)	Mg2—Ow6	2.051 (4)
D6a—Ow6	0.982 (4)	Mg2—Ow6^vi	2.051 (4)
D5b—Ow5	0.971 (4)	Mg2—Ow5	2.108 (5)
D5b—O4^i	1.753 (5)	Mg2—Ow5^vi	2.108 (5)
D5a—Ow7^iii	1.695 (6)	Mg2—Ow4	2.034 (4)
D5a—Ow5	0.980 (4)	Mg2—Ow4^vi	2.034 (4)
D4b—Ow4	0.976 (4)	Mg1—Ow3	2.056 (5)
D4b—O1^iv	1.822 (6)	Mg1—Ow3^ix	2.056 (5)
D4a—Ow7	1.784 (6)	Mg1—Ow2	2.075 (5)
D4a—Ow4	0.969 (4)	Mg1—Ow2^ix	2.075 (5)
D3b—O1	1.933 (6)	Mg1—Ow1	2.028 (5)
D3b—Ow3	0.964 (4)	Mg1—Ow1^ix	2.028 (5)
D9b—Ow9—D9a	102.1 (6)	O4—S1—O2	110.0 (3)
D8b—Ow8—D8a	103.7 (6)	O4—S1—O1	110.0 (3)
D7b—Ow7—D7a	103.8 (6)	O3—S1—O2	109.6 (3)
D6b—Ow6—D6a	106.1 (6)	O3—S1—O1	108.7 (3)
D5b—Ow5—D5a	107.1 (6)	O2—S1—O1	109.8 (3)
D5b—Ow5—Mg2	116.6 (5)	D3b—Ow3—D3a	102.2 (6)
D5a—Ow5—Mg2	109.9 (5)	D2b—Ow2—D2a	103.8 (6)
D4b—Ow4—D4a	106.5 (6)	D1b—Ow1—D1a	106.5 (6)
Ow6—Mg2—Ow6^vi	180.0	Ow3—Mg1—Ow3^ix	180.0
Ow6—Mg2—Ow5	90.57 (18)	Ow3—Mg1—Ow2	91.40 (18)
Ow6—Mg2—Ow5^vi	89.43 (18)	Ow3—Mg1—Ow2^ix	88.60 (18)
Ow6—Mg2—Ow4	90.21 (18)	Ow3—Mg1—Ow1	92.27 (18)
Ow6—Mg2—Ow4^vi	89.79 (18)	Ow3—Mg1—Ow1^ix	87.73 (18)
Ow6^vi—Mg2—Ow5	89.43 (18)	Ow3^ix—Mg1—Ow2	88.60 (18)
Ow6^vi—Mg2—Ow5^vi	90.57 (18)	Ow3^ix—Mg1—Ow2^ix	91.40 (18)
Ow6^vi—Mg2—Ow4	89.79 (18)	Ow3^ix—Mg1—Ow1	87.73 (18)
Ow6^vi—Mg2—Ow4^vi	90.21 (18)	Ow3^ix—Mg1—Ow1^ix	92.27 (18)
Ow5—Mg2—Ow5^vi	180.0	Ow2—Mg1—Ow2^ix	180.0
Ow5—Mg2—Ow4	89.9 (2)	Ow2—Mg1—Ow1	90.33 (19)
Ow5—Mg2—Ow4^vi	90.1 (2)	Ow2—Mg1—Ow1^ix	89.67 (19)
Ow5^vi—Mg2—Ow4	90.1 (2)	Ow2^ix—Mg1—Ow1	89.67 (19)
Ow5^vi—Mg2—Ow4^vi	89.9 (2)	Ow2^ix—Mg1—Ow1^ix	90.33 (19)
Ow4—Mg2—Ow4^vi	180.0	Ow1—Mg1—Ow1^ix	180.0
O4—S1—O3	108.8 (3)

Symmetry codes: (i) −x+1, y+1/2, −z+1/2; (ii) x−1, y, z; (iii) −x+1, −y+1, −z; (iv) x, −y+1/2, z−1/2; (v) −x, y−1/2, −z+1/2; (vi) −x, −y+1, −z; (vii) −x+1, y−1/2, −z+1/2; (viii) −x+1, −y, −z; (ix) −x, −y, −z.

(9_K_ice) ice Ih

Crystal data

(D2O)	V = 128.23 (1) Å^3
Mr = 20.03	Z = 4
Hexagonal, P63/mmc	Melting point: 277 K
Hall symbol: -P 6c 2c	T = 9 K
a = 4.496577 (8) Å	Particle morphology: fine powder
c = 7.32306 (2) Å	colourless

Fractional atomic coordinates and isotropic or equivalent isotropic displacement parameters (Ų)

	x	y	z	Uiso*/Ueq	Occ. (<1)
O1	0.33333	0.66667	0.06232 (4)	0.0037*
D1	0.33333	0.66667	0.19801 (7)	0.01715*	0.5
D2	0.45364 (4)	0.90718 (8)	0.01730 (5)	0.01931*	0.5