Characterisation, axial anisotropy, and formation conditions of celestine minerals from the Jabal Eghei (Nuqay) late Neogene – Pleistocene volcanic province, southeastern edge of the Sirt Basin, southern Libya: Constraints on the mineralogical geothermometer

Abstract

Five celestine crystals were sampled from the (palaeo)surface intervening between the late Miocene to Pleistocene basaltic sequences of the Jabal Eghei (Nuqay) volcanic province in southern Libya and then characterised by applying a combination of the SEM–WDS, ICP/OES, PXRD and IR methods. Colour variations and related minerogenetic frameworks were also investigated. Three samples have greenish-blue-to-blue colour (480.4–482.5 nm), whereas the other two samples have blue–green colour (cyan; 489.1–494.1 nm). The colour purity ranges from 1.36–7.16. Their composition is similar, end-member celestine, in which only 1.6–4.1 at.% of Sr²⁺ content was substituted by Pb²⁺ (0.7–0.9 at.%), Ba²⁺ (0.5–0.7 at.%) and Ca²⁺ (0.2–0.8 at.%). Three samples contained vacancies, from 1.0 to 1.9 at.%. The content of other chemical elements is minor. The resulting unit-cell parameters have the ranges: a₀ = 8.3578(9)–8.3705(6) Å; b₀ = 5.3510(5)–5.3568(4) Å; c₀ = 6.8683(7)–6.8767(2) Å and V₀ = 307.17(5)–308.34(4) ų. The PXRD and IR results are mainly in accordance with the SEM–WDS results, with a high level of correlation. However, a few discrepancies were found, producing several possible interpretations, the primary cause being a slight unit-cell axial anisotropy i.e. thermal expansion. As a consequence these results yield a new geothermometric tool that is based on the unit-cell axial anisotropy. The celestines investigated were formed during a Miocene intraplate volcanism with basaltic magmas, and associated brines lifted by the structural conduits (normal faults crosscutting the Sirt basin). The Sr-bearing fluids then poured into and over the faulted and fractured lagoon-type gypsum, anhydrite Eocene sediments. The celestine mineralisation formed within a ~368–430 K (~95–157°C) temperature range. The celestine formed at slightly elevated temperature and pressure conditions, close to the shallow subsurface environment (over 250 bars).

Overview of celestine-bearing environments

Celestine (SrSO₄; ideally containing 56.42% SrO and 43.58% SO₃) is a common sulfate mineral in the Earth's crust and an important raw material used frequently for the production of the critical element strontium (Palache et al., 1951; Deer et al., 2013; see Zhu et al., 2022), used for example in electronics, ceramics, pigments, cosmetics, papermaking and high-temperature solid lubricants (e.g. Li et al., 2008). Together with anglesite (PbSO₄) and baryte (BaSO₄), celestine belongs to the sulfates of the baryte group, MSO₄. These are isostructural minerals crystallising in the orthorhombic Pnma (Pbnm) space group (No. 62) and are characterised by the tetrahedral SO₄ anion and cations of considerably large size (M = Sr²⁺, Pb²⁺ and Ba²⁺), the latter falling into twelve-coordination with oxygen. The SO₄ tetrahedra and MO₁₂ polyhedral share their edges. Notably, the calcium sulfate anhydrite (CaSO₄) is not in the baryte group as it differs structurally (space group Cmcm, No. 63) from celestine, primarily because of the relatively smaller sized Ca²⁺ ion (Hanor, 2000). Twinning in the celestine could be of different structure types and crystals may contain many defects, at least if they are grown from gels (Patel and Bhat, 1972).

Available data on naturally occurring celestine samples show that the majority have compositions close to that of the end-members, with Ba and Ca generally present at only a few at.%. Globally, reported celestine occurrences are dominantly of sedimentary origin and related to evaporitic processes (Glazner, 1988; Tekin, 2001; Álvaro et al., 2006; Ariza-Rodríguez et al., 2022). These celestine-bearing evaporitic basins, of different age ranges, are frequently affected by hydrothermal episodes (e.g. Dai et al., 2014; Fristad et al., 2017; Dhote et al., 2021). Alternatively, celestine minerals can occur as fracture and rock-cavity fillings (including in caves) precipitated by migrating strontium-bearing groundwater or basinal brines (‘celestitisation’; e.g. Hou et al., 2007). These celestine-bearing brines pour into carbonate rocks, concretions (e.g. Reolid et al., 2019; Allouche et al., 2023) and nodules (including agate gemstones; Götze et al., 2020). In addition, celestine occurrences can be as hydrothermal veins associated with mafic volcanism (e.g. Baskina et al., 1978; Brigatti et al., 1997; Radivojević et al., 2015; Garcia et al., 2016; Tortelli et al., 2022) or basaltic dykes (Anenburg et al., 2014). Occurrences of younger Quaternary celestines comparable to those investigated in southern Libya have been recorded only rarely (e.g. McFadden et al., 1987; Garcia et al., 2016).

Celestines are associated commonly with gypsum, anhydrite, native sulphur (evaporites), strontianite, calcite, dolomite, anhydrite, gypsum and fluorite (cavities in carbonate rocks); and also occur jointly with analcime, natrolite, hydroxyapophyllite and celadonite (mafic volcanoclastics). The presence of celestine can be the mark of a redox environment (Hadji et al., 2019); however, anhydrite, celestine and analcime are rarely observed within borate deposits (Helvaci, 1995). Nakamura et al. (2010) describe an occurrence from the Bohemian massif, Czech Republic of eclogite as a carrier of Sr-sulfate celestine.

Celestine colour can range from colourless to blue or have white, yellow, orange, reddish, greenish, or brownish colours (e.g. brownish celestine from a shear zone in Egypt; El-Sheriff, 2013). The blue colour is typically distributed unequally within the mineral growth zones or face loci, and in some cases, the blue colour could have been produced by irradiation (Palache et al., 1951; Hanor, 2000; Deer et al., 2013).

Several celestine minerals were collected during a geological field-mapping campaign in the central and southern Libyan desert (Figs 1 and 2, Supplementary Fig. S1). The field-mapping campaign lasted from 2005 up to 2010, carried out in cooperation between the Industrial Research Centre (Tripoli, Libya) and the Geological Survey of Serbia (former Geoinstitute, Belgrade, Serbia; Supplementary Fig. S1b; Toljić and Turki, 2007a, 2007b; Toljić et al., 2007; Rundić and Dalub, 2007; Vasić and Sheriff, 2007; Toljić and Abu Agrab, 2014, 2016). The field-mapping campaign encompassed the areas in and around two magmatic provinces: (1) the central Libyan Al Haruj volcanic system (Fig. 1e and Supplementary Fig. S1c) that is interconnected spatially with a southern-positioned Neogene to Quaternary segment of the famous Sirt basin; and (2) the Tibesti–Sirt arch or Jebel Eghei area (Rusk 2001; Cvetković et al., 2010, 2022; Kovačević et al., 2013; Radivojević et al., 2015; Elshaafi and Gudmundsson, 2016, 2017). Spatially, the Al Haruj province of Central Libya is in contact with the investigated southern-positioned Jabal Eghei Volcanic Province (Cvetković et al., 2010; Radivojević et al., 2015; Elshaafi and Gudmundsson, 2017; Fig. 1c,d; Fig. 2). The occurrence of the celestine minerals in both Libyan volcanic provinces indicates the presence of the two mineral modes: (1) ‘sedimentary celestine’, typically embodied in abundant surfacing druses or clusters occurring either in the sedimentary formations or within a weathered surface crust (palaeosoil) (Fig. 2b,c; Toljić and Turki, 2007; Vasić and Sheriff, 2007; Abdusahmin, 2020), and (2) ‘magmatic celestine’, generally surfacing as fracture-related mineral infill (Radivojević et al., 2015).

