Shortwave infrared hyperspectral imaging as a novel method to elucidate multi-phase dolomitization, recrystallization, and cementation in carbonate sedimentary rocks

Carbonate rocks undergo low-temperature, post-depositional changes, including mineral precipitation, dissolution, or recrystallisation (diagenesis). Unravelling the sequence of these events is time-consuming, expensive, and relies on destructive analytical techniques, yet such characterization is essential to understand their post-depositional history for mineral and energy exploitation and carbon storage. Conversely, hyperspectral imaging offers a rapid, non-destructive method to determine mineralogy, while also providing compositional and textural information. It is commonly employed to differentiate lithology, but it has never been used to discern complex diagenetic phases in a largely monomineralic succession. Using spatial-spectral endmember extraction, we explore the efficacy and limitations of hyperspectral imaging to elucidate multi-phase dolomitization and cementation in the Cathedral Formation (Western Canadian Sedimentary Basin). Spectral endmembers include limestone, two replacement dolomite phases, and three saddle dolomite phases. Endmember distributions were mapped using Spectral Angle Mapper, then sampled and analyzed to investigate the controls on their spectral signatures. The absorption-band position of each phase reveals changes in %Ca (molar Ca/(Ca + Mg)) and trace element substitution, whereas the spectral contrast correlates with texture. The ensuing mineral distribution maps provide meter-scale spatial information on the diagenetic history of the succession that can be used independently and to design a rigorous sampling protocol.


Overview of the diagenetic features in the Cathedral Formation
In the southern Rocky Mountains, the Cathedral Formation consists of light-grey limestone, medium-grey to tan finely-crystalline replacement dolomite (RD), and white coarsely-crystalline saddle dolomite (SD). Dolomitization is most pervasive proximal to the Cambrian platform margin and dolomite grades laterally to limestone towards the northeast (Fig. 1d) 26 . Such dolomite bodies are typically non-stratabound (inclined-to-bedding) at their cores with stratabound (bedding-parallel) margins [24][25][26] . Cement-supported breccias and zebra textures are widespread in the Kicking Horse Rim area, with local occurrences of talc, magnesite, and Mississippi Valley-type (MVT) minerals 30,31,36 . These minerals are absent to the northeast, but zebra textures and cement-supported breccias are locally common 26,37 .
This study focuses on an outcrop, 240 m in width and 40 m in height, at Whirlpool Point that includes a faultcontrolled dolomite body in the Cathedral Formation. The outcrop contains a fault that is oriented at 028/52, has a normal offset of 30 cm, and intersects the formation 100 m from the east end (Fig. 1e). At the core of the dolomite body, coarsely-crystalline breccias extend 25 m into the hanging-wall, 5 m into the footwall, and are limited to the fault damage zone. Proximal to the fault, the hanging-wall (146/32) comprises cement-supported breccias and bedding-inclined zebra textures. The medial part of the hanging-wall includes fabric-retentive dolomitized microbial bindstone with bedding-parallel and rare bedding-inclined zebra textures. At the margin of the dolomite body, the hanging-wall includes fabric-retentive, finely-crystalline dolomitized peloidal wackestone with rare bedding-parallel zebra textures. The lower part of the Cathedral Formation includes a sharp, bedding-parallel contact with a 2 m thick bed of limestone. In the footwall (179/25), the upper part of the formation is similar to the hanging-wall, but cement-supported breccias and bedding-inclined zebra textures are rare. The medial part of the footwall includes fabric-retentive dolomitized microbial bindstone with bedding-parallel zebra textures that grade laterally to fabric-retentive, finely-crystalline dolomite at the margin (Fig. 1e).