Figure 1. (a) Geographical relief map of the wider investigated central-southern Libya. (b) The sketch map showing the distribution of the four main volcanic provinces in Libya and principal faults delineating the celestine-bearing Sirt Basin and Jabal Nuqau Volcanic Province (modified after Elshaafi and Gudmundsson, 2016). The Al Haruj Volcanic Province (#3) is situated on the southwestern margin of the Sirt Basin, central Libya. To the south is the Jabal Eghei or Jabal Nuqau Volcanic Province (#4). Other provinces are (#1) Gharayan Volcanic Province, and (#2) As Sawda Volcanic Province. (c) Panoramic view on Jebel Eghei rhyolite volcanic complex (photo taken from Kovačević et al., 2013; modified). (d) Basalts from Jebel Eghei area (Abdusahmin, 2020, modified). (e) Basalts from Al Haruj area (photos taken from Vasić and Sheriff, 2007; modified).

Figure 2. (a) Synthetic stratigraphic column constructed after Gumati (2021), Rundić et al. (2012) and Radivojević et al. (2015); (b) basalt fragments and sedimentary formations (from Abdusahmin, 2020, modified); (c) occurrences of celestine crystals at the Al Gata member, formation Wadi Thámat (Coordinates: Lat. 23^o53’17’’; Long. 19^o02’31’’). Photo by J. Kovačević; (d) selected 1–5 celestine samples used in the study. Photos by P. Tančić.

The following study of the well-preserved celestine minerals collected from the Jabal Eghei Volcanic Province (Fig. 2c,d and Supplementary Fig. S1d,e) provides constraints on the origin of the selected surface-exposed very young celestine crystals (from druses or clusters). The Jabal Eghei celestine minerals were investigated with scanning electron microscopy using wavelength dispersive spectroscopy (SEM–WDS), inductively coupled plasma optical emission spectroscopy (ICP/OES), powder X-ray diffraction (PXRD) and infrared spectroscopic (IR) methods. The objective of this study was to resolve their composition and further characterise the unit-cell parameters and vibrational spectra. In addition to the methods applied, we establish constraints on celestine colour variations, including colour purities. Few publications exist focusing on the origin of the North African (e.g. Anenburg et al., 2014) and in particular central-southern Libyan basalts and associated celestine occurrences (e.g. Elshaafi and Gudmundsson, 2017; Radivojević et al., 2015; Cvetković et al., 2022); hence our interpretation of the Jabal Eghei celestine provides new insight into the axial anisotropy of the celestine crystals investigated and their growth. The regional synthesis examines the origin of the Jabal Eghei celestine minerals, yielding its formation temperature, and a set of other subsurface conditions under which the investigated crystals were formed. Further, this investigation shows that the origin of the Jabal Eghei celestine is similar to that of the Al Haruj basalt area (Fig. 1b).

Regional geological setting: Al Haruj and Jabal Eghei (Nuqay) areas

There are four principal volcanic provinces in the Central Sahara Desert, where a Libyan segment stretches from the Mediterranean coast towards the southern borderland of Libya: (1) Gharayan; (2) As Sawda; (3) Al Haruj Al Aswad; and (4) Jabal Nuqai or Jabal Eghei (Fig. 1b). The youngest are the central Saharan intra-plate mainly basalt volcanic provinces referred to as the Al Haruj and the Jabal Nuqay or Jabal Eghei (central and southern Libya; Fig. 1b and Supplementary Fig. S1c,d; Elshaafi and Gudmundsson 2016; Radivojević et al., 2015). The extruded large magmatic provinces and their surroundings host a number of occurrences of near-surface celestine-bearing systems (Rundić and Dalub, 2007; Toljić and Turki, 2007a, 2007b; Toljić et al., 2007; Vasić and Sheriff, 2007; Toljić and Abu Agrab, 2014, 2016; Abdusahmin, 2020). In the field, the celestine mineral clusters are associated either with sedimentary rocks (e.g. gypsum and/or anhydrite; Rundić and Dalub, 2007; Toljić and Turki, 2007b; Vasić and Sheriff, 2007; Rundić et al., 2012) or are observed directly on the basaltic rock assemblages (Toljić et al., 2007; Cvetković et al., 2010; Radivojević et al., 2015; Figs 1c,d,e, 2a,b,c). Both provinces are accommodated along the central-southern realm of the NW–SE striking Sirt basin, though the Jabal Eghei province reaches an uncovered segment of the African basement referred to as the Tibesti massif (Toljić and Turki, 2007; Le Heron et al., 2013; Radivojević et al., 2015; Abdusahmin, 2020; Fig. 1b and Supplementary Fig. S1c,d,e).

Sirt basin and Al Haruj magmatic-sedimentary province, the northern celestine-bearing assemblage

The northern branch of the Cenozoic tectono-magmatic system corresponds to the top-sealing Tertiary segment of the famous oil and gas-bearing Sirt Basin (Rusk 2001; Abouessa et al., 2012; Fig. 1b and Supplementary Fig. S1c). The underlying Lower Cretaceous Tethyan rift system (Gumati and Nairn, 1991; van der Meer and Cloetingh, 1993; Hallett 2002; Abadi et al., 2008; Capitanio et al., 2009) triggered the Late Cretaceous regional faulting and subsidence contributing to (1) the elevation of magma into the upper crust and episodic infiltration of hydrothermal fluids (e.g. Toljić and Turki 2007a; Khalifa and Morad 2012; Elshaafi and Gudmundsson, 2016, 2017); of which subsidence produced a 4.5 km thick overburden (El-Khatri et al., 2015) allowing deposition of (2) a number of Palaeogene formations that are widespread across the Sirt basin (Toljić and Turki, 2007b; Rundić and Dalub, 2007; Vasić and Sheriff, 2007; Rundić et al., 2012) (onset of the subsidence stage at Cenomanian–Turonian; Rusk 2001). In terms of the oil exploration ‘language’, the syn-rift sedimentary loading of the newly developed accommodation space lasted from the Cretaceous to the Eocene, whereas the post-rift subsidence and deposition continued during the Oligocene and Miocene. The oldest exposed rocks are Eocene Ypresian–Lutetian, Lutetian–Bartonian and late Middle Eocene stages which were formed in shallow sea, lagoon-type environments, whereas the two overlying sequences are of the Oligocene and Lower to Middle Miocene age (Dur-at-Talah and Maradah formations; Rundić and Dalub, 2007; Toljić and Turki, 2007; Vasić and Sheriff, 2007; Abouessa et al., 2012; see Rundić et al., 2012, for details). These Eocene celestine sedimentary carriers (druses) are mainly formed of subhorizontal layers of lagoon-type mixed clastics, limestone, evaporates and evaporitic limestone sequences, stretching from the Al Haruj Basaltic complex further into the Jebel Eghei area (Abdusahmin, 2020; Fig. 1b; 2a,b,c and Supplementary Fig. S1c,d,e). The fibrous gypsum deposits fill out the bedding planes and vertical fractures (post-depositional feature or ‘secondary gypsum’; Abouessa et al., 2012).

Jebel Eghei magmatic-sedimentary province, the southern celestine-bearing assemblage

The Jabal Eghei basaltic volcanic province is situated in southern Libya (Fig. 1b and Supplementary Fig. S1d,e), representing an outcrop fragment of the north-western edge of the Tibesti Mts. The basalts extruded over a superposed metamorphic basement and sedimentary rocks (including the Tertiary formations similar to those in the Sirt basin; Toljić and Abu Agrab, 2014; Radivojević et al., 2015; Abdusahmin, 2020). The older part of this metamorphic–sedimentary assemblage is represented by the continental Neoproterozoic to Cambrian Mourizide Formation (Le Heron et al., 2013). A younger part of the sedimentary cover belongs to the Sirt basin, its southernmost periphery. It is composed of Palaeocene–Oligocene sediments, which include the Jabal Eghei basalts (Fig. 2a). In addition to the Cretaceous rifting, the extensional faulting across the Sirt basin probably contributed to the continental crustal opening and extrusion of the basalts across the southern sector of the Jabal Eghei province (Radivojević et al., 2015; Elshaafi and Gudmundsson, 2017; Cvetković et al., 2022; Abdusahmin, 2020).

The recent systematic K/Ar radiometric age determination of the Jabal Eghei volcanism pinpoints it to the Middle Miocene (~16 Ma) to mid-Pleistocene activity (<1 Ma; Radivojević et al., 2015). The three different volcanic phases include the oldest phase in the Middle to Late Miocene (basalts show K/Ar radiometric age ranges from 16.1±2.9 to 7.9±2.3 Ma); a second phase from 7.19±0.36 to 4.32±0.35 Ma; and the youngest extrusion stage at 3.1±1.1 to 0.97±0.68 Ma. The celestine samples investigated are mapped on top of the first volcanic extrusion phase. The samples were located in the weathered crust (Fig. 2b), intervening between the second and the third basalt outflow magmatic episodes (Lower Pliocene to Middle Pleistocene; Supplementary Fig. S1e). Considering the extensive celestine presence in the Eocene formations (Toljić and Turki, 2007b; Vasić and Sheriff, 2007; Abdusahmin, 2020) and the intense multi-phased basaltic volcanic extrusions (Cvetković et al., 2010; Radivojević et al., 2015; Elshaafi and Gudmundsson, 2017), the exact timing of the emplacement of Sr-rich brine remained uncertain.

Materials and methods

Materials

Large crystals (up to 12 cm long; Fig. 2c) were hand-picked during the field-mapping campaign from druses or clusters exposed on the surface. Among the crystals sampled, five well-formed euhedral to subhedral samples were chosen. The choice was based on their macroscopically visible different characteristics, e.g. colour, lustre, transparency, shape, size and cleavage (Fig. 2d).

Methods

The samples analysed (#1–5) show no radioactivity (measured in situ by the Saphymo-SRAT S.P.P.2 NF scintillometer, manufactured by Saphymo-PHY, Massy, France).

The electron microprobe analytical (EMPA) data and back-scattered electron (BSE) images were obtained on a polished sample surface, using a JEOL JSM-6610LV scanning electron microscope equipped with a wavelength-dispersive spectrometer (Oxford Instruments Wave 700). Operating conditions were at 30 kV accelerating potential, 10 nA beam current and 2 μm spot size dimension. Standards used for this analytical method were natural and synthetic mineral standards, including celestine (Sr), wollastonite (Ca), baryte (Ba and S) and galenite (Pb).

The contents of selected trace elements (Ag, As, Al, Be, Bi, Cd, Ce, Co, Cr, Cu, Fe, Ga, Ge, Hf, Hg, Mg, Mn, Mo, Nb, Ni, P, Pb, Sb, Sc, Se, Sm, Sn, Ti, Tl, W, Y, Zn and Zr) and rare earth elements (lanthanides and Th) were obtained using ICP/OES (Spectro Blue, Germany). Measurements were performed after total decomposition of the samples by fusion with NaKCO₃ at 1000°C, and dissolution with HCl.

Powder XRD was conducted on a Rigaku Smartlab X-ray Diffractometer, in θ-θ geometry (the sample in a horizontal position). We used parafocusing Bragg-Brentano geometry using a D/teX Ultra 250 strip detector, in the 1D standard mode with CuKα_1,2 radiation source (40 kV and 30 mA). The PXRD patterns were collected in the 5–90°2θ range, with a step of 0.01° and a data collection speed of 6.15°/min. A low background single-crystal silicon sample holder was used to minimise the background.

Unit-cell parameters were calculated and refined using the Le Bail (1988) method. The sequence of the full-matrix least-squares fittings was refined by changing the unit-cell parameters (a ₀, b ₀ and c ₀), 100 parameters for the background description, Eta (p-v) parameter, parameters for the description of the halfwidths (U, V and W) and zero point; until the refinements approached convergence.

The infrared (IR) spectroscopic analyses were performed using a Perkin Elmer 597 spectrometer, and the KBr pellet method for samples (2 mg of sample and 200 mg of KBr). Visible and near-infrared absorption spectrums were measured in the region 1600–200 cm⁻¹.

The dominant wavenumbers of colour for the raw celestine samples were measured with a diffuse reflectance apparatus (400–700 nm) by using a CCS200 spectrometer (Thorlabs). This spectrometer is equipped with an optical fibre wire and is calculated according to the Commission Internationale de l'Eclairage (1932).

The organic matter contents were detected by loss of ignition at 600°C.

Results

Wavelength dispersive spectroscopy

The sample-based SEM-WDS analyses yielded homogeneous composition results, revealing no zoning in Sr, Ca, Ba, Pb or S within the samples and various mineral phases (Table 1 and Supplementary Fig. S2). The results show similar contents for all the celestine samples analysed. The major SrO and SO₃ contents have the range of 54.33–55.21 (Δ_SrO = 0.88) and 43.23–43.93 (Δ_SO3 = 0.70) wt.%, respectively. Among minor oxides, PbO, BaO and CaO are in the range of 0.87–1.04 (Δ_PbO = 0.17), 0.42–0.62 (Δ_BaO = 0.20) and 0.07–0.24 (Δ_CaO = 0.17) wt.%, respectively. These concur with the major celestine minerals having minor anglesite, baryte and anhydrite substituted components.

Table 1. Compositions from WDS analyses of samples 1–5 (in wt.%). Atoms per formula units (hereinafter apfu; in at.%) are calculated on a 4 oxygen anions basis.

The atoms per formula unit (apfu) calculated based on the 4 oxygen anions (Table 1) show that samples 1 and 4 have ideal (or almost ideal) stoichiometry, between M and S ions; including the crystal-chemical formulas of (Sr_0.984Pb_0.009Ba_0.005Ca_0.002)_Σ1.000SO₄ and (Sr_0.978Pb_0.008Ca_0.008Ba_0.007)_Σ1.001SO₄, respectively. In contrast, samples 2, 3 and 5 show a minor non-stoichiometry, having a small deficit of M cations of ca. 1, 1 and 2 at.%, respectively. To account for the fact that the sulfur content shows a slight excess from the value of 1, we propose that such a discrepancy was induced by a possible presence in the total sum (hereinafter X) of H⁺ and/or (H₃O)⁺ and/or OH^– (from H₂O) and/or vacancies. These contents would not be detected using WDS. Accordingly, the crystal-chemical formulas of samples 2, 3 and 5 could be presented as (Sr_0.970X _0.010Pb_0.008Ba_0.007Ca_0.005)_Σ1.000S_1.003O₄, (Sr_0.973X _0.010Pb_0.007Ba_0.005Ca_0.005)_Σ1.000S_1.003O₄ and (Sr_0.959X _0.019Pb_0.008Ba_0.007Ca_0.007)_Σ1.000S_1.006O₄, respectively. Therefore, all the samples studied are mutually similar and are celestine near-end-members, in which only 1.6–4.1 at.% of Sr²⁺ content (Δ_Sr = 2.5 at.%) is substituted by Pb²⁺ (0.7–0.9 at.%), Ba²⁺ (0.5–0.7 at.%), Ca²⁺ (0.2–0.8 at.%) and the X (1.0–1.9 at.%; in three samples) component. The apfu content of Sr, between the samples, increases in the order 5<2<3<4<1, whereas the ratios of Σ(Pb+Ba)/Σ(Ca+X) increase according to 5<3<2<4<1.

By using the apfu 12-coordinated M ²⁺ cations, i.e. Ca²⁺, Sr²⁺, Pb²⁺ and Ba²⁺ derived from the WDS results (Table 1) and their ionic radiuses (r _M) of 1.34Å, 1.44Å, 1.49Å and 1.61Å, respectively (Shannon, 1976), we have calculated the theoretical ionic radiuses of these cations (Table 2). The results show that the calculated occupancies of the twelve-coordination site are as expected, being also mutually quite similar, further attesting to the previously suggested similarities among the samples analysed (Δ_occ = 1.99 at.%). This further includes the deficit of the M cations of ~1, 1 and 2 at. % for samples 2, 3 and 5, respectively. In this case, the increase of the calculated occupancies is in a slightly different order of 5<3<2<1<4.

Table 2. Calculated theoretical ionic radiuses (in Å) of the M cations, and calculated occupancies of the twelve-coordination site (in at.%) at the basis of the determined apfu (Table 1).

*Occupancies were calculated by the Σr_M/1.44 Å ratios.

ICP/OES

The content of the detected trace elements from ICP/OES is presented in Table 3. All other analysed elements are below a detection limit of 0.011–0.500 mg kg^–1, and consequently, have been disregarded.

Table 3. Determined contents of the trace elements (in mg kg^–1).

Powder X-ray diffraction

The resulting PXRD patterns [Supplementary Fig. S3 (column I) and Table S1] demonstrate that the observed inter-planar spacings (d _obs) are very similar (all samples), and are in line with the celestine ICDD–PDF 89-0953 and 05-0593 reference standards (from the International Centre for Diffraction Data, see Supplementary information for data references and sources). As no other minerals are detected, these samples should be exclusively monomineral phases and the data should corroborate the SEM-WDS results. However, the observed intensities (I _obs) show a distinct difference relative to the applied standards. For example, the strongest standard reflection (hkl 211) is only 25% for sample 1 (with the relative intensity ratio-RIR applied). In addition, a comparison of the reflection patterns for different Miller's hkl indices shows a significant variety between the samples analysed (Supplementary Fig. S4). Namely, samples 3 and 5 show distinctively higher RIRs that are directed towards the c ₀ axis, i.e. the [00l] direction (hkl: 002, 004 and 006; Supplementary Fig. S4a,g,j). The rest of the reflections, i.e. samples 1, 2 and 4 have higher RIRs (Supplementary Table S1). Samples 2 and 4 show a less prominent trend inclining towards the b ₀ axis, i.e. [0k0] direction (hkl: 020 and 040; Supplementary Figs S4e right and S4i left), as well as towards hkl 323 (Supplementary Fig. S4h). Such differences could be induced by the preferred orientation, accounting for (1) a variety in the cleavage properties (i.e. {001} perfect, {210} good and {010} poor; Palache et al., 1951; Deer et al., 2013). Differences can also be caused by (2) epitaxial growth of the celestine on anhydrite (Forjanes et al., 2020b). We also took into consideration the significant presence of calcium in all specimens (Table 1). In addition, differences could also be produced by the structural variation(s), which will be discussed later.

By using the equation of Goldish (1989), we have calculated the content of the substituted SrSO₄ mole fractions (Table 4). Taking into account that this equation is valid for the baryte–celestine solid solution (and hence obtained results should be treated cautiously), we presumed tentatively that from 0.2 (samples 1 and 5) to 2.4 (sample 4) mole fractions of SrSO₄ were replaced by BaSO₄. However, because the Pb²⁺ cation (Table 1) has a slightly larger ionic ratio size than Sr²⁺, yet smaller than Ba²⁺ (Shannon, 1976), these molar substitutions of SrSO₄ should exhibit a slightly larger presence. On the other hand, from 0.2 (sample 3) to 0.6 mol.% (sample 2) of SrSO₄ could be replaced by a smaller compound (such as CaSO₄, for example) or could be caused by vacancies at the M site. Accordingly, the maximal difference of the SrSO₄ content among the investigated samples is ~3 mol.%, having the substitutions increasing in the following order: 2<3<1=5<4.

Table 4. Calculated SrSO₄ mole fractions from the observed d ₍₂₁₁₎ inter-planar spacings (Table S1).

*Calculated from the equation by Goldish (1989): d ₍₂₁₁₎ = 3.1023–0.1317 × S, where S is the SrSO₄ mole fraction.

Unit-cell parameter calculations and comparisons

The unit-cell parameters (hereinafter UCPs) were calculated and refined by using the Le Bail profile fitting method [Supplementary Figs S3 (column II) and S5; and Table S2], and presented in Table 5. As expected from the observed inter-planar spacings (d _obs, Supplementary Table S1), all of the calculated UCPs are fairly similar and fit with the reference data standards. However, there are different relationships between the parameters of each sample for different characteristics. For example, the a ₀ and c ₀ axes lengths for samples 1–5 increase in order 2<3<5<1<4; however for the volume V ₀ it is 2<3<1≈5<4; whereas for the b ₀ axis, the order is 2<1<3≈5<4.

Table 5. Calculated unit-cell parameters of the investigated samples in the Pnma space group, compared to the celestine reference data standards; and also between each other.

*The numbers in parentheses in this and further Tables are the estimated standard deviations (hereinafter esd's) multiplied with Scor (Bérar and Lelann, 1991; Supplementary Table S2), and refer to the last significant number.

To examine these data and mutual correlations between the UCPs, we used MS Excel to calculate statistical parameters and the coefficients of the regression-R² (hereinafter correlations-C). These correlations, reported in our previous studies of different materials (Cvetković and Tančić, 2019; Tančić et al., 2021; Maksimović et al., 2023) are proven useful accessory tools, which allow the comparison between various parameters. According to the results, a polynomial regression between the a ₀, c ₀ and V ₀ parameters gives slightly better R ² coefficients than linear regression (Supplementary Fig. S6), suggesting deviations from Vegard's law, which is in concordance with the results of Goldish (1989) and Antao (2012). A very high positive correlation C(1) is present for the axes/volumes: a ₀ vs c ₀; a ₀ vs V ₀; and c ₀ vs V ₀ of 0.999; 0.985; and 0.990, respectively. The variations of the b ₀ axis length versus a ₀, c ₀ and V ₀ parameters (Supplementary Fig. S7) display slightly lower, yet likewise very high positive correlations, C(2) of 0.880, 0.951, and 0.956, respectively. The correlations for the various UCPs of the samples investigated (Supplementary Figs S6 and S7) were further compared with the celestine (ICDD-PDF 89-0953), anglesite (36-1461) and baryte (24-1035) reference standards in Supplementary Fig. S8. A lower C(2) correlation (Supplementary Fig. S8a,c,e) for the parameters a ₀ vs b ₀, c ₀ vs b ₀, and V ₀ vs b ₀ may be significant, accounting for the relationship with the larger angle declinations documented between the celestine–anglesite and the celestine–baryte solid-solution joins. Such differently angled declinations were induced by a different expansion behaviour of the UCPs between the Sr²⁺–Pb²⁺, Sr²⁺–Ba²⁺ and Pb²⁺–Ba²⁺ substitutions, i.e. between the various solid-solution series. Specifically, the calculated differences in this study (Supplementary Table S3) show that the highest expansion behaviour is for the a ₀ and c ₀ axes, and the volume V ₀. On the contrary, different substitutions do not significantly affect the b ₀ axis.

Infrared study

Major vibrational areas exist in the regions between 435–475, 602–650, 995–1005, 1080–1192 and 1200–1248 cm^–1 (Table 6 and Supplementary Fig. S9). Two broad bands were observed across all samples in the region between 435 and 475 cm^–1 and are most pronounced in samples 1 and 2. The vibrational band at 613–614 cm^–1 is highly pronounced in samples 1 and 2, whereas its position in sample 3 is slightly shifted towards 608 cm^–1. The most prominent shift is observed in samples 4 and 5, having the band peak at 602 and 604 cm^–1, respectively. A smaller doublet is visible between 610 and 622 cm^–1. This doublet is observed at 622 cm^–1 for samples 1 and 2, it is shifted to 612 cm^–1 in sample 3, whereas in samples 4 and 5, it is observed at 609 and 610 cm^–1, respectively. The vibrational features in the region 995–1005 cm^–1 are also more pronounced in samples 1 and 2. The region between 1080 and 1248 cm^–1 is a complex group of bands, wherein some of them have the form of a doublet band.

Table 6. The summarised IR results (in cm^-1) for the investigated samples together with the literature data.

References: 1* Ross (1972); 1a* Ross (1974); 2* Guirguis (1987); 3* Omori (1968); 4* Gadsden (1975), 5* Moenke (1962); 6* Kloprogge et al. (2001); and 7* Lane (2007).

Minerals from the baryte group are isostructural, resulting in vibrational spectra that are mutually similar, though the presence of different cations in the structures causes small offsets in the bands (Lane, 2007). According to literature data (Table 6), the internal vibrational features present in solid-state sulfates appear at ca. 1050–1250 (ν3), ~1000 (ν1), ~500–700 (ν4) and ~400–500 (ν2) cm^–1 as a result of the asymmetric and symmetric stretching and bending of the SO₄ anion (Herzberg, 1945; Nakamoto, 1986; Vassallo and Finnie, 1992; Bishop and Murad, 2005; Lane, 2007), including the lattice vibrations at ca. 550 cm^–1 and lower (Serna et al., 1986; Clark, 1999; Lane, 2007). The ν1 band is known to be significantly stronger and better pronounced than the ν2 band, which is generally not observable in the infrared spectra of sulfates (Hezel and Ross, 1966; Lane, 2007). Lane (2007) reports that a single ν2 band at 471 cm^–1 is visible in the spectrum of celestine, whereas Ross (1972, 1974) observed this band at ca. 491 cm^–1. According to Lane (2007), the spectral features for Sr as a major element in the celestine are evident at 1238 and 1138 cm^–1 for ν3; the ν4 band is evident at 648 and 614 cm^–1; whereas at 991 cm^–1, the ν1 band is observable. Moenke (1962) reports that the deformations of the SO₄ tetrahedra in the optically uniaxial sulfates occur as a doublet on the 620 and 1150 cm^–1 band.

In the baryte-group minerals, the isomorphous substitution of Sr and Ba is common, in particular, between celestine (SrSO₄) and baryte (BaSO₄), in addition, celestine may contain a limited amount of Ca²⁺ (Bernstein, 1979). As the Ca²⁺ ion is almost the same size as Na⁺, we also take into consideration that substitution between these two ions may also have been occurring. However, the replacement of Sr²⁺ is more likely by K⁺ ions than Na⁺, as K⁺ has a similar size to Ba²⁺. It has been hypothesised that any charge imbalance from this substitution could be compensated with additional Na⁺ or possibly H⁺ substitutions (Bernstein, 1979). The WDS data of the samples investigated (Table 1) show that the contents of K₂O and Na₂O are below the detection limit, indicating that these elements are not present in the substitutional positions in the structure. However, PbO, BaO and CaO are present, up to 1 wt.%. The vibrational features of Ba have a slight offset from the Sr vibration bands and are documented at 1220, 1128 (ν3); 641, 611 (ν4); and 981 cm^–1 (ν1) (Lane, 2007). A similar shift of the ν3 fundamental band is observed within the Ca-sulfate phases occurring at a lower-frequency (lower wavenumber) position (Lane, 2007). The spectral features of ν3, ν4 and ν1, are further systematically offset to smaller wavenumbers for the sulfates containing Pb (ν3 = 1183, 1054; ν4 = 632, 598; ν1 = 960 cm^–1) (Lane, 2007). The results compare well with the other applied methods, in particular the SEM–WDS and PXRD, indicating that the samples are of a monomineral composition.

Colour study

The spectrometer results gave various dominant wavelengths for samples 1–5 (Table 7). Namely, samples 1, 2 and 4 exhibit several nuances, ranging between 480 and 482 nm corresponding with a greenish-blue to blue colour (i.e. between 492 and 455 nm). Due to the observed colour purity, these shades are close to pure white (Supplementary Fig. S10). However, samples 3 and 5 exhibit a dominant wavenumber in the range corresponding to blue–green (cyan) colour. Sample 3 has the lightest shade of the colour (Table 7) and is nearest to pure white. Sample 5 has the most intense colour shade, and a purity value farthest from pure white, compared to the other celestine samples (see Supplementary Fig. S10).

Table 7. Dominant wavelength (in nm) and the colour purity of the investigated samples. The contents of the organic matter (in wt.%) are also presented.

The colour centres of the functional groups, such as SO₄^–, SO₃^–, SO₂^– and O^– show absorption in the visible light area, which effectively leads to the production of the blue colour of the samples (Bernstein, 1979). Nevertheless, this blue colour can also be influenced by an eventual presence of organic inclusions (Stromeyer, 1821) or inclusions of additional minerals such as vivianite, colloidal gold or sulfur (Wittstein, 1856; Doelter, 1915; Friend and Allchin, 1939, 1940). The blue colour can also be affected by the charge transfer between the elements in the structure, and/or their substitutions (Schulman and Compton, 1962; Bershov and Marfunin, 1967; Isetti, 1970; Bernstein, 1979). The chemical analyses of the celestine samples (Table 1) show that the content of gold is under the detection limits, fitting with the Bernstein (1979) report, which explains that the blue colour, in this case, cannot be influenced by the colloidal gold-theory, or with the content of Fe or Na in the samples [which are also disregarded (Table 3), or under the detection limit]. We underline, that no additional inclusions of minerals were detected by this study, thus excluding this theory as well. The difference in the formation temperatures indicates that the colour is influenced by the colour centres; however, there is no uniform association between them. Except for the colour centres of the functional groups, there is a correlation between the highest contents of Ca and the organic matter, exclusively in sample 4 (Tables 1 and 7), leading towards a deepest blue shade.

Discussion

Disagreements between WDS and PXRD results

A few revealing discrepancies occur between the WDS and the PXRD results (Tables 1, 2, 4 and 5). If the samples are ordered by data value for each parameter, then samples 3 and 5 vary with sample 2 [i.e. 5<3<2<1<4 (calculated occupancies of the twelve-coordination site; Table 2) vs 2<3<1=5<4 (molar substitutions of SrSO₄; Table 4) and 2<3<5<1<4 (axes a ₀ and c ₀ lengths; Table 5)], whereas sample 1 changes position for just a few of the parameters and only with sample 4 {i.e. 5<2<3<4<1 (apfu contents of Sr; Table 1) and 5<3<2<4<1 [ratios Σ(Pb+Ba)/Σ(Ca+X); Table 1] vs 2<3<5<1<4 (axes a ₀ and c ₀ lengths; Table 5)}. Thus, there is a poor correlation C(3) of the ionic radius (Table 2) with a ₀, b ₀, c ₀ and V ₀ UCPs (Table 5), presented in Supplementary Fig. S11 and Table S4 (which lists all the R ² values). The possible reason for these disagreements could be discovered by assessing the presumed structural and/or compositional circumstances, such as: (i) poorly calculated UCPs; and/or (ii) the wrong apfu basis; and/or (iii) the presence of other previously undetected mineral phases; and/or (iv) the incomplete calculation of ionic radiuses; and/or (v) the option that various structural variations within the samples could take place. Despite seemingly perplexing or even (in)significant, we believe that the listing of these possible reasons could be of the highest importance and that it is important to investigate each, and avoid any eventual misleading interpretations. They are discussed in turn below.

(i) Unit-cell parameter calculations

This investigation uses the Le Bail (1988) method because it has been applied previously to the baryte group of minerals (Kuang et al., 2017; Li et al., 2018; Girard et al., 2019; Ye et al., 2019). As already shown, very good profile parameters and relatively low-reliability factors are achieved by this method (Supplementary Table S2), a satisfactory level of quality for the resulting difference plots [Supplementary Figs S3 (column II) and S5], and very high positive correlations C(1) and C(2) among the UCPs (Supplementary Figs S6 and S7). Such results corroborate that the investigated samples are celestines and their Pnma structural model. In addition, considerably low estimated standard deviations (Table 5), which are mostly better than the reference standards used, provide validation of the UCPs calculation and refinement. Therefore, we rejected this option as a possible reason for the aforementioned disagreements between the WDS and the PXRD results.

(ii) Alternative apfu calculation

We calculated the apfu by using 4 oxygen anions (Table 1), which led to the estimated theoretical ionic radiuses of the M cations (Table 2). Another option is to calculate apfu on the basis of 2 (ΣM+S) ions. For example, Antao (2012) uses this option to give a celestine-derived (Sr_0.966Fe_0.001Ba_0.001)_0.968(SO₄)_1.027 formula, which fits the calculated (Sr_0.953Fe_0.001Ba_0.001)_0.955S_1.016O₄. As expected, there are no significant differences for samples 1 and 4; whereas samples 2, 3 and 5, show a slightly higher content of Sr, S and O ions compared to the 4 oxygen basis (Supplementary Table S5). The results demonstrate that the samples have mutual similarity, fitting into the celestine near-end members, having only 1.6–3.5 at.% of Sr²⁺ content (Δ_Sr = 1.9 at.%), substituted by the same content of Pb²⁺ (0.7–0.9 at.%), Ba²⁺ (0.5–0.7 at.%) and Ca²⁺ (0.2–0.8 at.%). The variations of the apfu content of the Sr and Σ(Pb+Ba)/Σ(Ca+X) ratios, i.e. 5<2<3<4<1 and 5<3<2<4<1, respectively, remain identical with the 4 oxygen basis. The only difference is a slightly lower content of the X component (0.7–1.3 at.%; Δ_X = 0.3–0.6 at.%) which in turn, leads to slightly higher calculated ionic radiuses Σr_M for the aforementioned three samples, and their 5<3<2<4<1 variation (Supplementary Table S6). However, the resulting correlations [C(4); Supplementary Fig. S12 and Table S4] among these ionic radiuses variations with the a ₀, b ₀, c ₀ and V ₀ parameters (Table 5), are significantly lower R ² values than those of C(3). Such a low correlation allows rejection of this possibility as well.

(iii) Presence of other phases

Although the observed composition is rather simple (Table 1), the X component could hypothetically indicate the presence of other mineral phases, such as gypsum (CaSO₄·2H₂O), anhydrite (CaSO₄), native sulphur (S) and calcite (CaCO₃). It is well known that these mineral phases frequently occur in the paragenesis with celestine (Palache et al., 1951; Deer et al., 2013). In addition, there is the possibility of the presence of numerous other minerals, such as anglesite, baryte, strontianite (SrCO₃), aragonite (CaCO₃), cerussite (PbCO₃), hydrocerussite [Pb₃(CO₃)₂(OH)₂], witherite (BaCO₃) or a SrSO₄·½H₂O phase (Takahashi et al., 1993), making the interpretation even more complex. To check for the possible presence of anhydrite (in samples 1 and 4) and gypsum (in samples 2, 3 and 5; including other minerals with the X component), we implemented additional recalculations of the WDS analyses (Supplementary Table S7). These recalculations are based on the CaO component, the apfu deficit of the M cations; including the sulfur excess (Table 1), and the theoretical composition of the aforementioned possibly involved minerals.

According to the results, it could be speculated that there is a tentative minor content of (1) anhydrite in samples 1 (0.17 mol.%) and 4 (0.58 mol.%); (2) gypsum in the samples 2 (0.49 mol.%), 3 (0.46 mol.%), and 5 (0.68 mol.%); and (3) other minerals with the X component in the samples 2 (0.10 mol.%), 3 (0.30 mol.%) and 5 (0.60 mol.%). However, these contents are below the detection limit for PXRD of 1–2% so no peaks were found, though we closely inspected the 2θ angle regions, which are the strongest reflections for these phases (Supplementary Fig. S3). In addition, the WDS method could have failed to separate these other minerals from the celestine, because of the similarities in the composition and its limitations in detecting some elements (previously mentioned). The apfu based on 12-coordinated M ²⁺ cations (Supplementary Table S7) was recalculated by applying the theoretical ionic radiuses of these cations (Supplementary Table S8). These correlations [C(5); Supplementary Fig. S13] of such recalculated ionic radiuses versus the a ₀, b ₀, c ₀ and V ₀ parameters, despite having a considerably better fit than C(3) (Supplementary Table S4), are still quite far from the ideal case. Therefore, this option does not provide a complete explanation of the previously discussed discrepancies and their link with the PXRD results, thus having a tentative character.

In contrast, the presence of other minerals with the X component in the paragenesis with the celestine should be rejected, due to their lower correlation [C(6); Supplementary Fig. S14] in comparison with the previously considered two options [i.e. C(3) and C(5); Supplementary Table S4].

(iv) The incomplete calculation of ionic radiuses

According to the apfu deficit of the M cations in samples 2, 3 and 5 (Table 1), uncoupling it with other chemical elements that should be disregarded (Table 3), we are of the opinion that the theoretical ionic radiuses could need another recalculation. This time by adding the appropriate sulfur excess of 0.3, 0.3, including 0.6 at.% of the S⁶⁺ taken from the tetrahedral site (0.12 Å in coordination IV; Shannon, 1976), respectively. Two options are chosen (Supplementary Table S9): (a) monomineral celestines (Table 2), and (b) major celestines that have barely any anhydrite or gypsum (Supplementary Table S8). A similar correlation of these ionic radiuses with the UCPs is observed in these two cases [C(7) and C(8); Supplementary Figs S15 and S16, and Table S4]. However, neither the first nor the second option reach any improvement, in regards to data from C(3). Consequently, this correlation should also be rejected as a possible explanation for the aforementioned disagreements, thus leaving one remaining possibility to be investigated, below.

(v) The option that various structural variations within the samples could take place

Although the crystal structure refinements are beyond the scope of this paper, our previous studies of diverse mineral species and solid solutions, have demonstrated that the different UCPs, polyhedral distortions, site occupancy factors, bond lengths, bond angles and valence units, could be influenced either by: (1) different inserted cations into the mineral structure; and/or by (2) diverse conditions of their formation, in particular, the temperature and pressure (Tančić 2017, 2018; Tančić and Kremenović, 2022; Tančić et al., 2012, 2020, 2023).

For situation (1), variations of the UCPs due to different inserted cations in the mineral structure, we calculated presumed values by the multiplication of the Sr, Ca, Pb and Ba apfu with the corresponding UCPs of celestine, anhydrite, anglesite and baryte reference standards (ICDD–PDF's: 89-0953, 37-1496, 36-1461 and 24-1035, respectively). For selection of the apfu, we chose the two previously discussed options: (a) the celestines as monomineral phases (Table 1); and (b) the celestine sample data calculated without Ca, which hypothetically belongs to minor amounts of gypsum and anhydrite phases (Supplementary Table S7). The results (presented in Supplementary Table S10) indicate that the largest differences (Δ) between these presumed values and those of the PXRD method (Table 5) are valid for both options, occurring in samples 2, 3 and 5. Such a correlation is expected and suited to the already emphasised vacancies and the here-used calculations. Accordingly, the ratios increased by ca. 1, 1 and 2%, respectively. In both cases (as monomineral and without Ca), we observe excellent correlations of the calculated ionic radii (Table 2 and Supplementary Table S8) with the a ₀, b ₀, c ₀, and V ₀ parameters [C(9) and C(10); Supplementary Figs S17 and S18, and Table S4]. Because both options are almost identical, this could not be used as a reliable parameter for the already distinguished factors. On the other hand, these results further indicate that the presumed UCPs should be theoretically near the realistic chemical composition. However, as discussed earlier, such a scenario is not the case. Namely, the presumed unit-cell volumes of samples 2, 3 and 5 do not fit to any case, because they have significantly lower values, inclusive of the comparison with the anhydrite (ICDD–PDF: 37-1496; V = 305.60 ų) reference standard. Consequently, we observe a very poor correlation between the presumed (Supplementary Table S10) and determined (Table 5) a ₀, b ₀, c ₀ and V ₀ UCPs [C(11); Supplementary Fig. S19 and Table S4].

Thus, we further reconsider the impact of the existing vacancies in the celestine structure on their UCPs. The available reference data (Supplementary Table S11) show that there is an example with a significant vacancy of 3 at.% [Antao, 2012; already discussed within part (ii)]. This celestine example has more or less similar [ICDD–PDF's: 05-0593 and 89-0953; Table 5; Miyake et al. (1978)], or even larger UCPs, comparative to those with no vacancy [Hawthorne and Ferguson, 1975; Jacobsen et al., 1998; Ye et al., 2019]. Likewise, three selected samples showed identical average <M–O> distances of 2.827(1) Å (Hawthorne and Ferguson, 1975; Jacobsen et al., 1998; Antao, 2012). Accordingly, it appears that the existing vacancies in the celestine structure have no significant impact on their UCPs. The causes for such behaviour could be the geometrical changes of the SO₄ and MO₁₂ polyhedral, i.e. a decrease of the average <M–O> distance, the charge on the O ions that have more subordinate values, including the average <S–O> distance that is lengthier in this case; or vice versa (Antao, 2012). In addition, the tetrahedral distortions are a function of the geometry of the structure rather than the chemistry of the twelve-coordinated site (Hawthorne and Ferguson, 1975; Brigatti et al., 1997).

To check the afore-described option, another two variances discussed earlier (ii and iv) yet without the existing vacancies (included for samples 2, 3 and 5) are taken into consideration: (a) celestines as monomineral phases, with previously characterised apfu (Table 1), including the theoretical ionic radiuses of the M cations with occupancies of the twelve-coordination site (Table 2) that are further recalculated (Supplementary Tables S12 and S13), and (b) the major celestines with the disregarded gypsum or anhydrite phases (including the previously determined apfu; Supplementary Table S7). The recalculation includes the theoretical ionic radiuses of the M cations and occupancies of the twelve-coordination site (Supplementary Table S8), which are given in Supplementary Tables S14 and S15. The correlations of the recalculated ionic radiuses (Supplementary Tables S13 and S15) with the a ₀, b ₀, c ₀, and V ₀ parameters show: (1) best fit [C(12); Supplementary Fig. S20 and Table S4] of all the available correlations associated with the problem of the ionic radiuses with UCPs variations [i.e. C(3)–C(8)], including (2) the lowest correlation among all of these correlations [C(13); Supplementary Fig. S21 and Table S4].

Finally, the newly calculated UCPs are characterised based on the apfu (Supplementary Table S12) and are given in Supplementary Table S16. The correlations between these a ₀, b ₀, c ₀ and V ₀ values, with those of the PXRD method (Table 5) have, in this case, an acceptable correlation [C(14); Supplementary Fig. S22 and Table S4], which is significantly better than the previous one associated with this problem [i.e. C(11)]. These observations substantiate the concluding remarks indicating that the existing vacancies in the structures of samples 2, 3 and 5 have no significant impact on their UCPs. Such a good correlation further leads to the conclusion that the complete set of 1–5 samples should be interpreted exclusively as monomineral celestine near-end-members, agreeing with the SEM–WDS, PXRD and IR results. Therefore, we can finally disregard ‘option iii' as a probable cause for the earlier discussed disagreements.

For possibility (b), the major celestines with the disregarded gypsum or anhydrite phases, the results in Supplementary Tables S10 and S16 demonstrate that there is a slightly different ratio between various crystallographic axes, such as c ₀ < a ₀ <b ₀ (samples 1 and 4), a ₀ < c ₀ <b ₀ (samples 2 and 3) and a ₀ = c ₀ < b ₀ (sample 5). Such diverse ratios indicate the presence of axial anisotropy, which has been observed previously and studied for the synthetic and natural end-members of the baryte group of minerals (e.g. Kuang et al., 2017; Li et al., 2018; Girard et al., 2019; Ye et al., 2019; and references therein). This phenomenon appears in conditions experiencing an increase in the temperature environment (axial expansion anisotropy) and elevated formation pressure (axial compression anisotropy) conditions. To the best of our knowledge, there is a lack of similar experiments so far, dealing with the solid solutions between minerals within the baryte group. Accordingly, we further recalculated the variations of the temperature dependence for the UCPs for Clt₉₈Ang₀₂, Clt₉₆Ang₀₄ and Clt₉₄Ang₀₆ celestine–anglesite (Supplementary Table S17) and Clt₉₉Brt₀₁, Clt₉₈Brt₀₂ and Clt₉₇Brt₀₃ celestine–baryte (Supplementary Table S18) solid-solution compositions [mineral name symbols such as Ang, Brt and Clt refer to anglesite, baryte and celestine, respectively (Warr, 2021)]. These recalculations use the combination of the initial experimental data at the ambient pressure for anglesite (Li et al., 2018), celestine and baryte (Ye et al., 2019); according to the suitable cell modification from the Pbnm (II) to the Pnma (I) space group, constrained as follows: a_{0 I} = b_{0 II}, b_{0 I} = c_{0 II}, and c_{0 I} = a_{0 II}. The resulting data (i.e. Supplementary Tables S17 and S18) are further used for the construction of Fig. 3. In addition, we also calculated the expansion degree of the UCPs (i.e. their ratio) of celestine, anglesite and baryte between 320 K and 520 K (Supplementary Table S19). By using its composition and UCPs, we strongly believe that this method can successfully be applied as a geothermometer.

Figure 3. Variations of the temperature (in K) at the ambient pressure for the Clt₁₀₀Ang₀–Clt₉₄Ang₀₆ (column I, left) and the Clt₁₀₀Brt₀–Clt₉₇Brt₀₃ (column II, right) solid-solutions series (in mol. %), according to the recalculated unit-cell parameters (Supplementary Tables S16 and S17): (a) & (e) axis a ₀ (in Å); (b) & (f) axis b ₀ (in Å); (c) & (g) axis c ₀ (in Å); and (d) & (h) volume V ₀ (in ų). Importantly, the proposed method can successfully be used as a geothermometer.

By using the recalculated theoretical ionic radiuses of the M cations (Supplementary Table S13), the excesses of the Sr²⁺, Ca²⁺ and Ba²⁺ apfu, reaching over 1.44 Å, are converted as the total Pb²⁺, i.e. as the anglesite 1.40–3.10 mol.% content in celestine. Similarly, the excesses of Sr²⁺, Ca²⁺ and Pb²⁺ apfu over 1.44 Å are also converted as the total Ba²⁺, i.e. as the baryte 0.41–0.91 mol.% content in celestine (Table 8). Subsequently, the Ang and Brt values (Table 8, plotted in Fig. 3) further allowed the estimation of the formation temperatures for the 1–5 celestine samples by using their determined UCPs (Table 5). The results show that a set of very similar temperatures characterises both the celestine–anglesite and the celestine–baryte solid-solution series. In that manner, we have the validation of the recalculated Ang and Brt values. The results demonstrate that sample 2 underwent the lowest temperatures reaching as low as ~368 K, whereas sample 4 was under the exposure to the highest temperature range reaching ~430 K. Samples 1, 3 and 5 were formed at remarkably similar average temperatures of ~387 K, ~384 K and ~387 K, respectively. Therefore, the sampled celestines were formed at the ~368–430 K (~95–157°C) temperature range, at the ambient pressure conditions. Different correlations between the crystallographic axes lengths are b ₀ < a ₀ <c ₀ (samples 1 and 4) and a ₀ <b ₀ < c ₀ (samples 2, 3 and 5).

Table 8. The estimated formation temperature (in K) of the studied samples based on their recalculated anglesite (Ang) and baryte (Brt) contents (Fig. 3).

Finally, we exemplify the discussed discrepancies between the WDS and the PXRD results. Supplementary Fig. S23 contains five possible variations marked by solid lines (plotted from Fig. 3h; should be taken into account only as an example, because of its validity for each of the other a ₀, b ₀ and c ₀ unit-cell axes, including the Ang contents): (1) volume increase by a temperature increase, including the Brt contents increase; (2) volume increase by a constant temperature, including the Brt contents increase; (3) volume increase by a temperature decrease, including the Brt contents increase; (4) constant volume by a temperature increase, including the Brt contents decrease; and (5) volume increase by a temperature increase, having a constant Brt content. We underline that each of these possible interpretations could be also vice versa. For the most favourable case, we conclude that our celestine 1–5 samples, are mostly with variation type (1), having partially the variation type (5), i.e. that their UCPs increased mainly due to a temperature increase, whereas the Ang and the Brt contents are reaching almost constant values (i.e. having a small content difference of Δ_Ang = 1.70 mol.%, and Δ_Brt = 0.50 mol.%; Table 8), but also following the increased status. Because samples 3 and 5 were exposed to higher temperatures than sample 2, whereas sample 4 was formed at a higher temperature than sample 1, the former records larger thermal expansions of the UCPs than the latter. This is the main reason for the average 2<3<1≈5<4 and 2<3≈1=5<4 UCPs variation behaviours (Tables 5 and 8, respectively). In addition, samples 3 and 5 have a larger expansion of the b ₀ axis, reaching values exceeding those of sample 1; another possible interpretation is that sample 1 was formed at a higher pressure by comparison with samples 3 and 5 because this axis has the largest compressibility (Kuang et al., 2017; Ye et al., 2019). Accordingly, we strongly believe that this argument provides sufficient evidence, to support the interpretation that the formation temperature is the primary cause for the UCPs variations of celestine, whereas different inserted cations into the structure should be interpreted as a secondary factor. To simplify, there is almost a perfect correlation C(15) of the average temperature dependence (Table 8) with the a ₀, b ₀, c ₀ and V ₀ UCPs (Table 5) of 0.965, 0.951, 0.991 and 0.997, respectively (Supplementary Fig. S24 and Table S4). The correlation is superior compared to the previous constraints on C(12) and C(14) correlations describing the celestine composition (Supplementary Table S4).

By using the extrapolated data (Supplementary Tables S17 and S18; Table 8 and Fig. 3) we can estimate the UCPs of the samples formed at a near-room temperature (23°C), and under the ambient pressure conditions (Supplementary Table S20). The resulting data are used exclusively for the characterisation of the relative UCPs (Table 9). In this case, a wide range of correlations between the crystallographic axes is observed: b ₀<c ₀<a ₀ (samples 1 and 4), a ₀<b ₀<c ₀ (sample 2), a ₀<c ₀<b ₀ (sample 3) and c ₀<b ₀<a ₀ (sample 5).

Table 9. Relative unit-cell parameters of the investigated samples, calculated from the ratio of the data outlined in Table 5 and Supplementary Table S20 (marked here with*).

The results have almost perfect correlation C(16) of average temperature dependence (Table 8) vs the a ₀, b ₀, c ₀ and V ₀ ratios (Table 9) of 0.979, 0.926, 0.988 and 0.998, respectively (Supplementary Fig. S25 and Table S4). These are, as expected, very similar to the C(15) correlations, observed relative to the temperature by the UCPs variations (Supplementary Fig. S24 and Table S4). In addition, the linear and the polynomial deviations have almost the same values in these two cases.

The synthesis of data (Tables 8 and 9; Supplementary Fig. S25), allowed us to additionally construct the variation diagram of the relative a ₀, b ₀, c ₀ and V ₀ UCPs (i.e. their ratios) versus the temperature dependence of the samples at the ambient pressure conditions (Fig. 4). Importantly, the resulting variation chart fits (i.e. for all of the five samples together) with the celestine c ₀<b ₀<a ₀ crystal growth expansive behaviour (see figures 11b and 12b in Ye et al., 2019; and Supplementary Table S19). Taking into account solely the crystal expansion induced by the elevated temperature (Table 8), samples 1 and 4 show variations that are similar to that of baryte, whereas samples 2, 3 and 5 show deviations analogous to that of celestine (Ye et al., 2019; and Supplementary Table S19). By interpreting the structural aspect (Antao, 2012), the first group of samples is less distorted, in comparison with the second, confirming the geometry change of the SO₄ and MO₁₂ polyhedral, influenced primarily by the vacancies occupied by the latter ones (already discussed before). On the other hand, taking into account solely the crystal compression by pressure, it appears that sample 4 was also developed at a considerably higher pressure (similarly to sample 1) relative to samples 2, 3 and 5 because of their identical axis variations (Table 9).

Figure 4. Temperature dependence by the relative unit-cell parameters (a ₀ = cyan; b ₀ = red; c ₀ = green; and V ₀ = purple) of the investigated 1–5 samples, at the ambient pressure conditions.

Celestine minerogenesis and formation conditions

Celestine as a low-temperature mineral has frequently been associated with magmatic-hydrothermal activity at temperatures below 300°C (Anenburg et al., 2014). To reconstruct minerogenesis, we take into consideration celestine formation conditions, in particular, the typifying shallow crustal marine, evaporitic palaeoredox environments (Hanor 2000, 2004). Such shallow depths of magma and hydrothermal fluid migration were facilitated by a cluster of faults as part of their conduits in the uppermost basinal segment of the crust (Abouessa et al., 2012; Elshaafi and Gudmundssen, 2016). The multistage volcanism contributed to the extrusion of primitive and evolving light fractions of basaltic magma (Elshaafi and Gudmundssen, 2017).

Our results indicate that the Jebel Eghei celestines were formed at a temperature range of ~368–430 K (~95–157°C), at the ambient pressure conditions (Table 8). Consequently, we tested the following options for the Jebel Eghei celestine minerogenesis (Fig. 5): (1) temperature remains constant by a pressure increase; (2) temperature increases by a pressure increase; and (3) temperature increases by a pressure decrease. Following the experimental data (Kuang et al., 2017; Li et al., 2018; Girard et al., 2019; Ye et al., 2019; Curzi et al., 2022), we tested whether, with the increase of subsurface pressure, the temperatures follow the rising trend. Taking into consideration the age of magmatic events, the proposed increasing temperature trend is somewhat aligned with the stage of minor subsidence during the Early Oligocene and Early Miocene (subsidence stage proposed by van der Meer and Cloetingh, 1993). The coupled temperature–pressure elevations have produced the axial compression (by pressure), which should be simultaneously compensated by the axial expansion (by temperature) so that the UCPs remain as already previously characterised (Table 5). Accordingly, we have chosen the second option (see yellow lines in Fig. 5), because the other two available options lead exclusively towards the smaller UCPs. The increasing slope is approximated accounting for the absence of the appropriate experimental data for the formation pressure. Nevertheless, the increase of the formation pressure allowed us to put subsurface constraints on probable subsurface conditions under which the investigated crystals could have been formed. For example, suppose we adjust the subsurface pressure to fit with the resulting temperatures. In that case, we get the pressure of ca. 250 bars or shallow crustal depths up to only a few kilometres under the Earth's surface (Fig. 5). Accordingly, the temperature range is shifted becoming slightly higher, ca. 372–434K (~99–161°C) fitting into the active hydrothermal fluids that were associated with the volcanic areas. To underline, the celestine-producing hydrothermal event probably occurred during the Miocene producing the brine which was poured initially into the faulted and fractured Eocene lagoon deposits (Abouessa et al., 2012). The subsequent process of precipitation and formation of celestine-bearing druses or nodules (e.g. Anenburg et al., 2014) was probably initiated once the sediments were exposed to the terrestrial conditions. Thus, the celestine druses mapped across the area investigated appear to be the product of the reaction of hypersaline Sr-bearing fluids with gypsum and anhydrite (Hanor 2000, 2004; Forjanes et al., 2020a, 2020b). Both gypsum and anhydrite have been mapped across the investigated areas (Al Gata, Tamed Al Kusur, and Krarat Al Jifah Formations; Vasić and Sheriff, 2007; Abdusahmin, 2020).

Figure 5. Pressure–temperature conditions: (a) Characterised average temperature ranges of ~368–430 K (~95–157°C) stand for all 1 to 5 investigated samples (green tick marks), having the same ambient pressure (Table 8). The sketch shows three scenarios: temperature remains constant by the pressure increase (orange interrupted lines); temperature increases by the pressure increase (yellow full lines with approximate slope); and temperature increases by the pressure decrease (pink dotted lines with approximate slope). The presumed temperature ranges at 250, 500, 750, 1000 and 1250 bars, are represented with blue, red, purple, cyan and magenta colours, respectively. (b and c) The presumed range within sedimentary basins (green full lines) of the investigated 1–5 samples formation (yellow full lines) ranges between (b) 500 and 1000 bars (pink interrupted lines); and (c) 250 and 750 bars (pink interrupted lines). Inset was taken and modified from Hanor (2000) and refers to the shallow crustal and marine environments.

However, subsidence and deposition of sedimentary systems, like these within the Sirt basin (Gumati and Nairn, 1991; van der Meer and Cloetingh, 1993), yield higher pressures than that observed within the surfacing volcanic areas. By presuming a subsidence-related pressure increase, reaching up to 500, 750, 1000 and 1250 bars, the temperature range may increase reaching ~382–444 K (~109–171°C), ~391–453 K (~118–180°C), ~429–482 K (~147–209°C) and ~452–514 K (~179–241°C), respectively (Fig. 5a). In that manner, we provide a reconstruction of the broader temperature range in which the celestine occurred. We further narrowed the celestine formation conditions and its formation depths. As a result, the Jebel Eghei (and probably Al Haruj) celestine minerogenesis occurred at the following subsurface conditions: (1) from ~382 K (~109°C; sample 2), over ~417 K (~144°C; sample 1), ranging to ~482 K (~209^oC; sample 4), at 500, 750 and 1000 bars, respectively (Fig. 5b); or (2) from ~372 K (~99°C; sample 2), over ~412 K (~139°C; sample 1), ranging to ~453 K (~180°C; sample 4), at 250, 500 and 750 bars, respectively (Fig. 5c).

It is to be noted that expected subsidence trends without overpressure zones can induce pressures of 340 bar or ca. 4931 psi, approximating shallow ca. 1500 m or 4921 foot of burial depth using:

$${\rm \rho }_{{\rm ob}}{\rm} = {\rm \rho }_{\rm l}{\rm} + {\rm \rho }_{\rm f}$$

where ρ_ob is total overburden pressure; ρ_l is lithostatic pressure; ρ_f is fluid pressure; as per Tiab and Donaldson (2012). However, the P-T emplacement paths of hydrothermal fluids can have more complex overpressure conditions (Curzi et al., 2022).

Conclusion

Five Late Neogene to Middle Pleistocene aged celestine crystals were analysed and characterised for minerogenesis by applying SEM-WDS, ICP/OES, PXRD and IR methods. Their chemical composition is very similar and fits into the monomineral celestine near-end members in which only 1.6–4.1 at.% of Sr²⁺ contents were substituted by Pb²⁺, Ba²⁺, Ca²⁺ or vacancies (in three samples). The results emphasise that the subsurface temperatures were related to the hydrothermal fluids associated with voluminous intermittent volcanic activities that are the primary cause of the transport of strontium into the Sirt basin. The celestines investigated were most probably formed as a secondary mineral, repositioned into the gypsum and/or anhydrite layer, at temperatures over ~373 K (~100°C) and pressures over 250 bars. The subsequent formation of geodes or druses and nodules is further attributed to the clustering of celestine across sedimentary formations, but under surface conditions (similar to Anenburg et al., 2014). The constraints on the differently inserted Pb²⁺, Ba²⁺ and Ca²⁺ cations into the structure are characterised as a secondary factor for the unit-cell parameters expansion, including the observed slight c ₀ < b ₀ < a ₀ axial anisotropy variations. We believe that these slight anisotropy variations could serve as effective new geothermometric tools capable of providing conditions on multiple mineral growth stages. This is of particular importance in exploring the enigmatic relationships between Sr mineralisation and igneous rocks (for a discussion of the complex relationships see Zhu et al., 2022).

Three celestine samples have a greenish-blue to blue colour, whereas the other two samples are of blue–green (cyan) colour. In addition, the celestine minerogenesis study of the Jabal Eghei area has confirmed that the coefficients of the regression-R² may be a very useful accessory tool for comparing mineral crystal-structure parameters, composition and genesis.

Except for the larger vacancy content in sample 5, the most prominent difference in the colour purity ranges between samples 3 and 5. These samples have very similar results from the various methods used in this study, so further in-depth investigations would be needed to resolve this problem.