Methods
Collection and processing of the infrared reflectance data. A set of four SWIR (930-2508 nm) spectral images were acquired on June 12, 2018 between 11 am and 2 pm using a Specim SisuROCK hyperspectral scanner (a linescan imager) that is mounted on a rotary stage for wall rock imaging. Integration time varied from 5 to 10 ms depending on the time of acquisition and it required 30 s for the stage to rotate 90°. The scanner contains a 256 spectral by 320 spatial pixels mercury-cadmium-telluride detector array that acquires data at a 6.3 nm sampling interval and a 10 nm spectral bandwidth. Two Spectralon panels of 2% and 99% reflectance were positioned in each scene. Data      www.nature.com/scientificreports/ the outcrop. For each scene, radiance data was obtained by applying appropriate gain and offset and conversion to reflectance, then an empirical line correction was applied based on the known reflectance of the Spectralon panels relative to their measured radiance spectrum. The latter were obtained as the mean radiance spectrum for a 20 × 20 pixel area over the panel (nominal pixel size of the acquired imagery = 5 cm). The empirical line method has the advantage of correcting for the influence of the atmosphere on the target radiance [38][39][40] . In contemporaneous studies with the same camera 7,41 , the wavelength position of SWIR absorptions of a National Institute of Standards and Technology referenced Mylar standard were accurate within 1 nm. Following calibration of the spectral data to reflectance, bands with the poorest signal to noise (e.g., 930-991, 1295-1435, 1735-1998, and 2461-2508 nm) were removed from the ensuing analysis. The four images were then spatially co-registered using tie points, resulting in a single image for further analysis. Next, an iterative spatial spectral filter was used to compare the spectral similarity of spatially adjacent pixels within a 3 × 3 pixel window [42][43][44] . When the spectral signatures were within a predefined similarity threshold, an average spectrum was substituted for the original data, thereby reducing the intra-class spectral variability. Lastly, areas in shadow and the calibration panels were masked.
Mineralogical and lithological information in the imagery was obtained by the extraction of endmember spectra and their distributions were mapped. To derive an image endmember set, spatial-spectral endmember extraction (SSEE) was used to divide the image into equal spatial subsets (each subset = 7 × 7 pixels) [42][43][44] . This method is designed to discern spectrally similar endmembers that occupy different portions of the scene. The endmember set derived from SSEE was clustered and labelled to derive final endmember sets for mapping. For clustering, we used a tree cluster that recursively merges a pair of clusters based on a similarity measurement. To start, each endmember was treated as an individual cluster and endmembers that are most similar were successively merged. In this study, the Spectral Angle (SA) between two endmembers was used as the measure of similarity. A minimum SA threshold was defined to stop the merging process that took place when all pairwise clusters had a similarity greater than the threshold. To address the spectral variability of the extracted endmembers, the tree cluster tool was applied twice on the given data. The first time, using all endmembers, a SA threshold of 0.2 radians produced clusters that capture the broad material classes, namely non-geological (e.g., panels, vegetation, weathering) and geological. The next level of clustering focused on the geological class to capture subclasses and define multiple geological endmember clusters. In this case, a smaller SA threshold (0.05 radians) was used because these endmembers are more spectrally similar. Clustered endmembers were then averaged to obtain an individual endmember that represents the given class, contributing to an endmember set of thirteen geological endmembers. This clustering process was data-driven.
After accounting for the spectral similarity between classes, these thirteen endmembers were condensed into four groups (Fig. 2a,b): limestone (Lst), two groups of replacement dolomite (RDa, RDb), and saddle dolomite (SD). Groups were labelled based on spectral interpretations that were supported by field and petrographical observations. Group Lst was defined by all pixels with a carbonate absorption-band position > 2330 nm and was validated in the field using dilute hydrochloric acid. Group RDa includes the spectral endmembers that correspond to light-and medium-grey replacement dolomite. Endmembers that corresponded to clasts, bedding, and bedding-parallel fractures within the RDa intervals were also included. Group RDb includes the endmembers that correspond to light-brown replacement dolomite and the alteration rims along the margins of the saddle dolomite intervals. Group SD includes three subgroups (SDa, SDb, SDc), labelled based on their paragenesis, that correspond to white, coarsely crystalline, saddle dolomite. Their paragenesis was determined by the relative positions of each endmember in macro-pores that were validated by petrographical analyses.
Distribution maps of the spectral endmembers. Mapping of the endmember spectra resulted in two image products. The first examined the distributions of limestone (Lst), replacement dolomite (RDa, RDb), and saddle dolomite (SD) and is suited for a synoptic view of the outcrop. The second, more detailed image product, examined the distribution of the saddle dolomite subgroups (SDa, SDb, SDc). In both instances, mapping the distribution of each endmember was conducted using a spectral angle mapper (SAM) algorithm that treats spectra as multidimensional vectors and computes the angle between spectral pairs 45 . For this purpose, the spectrum from each pixel of the image after processing was compared to that of each endmember. Spectra with the smallest SAM angles indicate the greatest similarity. In the first image product, the SAM results for RDa, RDb, and SD are presented as a red-green-blue composite where a higher color hue corresponds to a higher spectral similarity to the given endmember. If two or more of the endmembers predominate within a given pixel, then a color other than RGB is seen. In the second image product, the SAM results for SDa, SDb and SDc are classified in such a way that the endmember of highest similarity to that of the given pixel is assigned to that pixel. Consequently, only the colors that were assigned to the endmembers are seen in the second image product and when the SAM angle exceeds 5 degrees the pixel is not classified.
Sampling, petrography, and geochemical analyses. Fieldwork and sampling were conducted over two field seasons. Prior to obtaining the hyperspectral reflectance data, Stacey et al. 26 collected 72 samples from the Cathedral Formation in the Whirlpool Point area; 35 of these samples were systematically taken from the roadcut at ~ 2 m intervals along a 62 m logged section. After processing the hyperspectral data, an additional 33 samples were taken from the roadcut to support the analysis of the reflectance data and to evaluate specific features and trends within the mineral distribution maps. Samples were impregnated with blue epoxy and prepared as polished sections that were partially stained with alizarin red-S and potassium ferricyanide 46 . Polished sections were examined under plane-and cross-polarized light and then analysed using a CITL Mk5 cold cathodoluminescence (CL) system (operating conditions 15-20 kV and 350-450 μA) mounted on a Nikon Eclipse LV100N POL microscope. Dolomite crystal textures are described according to Sibley  www.nature.com/scientificreports/ Thirty-eight representative samples from endmembers RDa, RDb, and SD were analysed for mineralogical composition by powder X-ray diffraction (XRD) using a Bruker D8 Advance diffractometer. Quartz was added as a standard and samples were scanned at 40 kV and 30 mA from 5° to 70° 2θ in 0.02° increments. %Ca is calculated based on Lumsden 48 and degree of ordering is based on Goldsmith and Graf 49 . Two polished mounts that included each of the SD subgroups (SDa, SDb, SDc) were analysed by quantitative evaluation of minerals by scanning electron microscopy (QEMSCAN) and energy-dispersive X-ray (EDX) spectroscopy using an FEI Aspex eXstreme equipped with Bruker 5030 EDX detectors and an iExplorer software suite. An initial 10 × 10 mm mineral map was created for each sample at a stepping interval of 50 μm and 4 × 4 mm areas of interest were mapped at a stepping interval of 4 μm. An X-ray spectrum was generated for each point, matched against a standard library, and the map was constructed by assigning the library mineral to each point.
Twenty-four representative samples from endmembers RDa, RDb, and SD were analysed for trace elements by inductively coupled plasma mass spectrometry (ICP-MS) using an Agilent 7700× at the Advanced Isotope Geochemistry and Cosmochemistry Suite, The University of Manchester. Two polished sections that included each of the SD subgroups (SDa, SDb, SDc) were analysed by electron probe micro-analysis (EPMA) using a Cameca SX100. An initial 2 × 10 mm map was created at a stepping interval of 10 μm and 1.536 × 1.536 mm areas of interest were mapped at a stepping interval of 3 μm. Ca (Kα; PET) and Mg (Kα; TAP) were analysed at 15 kV, 10 nA, and a dwell time of 100 ms using calcite and magnesite as standards. Fe (Kα; LLIF) and Mn (Kα; LLIF) were analysed at 15 kV, 200 nA, and a dwell time of 200 ms using fayalite and tephroite as standards.  www.nature.com/scientificreports/ Results SWIR hyperspectral imaging. Each of the endmembers can be discerned based on their absorption-band positions and spectral contrast (the difference between the peaks and the valleys in the spectrum). Phase Lst presents an absorption-band position of 2333 nm, akin to calcite and unique among all endmembers. RDa, RDb, and SD have absorption-band positions of 2321, 2317, and 2314 nm, respectively, and display increasing spectral contrast (Fig. 2c). Phase SD includes SDa, SDb, and SDc that have absorption-band positions of 2316, 2313, and 2314 nm, respectively, and display increasing spectral contrast (Fig. 2d). The absorption-band asymmetry for Lst, RDa, RDb, and SD are 86, 150, 104, and 86 nm, respectively (Fig. 2c) and the asymmetry for SDa, SDb, and SDc are 96, 84, and 84 nm, respectively (Fig. 2d).
At the margin of the fault-controlled dolomite body, a sharp, bedding-parallel contact occurs between Lst, RDa, and RDb (Fig. 3a, b). Scattered RDa pixels are located below the limestone-dolomite contact but RDb and SD are absent. Phase Lst is absent above this contact and throughout the remainder of the hanging-wall and the footwall. The distal parts of the hanging-wall comprise a mixture of RDa and RDb, the contacts between which follow the bedding (Fig. 3b, c). Phase SD is rare in the distal parts of the hanging-wall and is typically restricted to isolated beds of bioturbated wackestone. The spatial distribution of RDb pixels correlate with the occurrence of SD pixels, whereas RDa pixels are typically not in contact with SD pixels (Fig. 3b, c). This isolated SD bed illustrates the paragenetic sequence of each of the constituent SD phases. SDa lines the margins of the SD intervals and is post-dated by SDb. SDc is rare, but this phase is located at the centers of the SD intervals (Fig. 3b, d).
At the core of the dolomite body, the hanging-wall largely comprises cement-supported breccias and zebra textures (Fig. 3e). Consequently, phase SD is more abundant at the core of the dolomite body relative to at the margins. Isolated clasts of RDa and RDb are suspended and fully surrounded by SD (Fig. 3e, f). Typically, RDb pixels are located adjacent to SD pixels, but contacts between RDa and SD are locally common (Fig. 3e, f). The spatial distributions of the constituent phases of SD correlate with fault proximity. SDa and SDb are located throughout the hanging-wall but increase in abundance towards the core of the dolomite body ( Fig. 3e, g). SDc post-dates these phases and is located at the centers of macro-pores such as breccias, fractures, zebra textures, and bedding planes. SDc increases in abundance towards the fault (Fig. 3e, g).
Petrography. Petrographical analysis of samples from the Cathedral Formation identified several diagenetic phases based on their crystal size, texture, fluid inclusions, CL properties, and mineral associations (Table 1). In addition to the host limestone, two phases of RD and three phases of SD were identified. RDa includes finelycrystalline (20-150 μm), planar-e to planar-s dolomite with dull-purple luminescent cores and bright-orange to dull-red luminescent rims (Fig. 4a, b). RDb includes medium-crystalline (100-400 μm), planar-s to non-planara dolomite with dull-red luminescent cores and dull-to moderate-red luminescent rims (Fig. 4a, b). The crystal size distribution in RDa has a normal distribution with a mode of 84 μm, a mean of 86 μm, and a standard deviation (±) of 24 μm (Fig. 4c). The crystal size distribution in RDb is broader than RDa and has a slight negative skew (− 0.3). The modal crystal size of RDb is 284 μm with a mean of 267 μm (± 59 μm; Fig. 4c). RDa includes trace clay minerals, detrital quartz, organic matter, and pyrite that are rare in RDb. Phase SD consists of three separate phases of non-planar (saddle) dolomite that have different petrographical characteristics (Table 1).
SDa includes medium-crystalline (250-550 μm), non-planar dolomite that typically grew syntaxially on RDb (Fig. 4d). SDb includes coarsely-crystalline (up to 4500 μm), non-planar dolomite crystals that are abundant in cement-supported breccias, zebra textures, and fractures (Fig. 4d, e). SDb crystals are elongate, oriented normal to the cavity wall, and are polymodal in size depending on the size of the cavity in which they were precipitated. SDa and SDb are unzoned with a dull-to medium-red luminescence. SDc includes coarsely-crystalline, non-planar dolomite that are largely indistinguishable from SDb based on their textural properties. SDc crystals, however, have a characteristic dull-to bright-red and dull-to bright-orange oscillatory zonation (Fig. 4e). SDc typically form rims (250 to 1250 μm thick) that are nucleated on SDb crystals in zebra textures and fractures. Individual SDc crystals (up to 4500 μm) are common in the central parts of cement-supported breccias.
Geochemistry. QEMSCAN indicates that the RD phases comprise 97.87% dolomite, 1.67% clay minerals (trace illite) and muscovite, 0.41% quartz, 0.02% calcite, 0.01% pyrite, and 0.02% other minerals (Fig. 5a, b). In contrast, phase SD comprises 99.91% dolomite, 0.06% clay minerals and muscovite, 0.02% pyrite, and 0.01% calcite (Fig. 5a, c). Although RDb and SDa are clearly discernible in outcrop and hand-sample due to their color and crystallinity, EDX spectroscopy indicates that there is minimal contrast in the abundances of Ca, Mg, and Fe between these phases (Fig. 5d-f). In contrast, SDb and SDc cannot be confidently distinguished in outcrop and hand-sample, but they are clearly distinguished by their composition (Fig. 5g-i).
Each of the constituent phases of SD were further analysed by EPMA (    (Tables 1, 2). Recrystallization is associated with a change from planar-e to planar-s dolomite in RDa to planar-s to non-planar-a dolomite in RDb (Fig. 4a, b). Given that RDa and RDb are cross-cut by low-amplitude, bedding-parallel stylolites, replacement dolomitization is interpreted to have occurred during shallow burial 26,50 . Saddle dolomite (SD) grows within pores in RDa and RDb and, therefore, post-dates them. SDa is consistently located at the margins of these pores and is overgrown and postdated by SDb (Fig. 4d). In thin-section, SDa forms a syntaxial rim that is in optical continuity with RDb; a common feature of cavity-filling cements with the same mineralogy as the cavity-wall 37 . SDa and SDb have similar compositions, gradational contacts, and were likely derived from a similar fluid-flow event. SDb crystals are an order of magnitude larger than SDa (Table 1) due to competitive crystallization; a process by which favourably-oriented crystals obstruct the growth of poorly-oriented crystals 51 . Although SDc is not present throughout the outcrop, it consistently nucleates on and, thus, postdates SDb. SDc is restricted to the central parts of vugs, fractures, zebra textures, and cementsupported breccias (Fig. 3).

Spatial distribution of each diagenetic phase. Dolomitization fronts in the Cathedral Formation and
the underlying Mount Whyte Formation at Whirlpool Point are interpreted to have "retreated" over time due to the occlusion of porosity from repeated fluid-pulses [23][24][25][26] . In this model, the core of the dolomite body is younger than the margins; with each successive fluid-pulse contributing to the recrystallization of earlier phases during the cementation of the dolomite body [23][24][25][26] . This retreating dolomitization front is associated with increased dolomite stoichiometry and ordering towards the core of the dolomite body 24 . The cement-supported breccias and associated zebra textures imaged in this study are largely restricted to the core of the dolomite body and are interpreted to have formed as a final event when the occlusion of porosity gave rise to high pore-fluid pressures and the rupturing of the formation during seismic valving 26 .
The mineral distribution maps, derived from the hyperspectral data, reveal outcrop-scale heterogeneities that provide evidence of how fault-controlled dolomitization in the Cathedral Formation progressed and terminated. Each paragenetic stage increases in abundance from the margin of the dolomite body to the core and their spectral signatures indicate that the %Ca decreases between each phase (Fig. 6). At the margin, the dolomitization front (contact between Lst and RDa) is sharp and bed-parallel (Fig. 3b, c). RDa grades laterally to RDb and the phases have a patchy distribution throughout the outcrop. Bedding planes can be identified by the distribution of RDa and RDb (Fig. 3b, c), which suggests that that they acted as permeability pathways that circulated the initial dolomitizing fluids during replacement dolomitization 52 . RDb has a spatial relationship to SD, which suggests that recrystallization occurred during the later fluid-pulses that also precipitated SD (Fig. 3). Consequently, there are textural changes from the nonstoichiometric RDa (concentrically zoned, poorly-ordered, planar-e to planar-s   www.nature.com/scientificreports/ dolomite) at the margins of the dolomite body to the more stoichiometric RDb (weakly zoned, well-ordered, planar-s to nonplanar-a dolomite) at the core. SD is present throughout the dolomite body, but it is more abundant at the core relative to at the margins and in the hanging-wall of the fault relative to in the footwall (Fig. 3a). Each of the constituent phases of SD, which are indistinguishable in visible light, also increase in abundance from the margin of the dolomite body to the core (Fig. 3a). SDa occurs as isolated pixels in RDa and RDb, and lines the margins of vugs, zebra textures, and cement-supported breccias (Fig. 3d, g). SDb has a gradational contact with SDa and is restricted to the central parts of these rock textures. Finally, SDc precipitated in macro-pores that were proximal to the fault (Fig. 3g), but occasionally occurs along fractures, stylolites, and bedding planes that were likely the few remaining permeability www.nature.com/scientificreports/ pathways after brecciation and cementation 50 . Consequently, the back-stepping and paragenesis of the SD subgroups could only be resolved using hyperspectral imaging because the endmembers are indistinguishable at the outcrop-scale and thin-sections are at too small a scale to capture their spatial distributions.
Compositional characteristics derived from the reflectance spectra. The position of the 2315-2335 nm absorption-band is commonly used to differentiate dolomite from calcite in the laboratory 10,11 and in remotely-sensed imagery 18,21,53 . Van der Meer 13 suggested that this absorption-band has a linear relationship with the degree of dolomitization and Zaini et al. 14,15 applied this approach to track the band positions of synthetic mixtures of calcite and dolomite. This methodology, however, has not been used to infer the %Ca of multiple phases of dolomite with varying compositions. RDa is the most calcium-rich phase of dolomite identified in the Cathedral Formation (%Ca = 50.21; Table 2) and has an absorption-band positioned at the longest wavelength (λ = 2321 nm; Fig. 2c). In contrast, SDb is the most magnesium-rich phase (%Ca = 49.67; Table 2), with an absorption-band positioned at the shortest wavelength (λ = 2313 nm; Fig. 2d). When SDc is excluded, the phases with intermediate compositions plot along a linear trendline with a slope of 0.07%Ca/nm and a correlation coefficient (r) of 0.94, p < 0.05, n = 5 (Fig. 6a). SDc deviates from this trend because it includes substantial trace element substitution for Mg ([Fe + Mn] = 13,320 ppm; Fig. 6b, c). Increased [Fe + Mn], each with a mass greater than Ca and a radius between Mg and Ca, shift the carbonate absorption-bands to longer wavelengths 11,16 . The [Fe + Mn], however, has a poor correlation (r = 0.50, p > 0.05, n = 6) with the position of the ~ 2315 nm absorption-band of each of the dolomite phases in the Cathedral Formation (Fig. 6d). In the laboratory, Gaffey 11 documented a broad Fe +2 absorptionband at 1200 nm and Mn 2+ bands from 300 to 800 nm, but these features are difficult to identify in the field due to spectral range limitations of the equipment 17 . Strontium, which substitutes for Ca in dolomite 54 , reveals the Although each of the RD and SD phases are > 97% dolomite, the asymmetry of the ~ 2315 nm absorption-band generally relates to the volume of non-carbonate grains identified in thin-section. RDa (S = 150 nm) and RDb (S = 104 nm) include trace clay minerals and detrital quartz that are absent in each SD phase (S = 84-96 nm). QEMSCAN indicates that RD includes 1.67% clay minerals and muscovite, whereas they are negligible in SD (Fig. 5b, c). Clay minerals in RD reflect the composition of the precursor limestone, whereas their absence in SD is consistent with it being precipitated as a void-filling cement. A muscovite absorption feature at ~ 2200 nm could not be distinguished, and it is unclear if non-carbonate minerals affect the absorption-band asymmetry of each phase.

Influence of textural characteristics on the reflectance spectra. Several laboratory-based studies
have demonstrated that the textural properties of minerals and rocks impact the surface scattering and volume scattering of light 11,13,14,22 . Gaffey 11 and Zaini et al. 14 , for example, studied the reflectance of powdered carbonate minerals that were sieved to various sizes, typically < 500 μm, and showed that finer powders have a higher overall reflectance with shallower absorption features, whereas coarser powders have a lower overall reflectance with deeper absorption features. The degree of compaction in rocks relative to powders, as well as the occlusion of porosity, can minimize these differences and it is not always possible to infer the grain or crystal size of rocks from their spectra.
In the Cathedral Formation, the endmembers present three broad groups based on their spectral contrast (Fig. 2). RDa, with a mean crystal size of 86 μm, has the lowest overall reflectance and absorption-band depths Absorption band position (nm) 14  www.nature.com/scientificreports/ (Fig. 2a, c). Organic matter (OM), which reduces the reflectance of a sample and attenuates its absorptionbands 13 , is locally common in RDa and rare in RDb (Fig. 3a, b). Consequently, the presence of OM in RDa may contribute to its low overall reflectance. RDb and SDa (mean crystal size = 267 μm and 400 μm, respectively) have similar overall reflectance and absorption-band depths (Fig. 2). These phases are more coarsely crystalline than RDa and their spectral contrasts are intermediate between RDa and the remaining SD phases. SDb and SDc have the largest crystals (mean size = 2000 μm), roughly five times larger than SDa, and they display the highest overall reflectance and absorption-band depths (Fig. 2b, d). In this case, it is likely that the specular reflectance from crystal facets dominates over volume scattering, thus, explaining their high reflectance 55 .

Discussion
The results of this work demonstrate that hyperspectral imaging is an invaluable tool that reveals mineralogy, composition, and texture in carbonate rocks at a macro-scale. Nevertheless, care must be taken to consider the scale of observation, or spatial resolution, because each imaged pixel comprises several phases that contribute to the reflectance 56,57 . Diagenetic phases in the Cathedral Formation, for example, comprise crystals that are below the image resolution (5 cm). Although each of the phases of dolomite are monomineralic, they are represented by an endmember because they are compositionally and texturally distinct. RDa and RDb co-exist in samples at the thin-section scale and their spectral signatures are, therefore, a linear average of their relative surface abundances. For this reason, applying the laboratory-based spectra of a single dolomite endmember to our field-based imagery would not have captured the compositional and textural diversity that is present at the outcrop-scale without petrographical and geochemical calibration. Similarly, Beckert et al. 18 cautioned against the over-reliance on published spectral libraries, based on pure, end-member, minerals (i.e., stoichiometric well-ordered dolomite), to discriminate carbonate rocks in the field. Given that most sedimentary dolomite is nonstoichiometric with 48.0-62.5%Ca 48,58 , consists of different crystal shapes and boundaries 47 , and includes different volumes of substituted trace elements, the application of a single dolomite endmember to field-or satellite-based imagery without calibration should be treated with caution. Initial observations of the Cathedral Formation suggests that the succession comprises only limestone, a single phase of RD, and a single phase of SD. Detailed petrographical and geochemical analyses 26 , however, revealed that the succession includes several diagenetic phases. Accordingly, this study examined the extent to which hyperspectral imaging can be used to map multiple phases of dolomite at the outcrop-scale and to determine their cross-cutting relationships; features that have traditionally only been revealed petrographically. The method presented here yields a map that illustrates the spatial distribution of each diagenetic phase in an outcrop; a feat that could not be achieved with photogrammetry or Lidar, even if supported by a dense sampling campaign 17,18 . As a result, hyperspectral imaging provides a geological toolkit that facilitates systematic sampling and improves confidence that all of the phases in the paragenetic history have been sampled, thereby allowing for macro-scale fluid flow pathways to be determined. This has economic implications because structurally-controlled deposits of magnesite and MVT-minerals (e.g., Mount Brussilof, Kicking Horse, Monarch) are associated with "hydrothermal sparry dolomite" in the Cathedral Formation 30,59,60 . Their effective exploitation, therefore, requires a robust understanding of the spatial distributions of each of the SD subgroups and their relationship to mineralization.
The methodological study presented here, with a wavelength accuracy greater than 1 nm 7,41 , tested whether hyperspectral imaging can be used to discriminate a narrow range of dolomite stoichiometries and textures in a Middle Cambrian succession. We distinguished five phases of dolomite, ranging from 49.67 to 50.21%Ca, that are approaching a stoichiometric, well-ordered endmember following 100's of millions of years of diagenesis and several kilometers of burial [23][24][25][26] . The progressive recrystallization of dolomite over geological time is driven by mineralogical stabilization during burial that increases dolomite stoichiometry and cation ordering 61 . Consequently, this methodology has profound transferability to other, younger, successions that have not been subject to such burial, recrystallization, and mineralogical stabilization. This includes the Cenozoic "island-type" dolomite bodies that have a wider range of dolomite stoichiometries [62][63][64] . With careful calibration, hyperspectral imaging can provide an upscaled view of the distribution of each of these diagenetic phases in an outcrop, akin to a geocellular model, that cannot be replicated by conventional geological methods.

Conclusions
Exposures of Middle Cambrian strata in the WCSB offer an unparalleled natural laboratory to unravel the extent and timing of dolomitization, recrystallization, and cementation in carbonate rocks. Hyperspectral imaging, in conjunction with the detailed analysis of samples from a fault-controlled dolomite body, has led to the following important conclusions: (1) The position of the ~ 2315 nm absorption-band is an effective analogue for the %Ca of dolomite and this can be used to map outcrop-scale variations in dolomite stoichiometry. (2) The relationship between %Ca and absorption-band positions is convoluted in dolomite with high trace element substitution (SDc; %Ca = 50.53; [Fe + Mn] = 13,320 ppm), therefore, non-stoichiometric dolomite can be distinguished using hyperspectral imaging. (3) The spectral contrast of the reflectance profile, which accounts for overall reflectance and absorption-band depths, correlates with the textural properties (e.g., crystal size, boundary-shape) of each dolomite phase, thus, enabling their discrimination and mapping. (4) The paragenesis and spatial distributions of the RD and SD phases in the Cathedral Formation support the prior interpretation that the dolomitization front "retreated" towards the fluid source during the ensuing recrystallization and cementation of the dolomite body. www.nature.com/scientificreports/ The results of this study demonstrate that SWIR hyperspectral imaging is capable of discerning subtle diagenetic heterogeneities in carbonate rocks, beyond the routine identification of calcite and dolomite. Consequently, robust multi-scale studies can be conducted through the targeted sampling of individual diagenetic phases for further petrographical and geochemical analyses. www.nature.com/scientificreports/