*Article* **Tracking the Origin and Evolution of Diagenetic Fluids of Upper Jurassic Carbonate Rocks in the Zagros Thrust Fold Belt, NE-Iraq**

**Namam Salih 1,2,\*, Alain Préat 3, Axel Gerdes 4, Kurt Konhauser <sup>5</sup> and Jean-Noël Proust <sup>6</sup>**


**Abstract:** Utilizing sophisticated tools in carbonate rocks is crucial to interpretating the origin and evolution of diagenetic fluids from the Upper Jurassic carbonate rocks along the Zagros thrust-fold Belt. The origin and evolution of the paleofluids utilizing in-situ strontium isotope ratios by high resolution laser ablation ICP-MS, integrated with stable isotopes, petrography and fieldwork are constrained. Due to the lack of information on the origin of the chemistry of the fluids, the cements that filled the Jurassic carbonate rocks were analysed from the fractures and pores. This allowed us to trace the origin of fluids along a diagenetic sequence, which is defined at the beginning from the sediment deposition (pristine facies). Based on petrography and geochemistry (oxygen-, carbonand strontium-isotope compositions) two major diagenetic stages involving the fluids were identified. The initial stage, characterized by negative <sup>δ</sup>13CVPDB values (reaching <sup>−</sup>10.67‰), involved evaporated seawater deposited with the sediments, mixed with the input of freshwater. The second stage involved a mixture of meteoric water and hot fluids that precipitated as late diagenetic cements. The late diagenetic cements have higher depleted O–C isotope compositions compared to seawater. The diagenetic cements display a positive covariance and were associated with extra<sup>δ</sup>13CVPDB and <sup>δ</sup>18OVPDB values (−12.87‰ to <sup>−</sup>0.82‰ for <sup>δ</sup>18OVPDB and <sup>−</sup>11.66‰ to <sup>−</sup>1.40‰ for δ13CVPDB respectively). The distinction between seawater and the secondary fluids is also evident in the 87Sr/86Sr of the host limestone versus cements. The limestones have 87Sr/86Sr up to 0.72859, indicative of riverine input, while the cements have 87Sr/86Sr of (0.70772), indicative of hot fluid circulation interacting with meteoric water during late diagenesis.

**Keywords:** origin of diagenetic fluids; strontium isotope-laser ablation ICP-MS; upper jurassic carbonate rocks; ZTFB; NE-Iraq

#### **1. Introduction**

The Upper Jurassic Barsarin formation is located along the Zagros thrust-fold belt (ZTFB) in NE Iraq, and is considered a giant undiscovered Jurassic source rock [1]. Several studies have documented the main petrographic features of the Barsarin formation in order to characterize the paleoenvironment, but without any attention to diagenesis. The formation comprises laminated limestones and dolomitic limestones, in places cherty, with autoclastically brecciated beds admixed with shaly and marly materials [2]. The microcrystalline quartz (chert) associated with carbonate rocks have different silica sources,

**Citation:** Salih, N.; Préat, A.; Gerdes, A.; Konhauser, K.; Proust, J.-N. Tracking the Origin and Evolution of Diagenetic Fluids of Upper Jurassic Carbonate Rocks in the Zagros Thrust Fold Belt, NE-Iraq. *Water* **2021**, *13*, 3284. https://doi.org/10.3390/ w13223284

Academic Editor: Domenico Cicchella

Received: 18 October 2021 Accepted: 12 November 2021 Published: 19 November 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

including a biogenic origin [3], silica-enriched seawater [4], and/or linked to silica input via a river source [5].

The origin and evolution of the diagenetic fluids that affected the Barsarin sediments are still largely unknown due to a lack of detailed studies. However, the fluids can be traced by geochemistry, and particularly by analyses of stable (carbon and oxygen) and radiogenic (strontium) isotopes. Oxygen isotope compositions can provide insights into the origin of mixed fluids or cross-formation water flows, and also to characterize paleo-temperatures and rainfall properties, such as the amount, seasonally, and moisture sources (e.g., [6–8]). Carbon-isotope compositions can provide information on effective rainfall and weathering processes [8]. The conventional 87Sr/86Sr isotope composition is classically used to infer the main sources of radiogenic strontium in sedimentary basins related to the pulses of midoceanic hydrothermal fluxes [9] or to infer overprinting by radiogenic dolomitizing fluids (e.g., [10]). It allows for the determination of continental riverine input [11]. However, the classical conventional technique to measure radiogenic isotopes is limited to the achievable spatial resolution of 87Sr/86Sr records.

To address the time specific gap in tracking the marine limestone and associated diagenetic fluids in the Barsarin formation, a new approach is developed here. It consists of coupling a high-resolution laser ablation (LA-) MC-ICP-MS analysis with fieldwork data, petrography, and geochemistry (oxygen-carbon stable isotopes). The remarkable high precision of laser ablation analysis is utilized to measure, for the first time in our sediments, the absolute radiogenic composition of the Jurassic carbonate rocks of the Barsarin formation, and evaluate the precise origin, involvement and evolution of different fluids during diagenesis that modified the original composition of the host limestones. This will also give information on the dynamic of the diagenetic fluids during the Zagros orogeny within the ZTFB.

#### **2. Geological Setting**

The Upper Jurassic high-folded zone, considered as a part of the Zagros fold-thrust belt, belongs to a NE-SW basin (Figure 1), inherited from the tectonic activity of the Arabian plate and the opening of the southern Neo-Tethys ocean [12]. This ocean was characterized by sea level fluctuations that affected evaporitic sabkhas in the basin during Jurassic-Cretaceous times. NE Iraq formed an euxinic basin separated from the Neo-Tethys by a rifted area where shallow water carbonates formed [12]. It is generally assumed that the Barsarin formation records an euxinic environment, subjected to a low subsidence that characterized the Neo-Tethys during this period.

The Barsarin formation is composed of laminated limestones and dolomitic limestones, and in places brecciated texture and crumbled and contorted beds have been also documented in the type locality and type section. The author also recognized a mixture of shaly and marly materials with melikaria structures, and the upper and lower contact boundaries of the Barsarin formation are identified by the Chia Gara formation and Naokelekan formation, respectively.

The Barsarin formation was previously described as an isolated lagoonal environment, mostly evaporitic, and based on stratigraphic position from fieldworks, it has a Kimmeridgian-Tithonian age [13]. The recent studied section of the Barsarin formation is located close to the Rawandus area where the series is considered as the reference section in northeastern Iraq (Figure 1).

**Figure 1.** General map illustrates the location of the studied area "green circle".

#### **3. Methods**

Twenty-three samples were collected from the outcrop of the Barsarin formation. All the samples were prepared in laboratory for thin sections and studied under the optical microscope to distinguish the different carbonate phases, with a particular attention on the geometric cross-cut relations of fractures and veins. Scanning electron microscopy (SEM) was also used to image the surfaces to provide information about the morphology or texture of the surface. Selected samples of broken fragments (chips) mounted on aluminum stubs were studied from samples coated with gold depending on the purpose of work. SEM and EDX are also utilized in the recent study, backscattered image mode was used to give the different in contrast between minerals with different atomic number. The low atomic number samples give low emissions of backscattered electrons, while high atomic number samples give high emissions of these electrons. The backscattered electrons have higher energies than secondary electrons—usually from approximately 8 × <sup>10</sup>−<sup>18</sup> J (50 eV) up to the energy of the primary beam electrons.

The oxygen and carbon isotopic compositions of (24) samples were analyzed from powders, after selective microdrillings of the different recognized carbonate phases. However, the co-occurrence of fracture-filling calcite and dolomite, and also evaporites in the host limestones, made it difficult to avoid a mixing between these phases. As a consequence, the dolomite and calcite are considered as one group, and host limestone and evaporite as a second group.

Carbonate powders were reacted with 100% phosphoric acid at 70 ◦C using a Gasbench II connected to a Thermo Fisher Delta V Plus mass spectrometer. All values are reported in per mil relative to V-PDB (Table 1). Reproducibility and accuracy were monitored by replicate analysis of laboratory standards calibrated to international standards NBS19, NBS18 and LSVEC. Laboratory standards were calibrated by assigning δ13CVPDB values of

+1.95‰ to NBS19 and −46.6‰ to LSVEC and by assigning <sup>δ</sup>18OVPDB values of −2.20‰ to NBS19 and −23.2‰ to NBS18. The analyses (23) were performed in the University of Erlangen.(Germany, M. Joachimski).


**Table 1.** δ13C (% VPDB), δ18O (% VPDB) values of selected samples from Barsarin formation (n = 24).

Strontium isotope measurements on carbonate samples were performed by Laser Ablation ICP-MS at Goethe University-Frankfurt, using a Thermo-Finnigan Neptune multicollector (MC)-ICP-MS system attached to a Resolution 193 nm ArF Excimer laser ablation system (ComPexPro 102F, Coherent), equipped with an S-155 two-volume (Laurin Technic, Australia) ablation cell. Square laser spots with diameters or edge lengths of around 235 μm were drilled with an 8-Hz repetition rate, and energy density of about 6–7 J cm–2, during 45 s of data acquisition.

The collector set-up included (1) 83Kr as a monitor to verify the successful correction for the isobaric interferences of 84Kr on 84Sr and of 86Kr on 86Sr by subtraction the gas blank measured before sample ablation; (2) mass 83.5 to gauge the production rate of doubly charged 167Er and calculate the production of doubly charged 164Er, 166Er, 168Er and 170Er, which interfere with 83Kr and 84Sr, 85Rb and 86Sr, respectively; (3) mass 86.5 to gauge the production rate of doubly charged 173Yb and calculate the production of doubly charged 168Yb, 172Yb, 174Yb and 176Yb, which interfere with 84Sr, 86Sr, 87Sr and 88Sr; (4) 86Sr and 88Sr to use 88Sr/86Sr for mass bias correction; (5) 84Sr to check the accuracy of the mass bias and interference (see above) correction using 87Sr/86Sr; (6) 85Rb to correct for the isobaric interference of 87Rb on 87Sr and to get the 87Rb/86Sr; and (7) 87Sr to use in the radiogenic isotope ratio 87Sr/86Sr.

At the beginning of the analytical session, a soda-lime glass SRM-NIST 610 was measured to three times for empirical determination of 87Rb/85Rb mass bias. The procedure yielded, after interference correction on 86Sr, 88Sr, and 85Rb from doubly charged Yb and Er, the 87Rb/86Sr ratio needed for accurate correction of the isobaric interference of 87Rb on 87Sr. An isotopically homogeneous in-house plagioclase standard (MIR-1) was measured throughout the analytical session to monitor accuracy and apply corrections to the measured unknowns, if necessary. This plagioclase is a megacryst from a lava of the

Dutsin Miringa Hill volcano (Northern Cameroon Line; [14]) that has been independently characterized for 87Sr/86Sr.

The reference materials BHVO-1 and BCR-2G were measured along with the in-house standard to monitor accuracy, in particular of the correction of the isobaric interference of 87Rb, and also the reproducibility of the results. The corrections for the presence of Rb were minor and the results for the unknowns, reported in Table 2, are considered to be accurate. Sr concentrations compared with MIR (3500 ppm), applying the same corrections as for MIR. All analyses together are given in Table 2.

**Table 2.** Strontium concentration values and 87Sr/86Sr ratios by laser spots analyses using ICP-MS (B9, B14, n = 35).


#### **4. Results**

*4.1. Field Observation*

The base of the Barsarin formation is characterized by decimetre-thick beds (<30 cm) of stromatolites, evaporites (Figure 2B,C) and a prominent cherty level, alternating with wellbedded dolomitic limestones and shaly limestones (Figure 2A). The top of the formation shows massive limestones and dolomitic limestones (Figure 3). The thickness of the formation is 8 m in the studied Barsarin section type.

**Figure 2.** (**A**) Photograph showing the exposed section of the well-bedded limestone deposited during Upper Jurassic, Barsarin formation. (**B**) Photomicrograph illustrates the micro-laminated features of Bositra-like pelecypod shells (see [15]), the micro-laminated features cross-cut by vein of evaporites; (**C**) collapsed mudstone microfacies, giving auto-brecciated sediment. The pore spaces of the breccia are completely occluded by evaporates, with no traces of late diagenetic cements.


**Figure 3.** Master log showing the general lithology of the Barsarin Formation (Not to scale).

The fracture- and void-filling dolomite and calcite are distributed throughout the Barsarin profile, in addition to several vertical calcite filling joints and few horizontal ones in the upper part.

#### *4.2. Petrography*

The host carbonates show three types of microfacies: (i) mudstones, (ii) radiolarian wackestones and (iii) stromatolite boundstones (Figures 3 and 4). Despite these facies having undergone dolomitization, they are still recognizable, particularly the fine Bositra-like pelecypod shells and stromatolitic fabrics in the lower and parts of the section (Figure 2B).

**Figure 4.** Photomicrographs illustrating characteristics of silicate mineral content. (**A**) Litho-clast "right most part of the photo" embedded by mega-quartz grains. (**B**) Micro-organisms of spherical shapes of Radiolaria, the chalcedony micro-crystalline quartz partially filled the pore spaces of Radiolaria. (**C**) Mega-quartz grains in places replaced by dolomite minerals, see the strong relief and inclusions of quartz in the core of some grains. (**D**) Late diagenesis dolomite cements "D**b**" replaced the microcrystalline quartz.

Abundant needles/lath crystals from evaporite minerals are found within the mudstone microfacies, which shows an in-situ brecciated fabric with irregular cracks and veins (Figure 2C). Fracture-filling evaporites are volumetrically less important than the needle shaped evaporites (Figure 2A). The in-situ collapse of the pristine facies leading to a breccia is related to the dissolution of the evaporitic minerals (Figure 2C). Sometimes, these breccias are still filled with evaporite minerals.

Dolomite appears as a cement in two forms. D**<sup>m</sup>** represents the initial precipitate, and it is characterized by dark coloured crystals of variable sizes and anhedral-euhedral shapes (Figure 5B,C). D**<sup>m</sup>** shows, in places, a typical rhombohedral shape, with micriticrich core and transparent cortex (Figure 5B). Larger transparent euhedral crystals are also observed. The second cement occurs in the fractures and vugs and consists of a dirty coarse euhedral dolomite (D**b**) with curved or non-planar surface. It can also be euhedral with a characteristic zonation and shows a sweeping extinction under crossed nicols (Figure 5C,D). D**<sup>b</sup>** only occurs in fractures and open spaces. Fracture-filling calcite crystals are always associated with the second dolomite (D**b**); it displays variable sizes, is transparent and consists of blocky, twinned crystals. Geometric relationships under petrographic analysis show that the blocky calcites were replaced by coarse dolomite (D**b**). In addition, the detrital grains, mainly quartz, are observed within the host carbonate. They show the traces of dissolution and erosion, and this could be due to the distance of grain transportation (Figure 6A,D).

**Figure 5.** Photomicrographs illustrate: (**A**) The needle shape of evaporites within the mudstone microfacies that filled the pore spaces, (**B**) early dolomitizing cement (D**m**) were mostly preserved the previous pristine facies "dark brown colour", while in the central part of the photo the rhombohedral shape of dolomite which precipitated later. The lower part (**C**) and upper part (**D**) of the photos show the predating cement formation "D**m**" and the coarse grains of dolomite "Db" from the late dolomitizing fluids.

**Figure 6.** (**A**–**D**) Photomicrographs illustrate the weathered clasts of quartz from outside of the Barsarin basin, close up the eroded outline of polygonal grains (**B**) and other detrital grains in (**D**,**C**) that embedded in host carbonate rocks.

SEM analyses reveal that the host limestone, dolomite and calcite cements are commonly associated with quartz grains, more specifically at the bottom and upper parts of the formation (Figure 7). The grain size of quartz ranges from micro-quartz to mega-quartz, and in places silica appears in its chalcedony form.

#### *4.3. Oxygen, Carbon and In Situ Strontium Isotopes*

The O–C isotope compositions of the pristine facies, calcite, dolomite and evaporite are listed in Table 1. The 23 samples from fracture-filling dolomite and calcite cements, matrix dolomite (D**m**) and pristine facies of the formation have oxygen–carbon isotopic values populated into two groups. Group I is isotopically heavier in composition than group II (Figure 8). The host limestone and evaporite oxygen–carbon values fall into group I, while the calcite and dolomite cement oxygen-carbon values mostly fall into group II.

Most of the carbon and oxygen marine isotopic compositions of the Barsarin facies are significantly lower than those of seawater limestone during Late Jurassic time reported by Brand and Anderson [16,17]. Isotope compositions of Late Jurassic belemnites led these authors to define the isotopic composition of seawater during this time, with δ18OVPDB values ranging from 0.04‰ to −0.99‰ and <sup>δ</sup>13CVPDB values from 0.97‰ to 1.51‰. In group I of the Barsarin samples, <sup>δ</sup>18OVPDB values are between −5.09‰ to −0.34‰, and <sup>δ</sup>13CVPDB between −5.23‰ to −11.66‰. In group II, <sup>δ</sup>18OVPDB values are −12.87‰ to −5.90‰, and <sup>δ</sup>13CVPDB values are −8.55‰ to +1.47‰. However, a few host limestones that are free of evaporites have lighter oxygen isotopic values (up to −4.34‰) than those containing evaporites (up to −2.33‰, see the in Table 1).

**Figure 7.** SEM/EDX images of the evaporite carbonate rock sample obtained from freshly gold coated surface, the images illustrate: (**A**) the carbonate grains cemented by evaporitic minerals, the circle marks the aggregate grains of gypsum crystals, the rectangle marks two phases of minerals and this part, enlarged in (**B**,**C**), where the central part, composed of euhedral, eroded the surfaces of calcite crystals where cemented by evaporites. (**D**,**E**) Mineralogical mapping of the same sample location (**C**); (**D**,**E**) the majority of sodium element was distributed around the euhedral crystals of calcite.

Two samples (B9, B14, Table 2) have been analysed by laser ablation to obtain an absolute high-resolution radiogenic 87Sr/86Sr ratios on fracture-filling cements and the carbonate matrix. In-situ measurement of 87Sr/86Sr signatures in calcite and dolomite are provided by laser ablation to obtain the absolute radiogenic composition of many measurements on a micrometre-sized scale (MSS) within a single crystal (Figures 9–11). Laser spots analyses by ICP-MS measured several precise points on each crystal in two samples (B9, B14, n = 35). The plotted87Sr/86Sr ratios display two positive excursions (Figures 10 and 11). The first excursion is characterized by the highest radiogenic signatures, ranging from 0.70789 to 0.72859. Dolomite and calcite crystals fall into the second positive excursion, with strontium isotope ratios varying from 0.70721to 0.70772. They are less radiogenic than the host limestone values. Most of the laser spot measurements from the host limestone show a significant Sr radiogenic ratio, and mostly fall into the first positive excursions. The absolute strontium isotope data of the host limestone is radiogenically very high compared to the data from literature [9], and do not fit with the late Jurassic seawater strontium isotope curve, which varied between 0.70671 and 0.70713.

**Figure 8.** Oxygen and carbon isotope compositions in the Barsarin formation (n = 24). The values are populated on two groups: the first, group I (evaporitic mudstones and wackestones), has heavier values than the second, group II (calcite and dolomite cements).

**Figure 9.** The seawater Sr-isotope curve shows minima in the Late Triassic, Pliensbachian-Toarcian, Callovian-Oxfordian, Aptian-Albian, and Cenomanian-Santonian. The recent study data "blue plots" from our study are shown for comparison. The recent data "this study" populated around two positive excursions: the first host carbonate positive excursion, and the second diagenetic carbonate positive excursion. The host carbonate value is highly radiogenic, higher than the whole seawater Sr isotope curve, while the value of carbonate diagenesis samples is still high and less radiogenic than host carbonate samples (modified after Jones and Jenkyns [9]). OAE = Oceanic Anoxic Event.

**Figure 10.** The absolute radiogenic strontium isotope values v. the location of the spot analysis "sample No. B14". The location of analysis is represented by cubic blue colour where the ablated spot size used is 213 μm and depth of crater is ~20 μm.

**Figure 11.** The absolute radiogenic strontium isotope values v. the location of the spot analysis "Sample No. B9". The location of analysis is represented by cubic blue color where the ablated spot size used is 213 μm and depth of crater is ~20 μm.

#### **5. Interpretation and Discussion of Conventional and Non-Conventional Isotopic Signature Trends**

#### *5.1. δ13CVPDB-δ18OVPDB Isotopes*

The estimation of isotope values of marine and diagenetic sediments is useful in providing a better understanding of seawater compositional changes during deposition or post-deposition due to any chemical or biological alteration. The oxygen isotope reconstruction in Mesozoic marine sedimentary rocks shows a general decrease of δ18OVPDB values of the global seawater throughout the Late Jurassic [16–18]. Brand and Anderson proposed a past seawater <sup>δ</sup>18OVPDB from around +2.5‰ to −3‰ during Late Jurassic. Furthermore, the general trends of a warming through the Late Jurassic yields a paleo-temperature value in the range of 20–30 ◦C [19]. <sup>δ</sup>18OVPDB values between −5.09‰ and −0.34‰ of our host limestone (Group I; Table 1) do not fit within the previously documented oxygen isotope values of marine carbonate during Late Jurassic and is probably related to paleo temperature conditions.

Paragenetically, our host carbonate rocks from the Upper Jurassic have significant negative values of oxygen and carbon isotopic compositions, as low as −5.09‰ and −11.66‰, respectively. The host carbonates predated evaporite formation (see Figure 7 and Group I from Figure 8). Host carbonate samples that contain relicts of evaporites have <sup>δ</sup>13CVPDB <sup>δ</sup>18OVPDB values (from 0.34‰ to −2.33‰ and −6.25‰ to −9.30‰, respectively) lower than those of the previous reported marine domain for the same period [20], but higher than those of the pure host carbonates in our study. This trend of wide variation in oxygen and carbon isotope compositions from the host carbonates might be due to fluctuations of sea level, possibly related to climate change, at a local, regional scale (e.g., [21]), if not at a global scale [9], and/or could be linked to a secondary imprint on marine rock composition.

Floating angular clasts of host carbonates are cemented by evaporite crystals, suggesting the evaporative episode (mainly gypsum and halite) post-dated the host carbonates (Figures 2C and 7). This led to the formation of collapse breccias from the host carbonates, probably related to sea-level fall and climate change, as measured previously by the ratio of precipitation to evaporation (e.g., [22]). Hence, the large depletion in oxygen and carbon values, together with systematic occurrence of evaporates, suggest that dry and warm climate conditions drove the primary oxygen signal to negative values, where oxygen isotopic composition of seawater in the study basin increased during precipitation (e.g., [23]).

Evaporite minerals reflect the hypersaline evolution of the basin that led to heavier oxygen isotopic values than those of marine limestone (e.g., [24]) [25] in contrast to the carbon isotopic composition [26]. Nevertheless, our 18O depleted values in Group I are lighter than those of the oxygen isotopic composition of the marine curve. δ18OVPDB depleted values in host carbonates (Group 1) reveal that the former evaporative episode was affected by climatic changes during Late Jurassic, possibly through river input during a humid season [26]. Therefore, the depleted oxygen values from our study, combined with published data for the Late Jurassic-Early Cretaceous interval, suggest that global and gradual warmer palaeotemperatures occurred during this interval when the Barsarin formation was deposited. A similar scenario has previously been reported by Lécuyer and Allemand [27] and Vickers [28], from clumped isotopes and conventional oxygen isotope data, to interpretate the climatic change throughout Jurassic Cretaceous times. Furthermore, the covariance between oxygen and carbon isotopic compositions in marine carbonates (Figure 8) is probably linked to period boundaries with climate change [29,30].

The O–C isotopic values of second group (Group II) are lighter in their oxygen and heavier in their carbon isotope composition than those in Group I, although some of the carbon values overlap the marine signature. Group II has distinct and wide variations of oxygen–carbon isotopic compositions, compared to Group I. Group II is influenced by late diagenetic fluids that had a role in precipitation of calcite and dolomite cements (Figure 5). The δ18OVPDB values of diagenetic cements exhibit more depleted values than those not influenced by diagenetic fluid alteration (Figure 8). This could be due to the long-term

interaction between the host rock and the diagenetic fluids, promoted by fluids heated up to 80 ◦C or more (e.g., [31]). Therefore, the considerable and significant variation of δ18O–δ13C from Group II could be linked to the evolution of fluids (e.g., [31]), with slow precipitation rates under subsurface conditions (e.g., [32]), or rapid precipitation of dolomite and calcite due to abrupt change in pressure and temperature [8,33]. The extremely negative δ18O values in the coarse-sized dolomite and calcite crystals (Figure 8) highlight the interaction of the host limestones with diagenetic fluids (e.g., [8]) as indicated in Group II. Since the oxygen isotope values are extremely negative in deep diagenetic conditions, where the mechanism for dolomite and calcite precipitation is slow, hot circulation is the most probable source for samples belonging to Group II. Nevertheless, Group II is populated in two δ18O assemblages. One follows the J trend of meteoric water with extreme depleted δ13C values, and the second shows lighter oxygen isotope compositions with slight changes in δ13C values. These are quite consistent with the mixing of hot and meteoric fluids.

The carbon isotope values in Group II also exhibit a wide variability in their composition, although some values still possess host rock signatures (Group I: Figure 8). This could be related to the buffered system during deposition of diagenetic cements (e.g., dolomite and calcite) where the diagenetic fluid mixed with rock-derived carbon of the previous carbonates (Group I).

Considering the quartz grains in the host carbonate, the lower and upper parts of Barsarin formation are richer in silica than those of the middle part. Since the oxygen isotope composition of carbonates rich in silica is lighter than seawater, the chemistry of water facies rich in silica is not related to a marine source. Therefore, the origin of the quartz grains is not related to diagenetic evolution, as previously reported in a similar case by Bustillo [34], despite the outer rims of the chert or quartz having been often replaced by dolomite/calcite in a late diagenetic event (Figure 4D). This case has been reported also by Bustillo [35].

#### *5.2. Sources of Variation of Absolute Strontium Isotopes from the Late Jurassic and Onwards, Utilizing Non-Conventional Laser Ablation ICP-MS*

Seawater 87Sr/86Sr ratios varied homogeneously through geological time, thus making this isotope proxy an important chemostratigraphic tool [36,37]. The evolution of the seawater 87Sr/86Sr ratio curve during the Jurassic-Cretaceous is widely driven by variation in continental weathering rates (riverine) and predominant oceanic fluxes of non-radiogenic strontium sporadically during the Jurassic [5]. Silica hosted carbonates are another source that impacted the seawater radiogenic composition curve [38]. These sources of radiogenic strontium isotopes had a different effect on the seawater strontium isotope record (Figure 8, [5,9]). Our data provides, for the first time ever, a high-resolution strontium isotope composition from pristine facies and associated diagenetic carbonates using non-conventional direct laser ablation (Figures 10 and 11) in Upper Jurassic carbonate source rock.

The diagenetic imprint has to be precisely considered since the diagenetic fluids are impacted by the uptake or loss of Sr radiogenic isotope and Sr concentration. Therefore, the petrographical characteristics and paragenetic sequence for the entire section is critical before interpreting the strontium isotopic analysis. In the Barsarin formation early silicification and evaporite minerals formation is possibly syn-depositional [39] or occurred during a very early diagenetic event (e.g., [35]). The collapse and brecciation of the mudstone facies suggest that the host sediments and evaporite minerals were coincident with the depositional time of Barsarin formation. Silica infilling the intragranular porosity is also observed within dolomite and calcite lattice and suggests a late diagenesis process. The co-occurrence of silica and diagenetic carbonates bring the idea of silica recycling and/or continuous flux "input" of silica during and after the depositional setting.

The 87Sr/86Sr ratios are populated into two positive excursions. The first excursion gives a radiogenic isotope value of up to 0.72859 in the host carbonate rocks; these values are three times higher than those reported during Late Jurassic times [9]. The second excursion observed in the calcite and dolomite cements is less radiogenic than host carbonates. The core of carbonate cement shows impurities due to growth of quartz grains in their crystal framework (Figure 4C). These two positive excursions of absolute strontium radiogenic data could be explained in the following two ways:

(i) The silicate sediments that precipitated from seawater were converted into fibrous quartz and/or micro- (i.e., the biogenic silica-rich sediments) and mega-quartz via dissolution–reprecipitation processes [40]. The first excursion in pristine facies, that often contains silicate minerals in the studied samples, has lower Sr concentration and higher radiogenic strontium isotope values (high silicate minerals). The silicates were observed through the whole thickness of the carbonate facies, and as well as solid inclusions within the core of the calcite and dolomite crystal lattices. These inclusions could be the reason for the enrichment of the strontium ratio in the carbonate cements. Therefore, silica cycling via dissolution–reprecipitation processes could be one possible scenario for the observed radiogenic strontium, mainly for the first positive excursion (Figures 10 and 11). The silicate sediments in host carbonate show no corrosion or disturbed grains due to dissolution of the siliceous sediments or replacement, but in places the groundmasses are embedded with silicate grains totally different from the host carbonate facies, as highlighted frequently by their different colours from those of the facies (Figure 4). This indicates that some of the silicate sediments likely derive from outside of the Barsarin depositional basin, i.e., they are siliceous sediments continental detritus (Scenario II).

(ii) The second positive excursion relates to the diagenetic cement, where the strontium isotope ratio and strontium contents formed two assemblage groups: the highest strontium ratios with lowest strontium concentration (first group) and the lowest strontium ratios with highest strontium contents (second group). In addition, the weathered and eroded particles embedded in host carbonate (Figure 6A–D) reveal that this phase it could be derived from the mixing of detrital silicate minerals transported by rivers from the catchment "low Sr content" with high diagenetic Sr content "hot fluids".

The strontium isotopic composition of seawater is predominately controlled by inputs at mid-ocean ridges, through hydrothermal exchange (87Sr/86Sr = ~0.703), and from rivers, through continental erosion (87Sr/86Sr = 0.705 to > 0.800; mean ~0.712; Palmer and Edmond, 1989). Variations in seawater (87Sr/86Sr reflect changes in the concentration of Sr and the 87Sr/86Sr of dominant sources, [18]), with the hydrothermal source being less radiogenic (~0.703), generally account for less of the total Sr input in the ocean system [41]. The concentration and isotopic value of the riverine flux is highly variable due to heterogeneities in the isotopic composition of continental crust, the concentration of Sr in source rocks and weathering rates.

The increased input of hydrothermal strontium at oceanic ridges had a major effect on seawater strontium isotope composition during the Jurassic [11]. Hydrothermal sources decreased seawater 87Sr/86Sr ratios sporadically during the Bajocian–Callovian period, while, during mid-Bajocian, elevated submarine volcanism increased the 87Sr/86Sr ratios (cf. [36]; Figure 9). The enhanced hydrothermal strontium input culminated with the global sea-level rise at the Middle-Late Cretaceous [42]. However, more published data showing a degree of scatter confirm the uniform rise of seawater strontium isotopic curve from Late Jurassic–Early Cretaceous times [23]. During the Late Jurassic, major events, like early Mesozoic rifting and Jurassic subduction, have been suggested to occur during the Zagros Orogeny [43]. Stratigraphic and geochronological keys suggest that a large volume of detrital continental crust contaminated by mantle-derived magmas impacted Late Jurassic times [44]. These authors concluded that the strontium enrichment was either derived from an enriched mantle source or acquired by crustal assimilation. Although, by the Late Jurassic the highest 87Sr/86Sr ratio reported by Lechmann [44] was as large as 0.70697, which is very low compared to our study (values reaching 0.72859).

The widespread values along one host limestone sample's "first excursion" has a higher radiogenic isotope ratio, of as much as 0.72859 (Figure 10), although the ratio in the seawater reference curve did not exceed 0.70720 during Late Jurassic times (Figure 9). Davies and Smith [45] documented that the strontium isotope ratio increased relative to

pristine values by releasing strontium through the interaction of hydrothermal "hot" fluids with siliciclastic basement, or by release from continental siliciclastics [9]. The latter source is considered as radiogenic when evolved from riverine input [46]. Linking the continental siliciclastic source with our petrographical observations, the host carbonate rocks of the Barsarin formation show that silicification came from radiolarian cherts (Figure 4B). Yet, Radiolaria are not sufficiently abundant to be a sink for strontium isotopes. The early diagenetic siliceous cement is related to the source of silica from Radiolaria [47], and the primary intraparticle porosity was filled by chalcedony chert, as occasionally observed in our study area. Different-sized lithoclasts are embedded in the host carbonate rocks and contain variable sized-silicate grains (Figure 4A). These embedded silicate grains exclude that the radiogenic strontium was provided from the radiolarian cherts. Hence, at least the first radiogenic excursion was controlled by transport-limited nature of the weathering reaction, the same case was reported by Palmer and Edmond [41].

In SE France, hydrothermal dolomite, hosted in Upper-Jurassic sediments, were discovered with similar radiogenic values (ratios reaching 0.71120, Salih et al., 2020b). In this case the dolomite was related to radiogenic sources derived from riverine input. Furthermore, uplift in the latest Jurassic initiated detrital inputs related to the erosion along the Zagros Mountains [43]. Thus, the first positive excursion of an absolute strontium ratio would indicate an increase in weathering rates due to flux of radiogenic strontium.

The second positive excursion from the diagenetic cement samples evolved from diagenetic fluid sources. The dolomite and calcite cements preserved the inclusions and traces of silicate materials (Figure 4C). The remanent of solid inclusions inside the dolomite and calcite thereby limits the strontium isotope ratio. The remanent of silicate minerals either originated from host carbonates during replacement process or from interaction of hot fluids with radiogenic basement [5]. The origin of dolomite from hot fluids would preferentially record the interaction of hot fluids interaction with a radiogenic basement (cf., [5]), or the radiogenic source together with negative δ13C–δ18O linked to seepage of low-temperature meteoric water in subsurface carbonate rocks [48].

#### *5.3. Origin of δ13CVPDB δ18OVPDB-Light and 87Sr/86Sr-Rich Sediments*

The enrichment and depletion of strontium and oxygen–carbon isotopes in carbonates is mainly controlled by sea-level fluctuation, the erosion of continental sediment and the recycling of ocean water through mid-oceanic ridges [49,50]. The absolute 87Sr/86Sr ratios of the host rocks and diagenetic cements in Barsarin formation show two positive radiogenic excursions, which are higher than those reported in the Upper-Jurassic marine carbonates [9]. The two positive excursions that have a linear extrapolation with Sr concentration, suggesting that 87Sr/86Sr ratios could be released from meteoric (87Sr/86Sr ratios of up to 0.72859 and Sr concentration of 75 ppm) in the host carbonate rocks (Figure 12).

The deviation of our Sr isotopes values from the Sr isotope-seawater curve could be linked to sea-level fall during Late Jurassic, with a relatively warmer climate in the Kimmeridgian [51]. Sea-level fall is further supported in our study by the co-occurrence of evaporitic minerals floating inside the mudstone and in-situ brecciated mudstones. These observations are consistent with a first positive excursion of the strontium isotope ratio (up to 0.72859) and the heavier oxygen-isotope composition in those samples containing evaporites compared to those of host carbonates free of evaporites. These arguments suggest the beginning of sea-level fall and dry climatic conditions during the deposition of the Barsarin formation in Late-Jurassic times.

**Figure 12.** Strontium concentration v. 87Sr/86Sr ratios in host carbonate samples; note the variable data of Sr concentration and wide range of Sr isotope ratios.

This discussion can be further evaluated by oxygen and carbon isotopes and petrographic examination: if we considered that high Sr-isotope ratio and low Sr concentrations represent fresh water mixed with evaporated seawater in the host limestones, then δ13C and δ18O values would be depleted in their composition. This is the case in Group I, where oxygen-isotope composition shows a wide range, being heavier in the host limestones lacking evaporite. This wide range in Group I suggests the involvement of continental input by a fresh-water river source during the deposition of the Barsarin formation. As a consequence, δ13C and δ18O values from the host limestones that closely match the first positive excursion of strontium isotopes could be linked to a short-lived arid episode during Late Jurassic. A global warm climate will increase the rate of chemical weathering [18]. Therefore, warmer climates must increase both riverine Sr flux and its 87Sr/86Sr ratio [9]. This is well-supported by the dissolution and brecciation of the mudstones cemented by evaporites (Figure 2C).

The nature of the fluids involved in and after deposition of any carbonate is influenced by the composition of precipitated carbonates, therefore, other environmental conditions need to be considered. The late carbonate cements that contain silica showed relatively high radiogenic signatures (Figures 10 and 11) matching closely the radiogenic isotope values of hot fluids that interacted with the basement or with bedrocks rich in silicate minerals [5,45]. The Sr-isotope values of the dolomite/calcite cements fell into the second positive excursion of 87Sr/86Sr values (up to 0.70772) and could be linked in the Barsarin formation to the interaction of hot circulation fluids with the basement recrystallized rocks or interaction of hot fluids with evaporites or siliciclastic rocks (e.g., [52,53]). The negative values of oxygen isotope compositions from dolomite and calcite cements suggest a hot diagenetic fluid involvement [33]. However, the calcite precipitation from meteoric water could also be associated with an enrichment of strontium isotopes and the depletion of δ13C and δ18O values, but the negative covariant trend of δ13C and δ18O values exclude the involvement of meteoric water [30]. Furthermore, the lower limit of marine isotope concentration is about 200 ppm [54], while the strontium concentration from diagenetic carbonate is populated in two assemblage groups between 148–195 ppm and 362–1029 ppm (Figure 13). The lowest

Sr content, with highest strontium isotopes, indicate a meteoric-derived water, while the highest strontium content with lowest strontium isotopes indicates a possible influence of hot fluid in the diagenetic system. The co-occurrence of dolomite and calcite crystals better supports our scenario, in which dolomite and calcite are attributed to mixing of two different fluids. Therefore, the first positive excursion of direct 87Sr/86Sr ratio by laser ablation from the host limestones is partly linked to the weathering riverine strontium input during the depositional time of the Barsarin formation along the ZTFB. The second positive excursion of 87Sr/86Sr ratios is most likely related to post-depositional processes, involving two mixed sources of diagenetic fluids in the subsurface setting.

**Figure 13.** Strontium concentration v. 87Sr/86Sr ratios in diagenetic samples, the high strontium concentration is consistent with wide range of 87Sr/86Sr ratios, and the low strontium concentration is consistent with high 87Sr/86Sr ratios.

#### **6. Conclusions**

1. The origin and evolution of host carbonates and the diagenetic fluid history from Barsarin formation were identified on the basis of fieldwork, petrography and geochemistry (δ13C–δ18O and high resolution 87Sr/86Sr isotope ratios by laser ablation on a micrometersized scale -MSS).

2. The Barsarin carbonate sequence contains silica infilling the intragranular porosity. Silica is also very abundant, as are inclusions in the dolomite and calcite cements that filled the fractures and pore spaces that affected the formation. The presence of silica indicates that simultaneous diagenetic processes affected the host rocks and cement formation.

3. δ13C–δ18O isotope values were populated into two groups. The host carbonates of Group I had isotopically heavy oxygen compositions compared to the diagenetic cements of Group II. Group II showed two different assemblages of oxygen-carbon isotope values, suggesting that mixed diagenetic fluids were involved in the precipitation of the carbonate cements.

4. Two positive excursions of the strontium isotope-ratio curve have been identified. The first excursion is related the host carbonates, and the second is the fracture-filling carbonate cements.

5. Two scenarios were proposed for the first positive excursion of the strontium ratios. The first suggests that a silica cycling occurred through dissolution-reprecipitation processes, while the second scenario linked the high values of strontium ratio to weathering of detrital silicate rocks in the headwaters of a fresh river mixed with evaporated seawater. The carbon and oxygen isotopes, and the two assemblage groups of the first excursion of strontium isotopes and strontium contents exclude the first scenario and indicate that the host carbonates formed in a marine evaporative environment mixed with fresh water originated from riverine-rich silicate detrital "Second scenario".

6. The closely associated silica contents with the first positive excursion of absolute 87Sr/86Sr ratio in the host carbonates are related to the continuous flux of riverine input during Late Jurassic.

7. Silica and evaporites with low Sr concentration, high 87Sr/86Sr ratios and wide range of oxygen and carbon isotope compositions point to a riverine weathering input in the evaporated seawater during the sea level fall which occurred during the Scenario II event.

8. The wide range of δ13C–δ18O values in Group II and their positive covariance suggest a late fluid mixing during the second excursion of the 87Sr/86Sr ratio. In this case mixed fluids, consisting probably of low-temperature meteoric and hot diagenetic fluids, were inferred from the stable isotope compositions of the post-depositional carbonates "calcite and dolomite cements".

**Author Contributions:** Writing and preparation the original draft, N.S. and A.P.; review and editing, K.K., J.-N.P., N.S. and A.P., Methodology and software, A.G., A.P. and N.S. All authors have read and agreed to the published version of the manuscript.

**Funding:** The study benefited from research funds of the Université Libre de Bruxelles (ULB)- Belgium.

**Acknowledgments:** This work was supported by Université Libre de Bruxelles-Belgium and University of Alberta-Canada.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


## *Article* **Dolomitization of Paleozoic Successions, Huron Domain of Southern Ontario, Canada: Fluid Flow and Dolomite Evolution**

**Ihsan S. Al-Aasm 1,\*, Richard Crowe <sup>2</sup> and Marco Tortola <sup>1</sup>**


**Abstract:** Integrated petrographic, isotopic, fluid inclusion microthermometry, and geochemical analyses of Paleozoic carbonate successions from multiple boreholes within the Huron Domain, southern Ontario were conducted to characterize the diagenetic history and fluid composition, on a regional scale, and evaluate the nature and origin of dolomitized beds. Multiple generations of non-stochiometric dolomite have been observed. These dolomites occur as both replacement (D1 and D2) and cement (saddle dolomite; SD) and formed either at near-surface to shallow burial zone (D1) or intermediate burial (D2 and SD). Petrographic and geochemical data of dolomite types and calcite cement suggest that these carbonates have experienced multiple fluid events that affected dolomite formation and other diagenetic processes. Cambrian and Ordovician strata have two possibly isolated diagenetic fluid systems; an earlier fluid system that is characterized by a pronounced negative shift in oxygen and carbon isotopic composition, more radiogenic Sr ratios, warm and saline signatures, higher average ∑REE compared to warm water marine brachiopods, negative La anomaly, and positive Ce anomaly; and a later Ordovician system, characterized by less negative shifts in oxygen and carbon isotopes, comparable Th, hypersaline, a less radiogenic, less negative La anomaly, and primarily positive Ce anomaly but also higher average ∑REE compared to warm water marine brachiopods. Ordovician, Silurian, and Devonian Sr isotopic ratios, however, show seawater composition of their respective age as the primary source of diagenetic fluids with minor rock/water interactions. In contrast, the isotopic data of the overlying Silurian and Devonian carbonates show overlaps between δ13C and δ18O values. However, δ18O values show evidence of dolomite recrystallization. D2 shows wide Th values and medium to high salinity values. Higher Th and salinity are observed in SD in the Silurian carbonates, which suggest the involvement of localized fluxes of hydrothermal fluids during its formation during Paleozoic orogenesis. Geochemical proxies suggest that in both age groups the diagenetic fluids were originally of coeval seawater composition, subsequently modified via water-rock interaction possibly related to brines, which were modified by the dissolution of Silurian evaporites from the Salina series. The integration of the obtained data in the present study demonstrates the linkage between fluid flux history, fluid compartmentalization, and related diagenesis during the regional tectonic evolution of the Michigan Basin.

**Keywords:** fluid flow; dolomitization; Paleozoic; Huron Domain

#### **1. Introduction**

As part of an ongoing investigation into the occurrence of strata-bound dolomite layers within sedimentary, Paleozoic-age formations of southern Ontario an initial study using core acquired from the site of a formerly proposed deep geological repository (Bruce nuclear site) was conducted [1]. This initial study provided a site-specific analog to examine fluid migration and rock-formation barrier integrity at geological timescales. Results showed that hydrothermal alterations were minimal, burial-related, and lacked isotopic signatures associated with fault-controlled diagenesis, as observed in other regions

**Citation:** Al-Aasm, I.S.; Crowe, R.; Tortola, M. Dolomitization of Paleozoic Successions, Huron Domain of Southern Ontario, Canada: Fluid Flow and Dolomite Evolution. *Water* **2021**, *13*, 2449. https:// doi.org/10.3390/w13172449

Received: 3 August 2021 Accepted: 30 August 2021 Published: 6 September 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

of southern Ontario. This program was then expanded to a regional scale, to determine if the conditions observed at the Bruce nuclear site were consistent across the Huron Domain of southern Ontario. Core samples from multiple deep boreholes within the Huron Domain of southern Ontario were analyzed for petrographic, stable, and Sr isotopic composition, fluid inclusion microthermometry, and major, trace, and REE to characterize diagenetic history, fluid composition, and sedimentary provenance (e.g., [2,3]).

Extensive fluid flow occurs during tectonic thrusting, sedimentary loading, uplift, and compression in many sedimentary basins [4]. These fluids can affect the redistribution of mass, heat, petroleum migration, and formation of sediment-hosted ore deposits [5–7]. A fundamental problem in geology is to constrain the composition of the fluids within the context of thermal and tectonic evolution (i.e., changes of stress regimes) of hosting sedimentary basins. The composition and evolution of ancient sedimentary fluids have been successfully reconstructed using recent advances in stable and radiogenic isotopes, trace and rare-earth element geochemistry, fluid inclusion analyses, and paleomagnetic techniques (e.g., [7–13]). The fluids identified in these studies ranged in composition from evaporative brines, marine, mixed marine-meteoric, to deep brines and hydrothermal waters.

Dong et al. [14] showed that strata-bound dolostones form in a broad range of geologic and geochemical processes and generally result from large-scale fluid flow. Important aspects of dolomitization that need further thorough investigation are: (1) the source and supply of Mg [15,16]; (2) the plumbing mechanism necessary to transport Mg-bearing solutions to the site of dolomitization [17]; (3) the stoichiometry and recrystallization of dolomite [18–23], and (4) the effect of tectonics on the fluid flow [5,24–27]. An understanding of the genesis of ancient dolomites is only possible through detailed petrographic and geochemical studies of dolomitized sequences where the geologic history of the rocks is also reasonably well constrained.

There are several models suggested explaining the formation of dolomite over a broad range of geologic and geochemical processes. These models include [28–31]: (1) sabkha dolomitization model, (2) seepage-reflux dolomitization model, (3) mixing zone dolomitization model, (4) seawater dolomitization model (5) burial-dolomitization model, and (6) structurally-controlled hydrothermal dolomitization model.

Geological and geochemical processes lead to a change in formation density that can alter the physical properties of the bedrock, including increased porosity and permeability. In this sense, dolomitization can influence groundwater system evolution through the creation of permeable pathways following lithification, especially along localized geologic structures suitable for the accumulation of hydrocarbons. These models are dependent on three criteria [30]: (1) Thermodynamic: the fluids must be supersaturated with respect to dolomite and undersaturated with respect to calcite. Thus, dolomite cement will precipitate, and calcite cement will dissolve. (2) Kinetic: the rate of dolomite precipitation must be equal to or greater than the rate of calcite dissolution. (3) Hydrologic: Mg-rich fluids must be pumped continuously in the limestone pore system, faults, and fractures [30].

Recent studies by Haeri–Ardakani et al. [32,33] identified three distinct dolomite types within southwestern Ontario: a microcrystalline dolomite (D1) formed during early diagenesis from Middle Ordovician seawater and recrystallized during progressive burial, a medium crystalline dolomite (D2) formed by hydrothermal fluids (68 to 99 ◦C) and a late-stage saddle dolomite cement (D3) related to fault controlled, high temperature fluids (>125 ◦C).

Most recently, Al-Aasm [1], Al-Aasm and Crowe [2] and Tortola et al. [3] assessed the nature and origin of strata bound dolomitization occurring in rock samples from Cambrian through Devonian carbonates and siliciclastics beneath the Bruce nuclear site and adjacent areas within southern Ontario. These latter studies better supported the notion of fluid compartmentalization between different stratigraphic successions.

Hence, the main purpose of this study is to thoroughly examine all data gathered in both previous and new studies by the authors to further evaluate, on a regional scale, the

nature and origin of dolomitized beds within Paleozoic carbonates of the Huron Domain with the intention on deciphering the nature of dolomitization processes, evolution of diagenetic fluids, the nature of paleohydrogeological system(s) and whether they are connected on not, the source of heat for hydrothermal fluids (if any) and the driving force for diagenetic fluid(s) migration, in contrast to previous studies in southwestern Ontario.

#### **2. Regional Geologic, Tectonics, and Stratigraphic Framework**

The study area is located on the northeastern margin of the Michigan Basin (Figures 1 and 2). It represents a portion of the northwestern flank of the Algonquin Arch that consists of Paleozoic sediments overlying the subsurface basement [34]. The thickness of the Paleozoic sequence is estimated to range from a maximum of 4800 m at the center of the Michigan Basin to 850 m on the flank of the Algonquin Arch.

**Figure 1.** (**A**) Generalized Paleozoic bedrock geology of southern Ontario (modified from [35,36]); (**B**) Tectonic history of the Phanerozoic successions in southern Ontario (modified from [37]).

**Figure 2.** Bedrock geology of southern Ontario showing the location of boreholes used for this study within the Bruce Nuclear site and other areas in southern Ontario (modified from [38]).

Southern Ontario separates two major Paleozoic sedimentary basins, the Appalachian Basin to the southeast and Michigan Basin to the west (Figure 1A). The Algonquin Arch, which separates the two basins, extends from northeast of the Canadian shield to the southwest near the city of Chatham. Near the Windsor area, another structural high named the Findlay Arch begins and extends to the southwest into Michigan and Ohio states [39]. The Michigan basin is a nearly circular intracratonic basin dominated by carbonate and evaporite sediments, whereas the Appalachian Basin is an elongate foreland basin developed because of the collisional tectonic event along the eastern margin of the North American continent during the Paleozoic and is mainly characterized by deposition of siliciclastic sediments [40].

During Precambrian to Cambrian, a rifting event that separated North America and Africa initiated subsidence (approximately 1100 Ma) and deposition within the Michigan Basin (Figure 1B; [37]). Thermal subsidence of the Precambrian basement followed approximately 580 Ma to 500 Ma after the cooled lithosphere increased in density [41]. The subsequent and continuous deposition of sediment into the Michigan Basin caused the bending of the basement [42–44]. However, the development of the basin is not the result of a continuous subsidence but represents the outcome of a series of tectonic events that characterized the Paleozoic [45,46].

The Taconic and Acadian Orogenies (Figure 1B) had a dominant control on the Paleozoic strata of the eastern flank of Michigan [47]. In fact, during the Taconic Orogeny (Ordovician), a large-scale eastward tilting of the Laurentian margin has been produced with the consequent disappearance of the concentric depositional geometry of the Michigan Basin and with the subsequent deposition of Ordovician limestone and overlying shale successions [46,48,49].

The Middle Devonian Acadian Orogeny (Figure 1B) determined the reactivation of the Algonquin-Findlay Arch system with consequent return of the concentric geometry of the Michigan Basin [45]. During the Upper Silurian where restricted marine conditions dominated, the deposition was mainly characterized by the evaporites of the Salina Group [50,51]. In addition, the presence of an unconformity reveals emergent conditions at the end of the Silurian [47].

The Caledonian (Devonian) and Alleghenian (Carboniferous) Orogenies (Figure 1B) played a major role in terms of diagenetic fluid migration [5]. The regional maximum principal stress at that time was northwesterly oriented in accordance with the main direction of the thrust motion along the Appalachian tectonic front [47,52].

The Paleozoic sequence is separated from underlying Precambrian rocks of the ca. 1.1 Ga old Grenville orogeny by an unconformity overlain by a relatively thin Cambrian sandstone, which represents an aquifer [53,54]. The depositional environment of the Ordovician and Cambrian formations encompasses various facies assemblages, including supratidal, lagoonal to subtidal environments (e.g., [55,56]). The diverse facies changes observed in these formations were controlled by relative sea-level changes and sources of detritus, related to the regional tectonic setting during the Paleozoic [37,46]. In the study area, many of the Cambrian units along the Algonquin Arch (Figure 1) were absent due to arch rejuvenation and uplift during early Paleozoic times (e.g., [56]). Ordovician formations consist of thin, decimeter-scale, alternating shale, and carbonate beds; some of which have been dolomitized [47]. The Michigan Basin has been characterized by a typical tropical climate during mid to late Silurian because it was located at 25◦ south of the equator [57] leading to the development of shallow waters reef complexes. Moreover, the Michigan Basin has been isolated and is becoming an evaporative basin in which alternation of evaporites and carbonates is deposited during the rest of the Silurian [58]. Conversely, in the Middle Devonian, the Michigan basin has been characterized by an extremely arid climate with consequent deposition of shallow intertidal and supratidal lithofacies such as Amherstburg and Lucas formations [59]. The Devonian sequences are represented by limestones, shales, dolostones, and sandstones [39] with arid climate and deposition of shallow intertidal and supratidal lithofacies (e.g., Amherstburg and Lucas formations).

#### **3. Sampling and Analytical Methods**

This section may be divided into subheadings. It should provide a concise and precise description of the experimental results, their interpretation, as well as the experimental conclusions that can be drawn.

Cores (n = 91) from specific sections from deep boreholes (n = 10) from Bruce Nuclear Site (Tiverton, Ontario) and across the Huron Domain of southern Ontario were sampled (Figure 2). These samples encompass formations ranging from Cambrian to Devonian age successions and representing a wide range of rocks including dolomitized limestones, dolostones, sandy dolostones, and sandstones (Figure 2). These core samples were examined and photographed prior to their selection for making thin sections and fluid inclusion wafers. Petrographic examination of 145 thin sections was performed under a standard petrographic microscope to compliment core descriptions and by cathodoluminescence microscopy (CL) using a Technosyn cold cathodoluminescence stage with a 12–15 kV beam and a current intensity of 0.42–0.43 mA.

Samples for oxygen and carbon isotopic analyses were extracted from calcite and dolomite samples (n = 259) from polished slabs using a microscope-mounted drill assembly. The samples were reacted in a vacuum with 100% pure phosphoric acid for four hours at 25◦ and 50 ◦C for calcite and dolomite, respectively. CO2 gas from samples that contained a mixture of calcite and dolomite was extracted using the chemical separation method of Al-Aasm et al. [60]. The evolved CO2 gas was analyzed for isotopic ratios on a Delta Plus mass spectrometer. Delta (δ) values for oxygen and carbon are reported per mil (‰) relative to the PeeDee Belemnite (VPDB) standard. Precision was better than 0.05‰ for both δ18O and δ13C.

Twenty-six core samples were selected for Sr isotopic analyses from various calcitic and dolomitic components representing host rock, skeletal components, and diagenetic phases. Prior to 87Sr/86Sr analysis, powdered samples (n = 78) were extracted employing a microscope-mounted drill assembly to extract 1–10 mg of powdered calcite or dolomite from polished slabs. Care was taken to avoid contamination from other components such as later cement fill. These samples were initially dissolved in 2.5 N HCl suprapure for 25 h in sealed PFA vessels at room temperature. They were then passed through cation exchange quartz columns. To remove the eventual pore salts that may contain Sr, these samples were, after crushing and pulverization were washed twice with MilliQ water. Strontium isotopic ratios were measured on a Finnigan MAT 261 mass spectrometer. All analyses were performed in the static multi-collector mode using Re filaments. NBS and ocean water were used as standard references and 87Sr/86Sr ratios were normalized to 87Sr/86Sr = 8.375209. The mean standard error was 0.00003 for NBS-987.

Thirty-one samples were analyzed for fluid inclusion microthermometric studies. Petrographic work was done with a conventional petrographic microscope. Careful and detailed observations of size, occurrence habit, and paragenetic contexts of the fluid inclusions were conducted in the double polished wafers. Distinctions were made between primary, secondary, and pseudo-secondary fluid inclusions. This was followed by a determination of the number of phases (i.e., the liquid to vapor ratio) and their ratios entrapped in each of the fluid inclusions. Then determination of the composition (aqueous or hydrocarbon fluids) of the inclusion using fluorescence microscopy; most oil-filled inclusions fluoresce, while the aqueous inclusions do not. Finally, microthermometric measurements of the homogenization (Th) and final ice melting temperatures (Tmice) of fluid inclusions were measured. These measurements were conducted using a Linkam THMGS600 heating-freezing stage, which was calibrated using synthetic fluid inclusions of known compositions. Homogenization temperatures (Th), and ice-melting temperatures (Tm-ice) were measured with a precision (reproducibility) of ±1 ◦C and ±0.1 ◦C, respectively. The salinities for inclusions in which ice was the last melting phase were calculated using the program of Chi and Ni [61], and those for inclusions with hydrohalite as the last melting phase was done with the program of Steele-MacInnis et al. [62].

Deformation of inclusion walls, necking-down of single inclusions into smaller inclusions, or leakage of fluid can all lead to changes in homogenization temperature. To limit the possibility of measuring deformed aqueous inclusions, only primary inclusions from the same field of view were measured during a single heating or freezing run. By restricting measurements to inclusions within the same field of view, any sudden changes in liquid/vapor ratios due to inclusion deformation could be observed and removed from consideration. This does not exclude the possibility of subtle deformation in the z-axis, which is difficult to detect. Heating runs were conducted before freezing runs to reduce the possibility of inclusion stretching by freezing [63].

Major (Ca, Mg), minor, trace elements (Sr, Na, Mn, and Fe) and rare-earth element (REE) analyses (n = 95) were performed using ICP-MS. The procedure used for the sample digestion has been adapted to digest carbonates only. Around 50–100 mg of powdered carbonate samples from different calcite and dolomite phases were digested utilizing 5% v/v acetic acid. Measurements were carried employing an Agilent 8800 inductively coupled mass spectrometer at the University of Waterloo. All results we reported in part per million (ppm).

For rare-earth (REE) element analysis, samples were analyzed using external calibration, internal standard, and standard addition methods, following the procedure described in Jenner et al. [64]. Measured REE values were normalized to Post Archean Australian Shale (PAAS). Rare-earth element anomalies such as CeSN [(Ce/Ce\*)SN = CeSN/(0.5LaSN + 0.5PrSN)], LaSN [(Pr/Pr\*)SN = PrSN/0.5CeSN + 0.5NdSN)], EuSN [(Eu/Eu\*)SN = EuSN/(0.5SmSN + 0.5GdSN)], and GdSN [(Gd/Gd\*)SN = GdSN/(0.33SmSN + 0.67TbSN)] were calculated employing the equations proposed by Bau and Dulski [65]. The proportions of LREE(La-Nd) over HREE (Ho-Lu) were calculated via (La/Yb)SN ratios as proposed by Kucera et al. [66]. The proportions of MREE (Sm-Dy) was calculated using (Sm/La)SN and (Sm/Yb)SN ratios [67] over LREE and HREE, respectively.

An environmental scanning electron microscope (ESEM) coupled with an EDAX detector was utilized to analyze selected carbon-coated thin section samples to investigate the nature of dolomitization and subsequent recrystallization.

#### **4. Results**

#### *4.1. Petrographic Analysis*

Cambrian strata are dominated by mixed carbonate-siliciclastic lithofacies (Figure 3). The carbonate lithofacies is composed of dolomitized and silicified crinoidal and bryozoan components. The siliciclastic facies consist of silicified ooids and other coated grains as well as detrital carbonates with large detrital quartz, cemented by carbonates.

The Middle Ordovician carbonates in the subsurface are divided into the Black River Group and the overlying Trenton Group ([47]; Figure 3). The original depositional facies of the investigated Ordovician formations (Shadow Lake, Gull River, Coboconk, Kirkfield, Sherman Fall, Cobourg, and Collingwood) consists of nodular, bioturbated, micritic, and bioclastic facies—now partly dolomitized. In the Coboconk Formation, a thin fossiliferous layer (grainstone facies) occurs with a sharp contact with the nodular facies. This layer is also pervasively dolomitized. Skeletal components include bryozoans and echinoderms. In the Gull River Formation, the matrix is highly dolomitized but fossils such as brachiopods, crinoids, and trilobites are dolomitized but their textural integrity is preserved.

The Silurian formations in the Huron Domain are characterized by several lithofacies of partially to fully dolomitized carbonates, such as dolommudstone with abundant fenestrae, sometimes brecciated, and coated and bioclastic grainstone lithofacies. These carbonate lithofacies alternate with evaporite units, such as A-1, A-2, and B (Figure 3).

**Figure 3.** Paleozoic stratigraphy of the studied successions also showing the listing of the studied cores (modified from AECOM [47]).

The Devonian carbonates contain nodular and bioturbated lithofacies, partially dolomitized (e.g., Lucas Formation; Figure 3), algal laminated dolomudstone, dolomitized bioclastic grainstone, undolomitized bioclastic grainstone, and brecciated, cherty dolostone lithofacies.

The diagenetic history of the Cambrian, Ordovician, Silurian, and Devonian successions in the Huron Domain encompasses multiple episodes of calcite cementation, dolomitization, anhydrite cementation, mechanical and chemical compaction, dissolution, and silicification. It is understood (e.g., [2,3]) that these diagenetic processes occurred during deposition and continued upon burial. The paragenetic sequence shown in Figure 4 is based on petrographic relationships, spatial distribution, and geochemical evidence. Due to the significance of dolomitization in these carbonates, the discussion in this paper will be concentrated on the most significant and relevant diagenetic events.


**Figure 4.** Paragenesis of the Cambrian-Ordovician formations (modified from Al-Aasm and Crowe [2]) and paragenesis of the Silurian-Devonian formations (modified from Tortola et al. [3]).

#### 4.1.1. Dolomite Types

In this study, multiple generations of dolomite have been observed (Figure 4).

The dolomite types occur as both replacements (D) and cement (SD) (cf. [2,3]). According to Sibley and Gregg's classification [68] and based on petrographic characteristics (e.g., crystal size, shape, extinction, fabric, and paragenetic relationships), several types of dolomite development are distinguished. Table 1 summarizes the main types of dolomite and their petrographic and CL characteristics.

**Table 1.** Summary of the main petrographic characteristics of the main diagenetic minerals in the Paleozoic successions.



**Table 1.** *Cont.*

Cambrian:

In the mixed carbonate-siliciclastic lithofacies of the Cambrian strata, the presence of dolomite is restricted to the dolomitized bioclastic lithofacies and scattered dolomite crystals within the silicified lithofacies (Figure 5). Two replacive dolomite types occur in the Cambrian rocks (microcrystalline subhedral-D1 and medium to coarse crystalline subhedral to anhedral, zoned- D2; Figure 5) and minor amounts of saddle dolomite (SD) cement present in vugs (Table 1; Figure 5). This dolomite cement is succeeded by large blocky calcite cement (Figure 5I) in fractures and vugs. D2 mimically replaces allochems, such as fossils and ooids.

**Figure 5.** (**A**) Photomicrograph (XPL) of fracture-filling coarse saddle dolomite cement (SD) postdating anhedral, replacive medium to coarse crystalline dolomite matrix (D2). Sample: DGR2-CR279- 840.20 Cambrian. (**B**) Photomicrograph (PPL) of fracture-filling coarse saddle dolomite cement (SD) postdating planar-s, pervasive replacive medium to coarse crystalline (D2) and fine crystalline (D1) dolomite matrix. Sample: DGR2-CR139-2-844.87 Cambrian. (**C**) Photomicrograph (XPL) of fracture-filling coarse saddle dolomite cement (SD) postdating anhedral, replacive medium to coarse crystalline dolomite matrix (D2). Sample: DGR2-CR139-2-844.87 Cambrian. (**D**) Photomicrograph (CL) of fracture-filling, bright luminescent blocky calcite cement (BKC) postdating non-luminescent to red luminescent, zoned medium to coarse crystalline dolomite matrix (D2). Sample: DGR4-CR274- 2-844.39 Cambrian. (**E**) Photomicrograph (CL) of vug-filling, bright luminescent blocky calcite cement (BKC) postdating non-luminescent saddle dolomite cement (SD) and replacive medium crystalline dolomite matrix (D2) showing dull to red luminescence. Sample: DGR2-CR139-2-844.87 Cambrian. (**F**) Photomicrograph (UV-light) of planar-e pervasive replacive medium to coarse crystalline dolomite matrix (D2) showing dull fluorescent cores and bright rims. Sample: DGR4-CR274-2-844.39 Cambrian. (**G**) Photomicrograph (PPL) of fracture-filling coarse saddle dolomite cement (SD) in an oolitic sandstone. Sample: DGR2-CR134-1-844.54 Cambrian. (**H**) Photomicrograph (CL) of medium to coarse crystalline dolomite matrix (D2) partially replacing coated grains in an oolitic sandstone. Sample: DGR2-CR134-2-844.54 Cambrian. (**I**) Photomicrograph (PPL) of fracture-filling coarse blocky calcite cement (BKC) postdating planar-s, pervasive replacive medium to coarse crystalline (D2) and fine crystalline (D1) dolomite matrix. Sample: DGR2-CR133-1-843.96 Cambrian.

Ordovician:

Two types of replacive dolomite are recognized in the Gull River, Coboconk, and Cobourg formations (Table 1; Figure 6): (1) microcrystalline dolomite replacing skeletal components and matrix (D1), and (2) medium crystalline replacive matrix dolomite (D2), which is fabric destructive, partly zoned and replaces the intranodular micritic matrix (Figure 6). D2 dolomite is volumetrically more abundant than D1 and it is also observed along dissolution seams in partially dolomitized limestones, usually characterized by cloudy-dark cores and clear rims (Figure 6C). A minor presence of saddle dolomite (SD) occludes pores and fractures and postdates D1 and D2 (Figure 6D,E). Anhydrite cement postdates saddle dolomite (Figure 6E,F). In other Ordovician formations, such as the Shadow Lake, Kirkfield and Sherman Fall, Silurian, which are mostly dominated by limestones of varied lithologies, no significant dolomitization is observed. However, the scattered occurrence of D1 dolomite is observed in Sherman Fall Formation replacing the calcitic muddy matrix.

**Figure 6.** (**A**) Photomicrograph (PPL) of fracture-filling coarse saddle dolomite cement (SD) postdating planar-e to planar-s, replacive medium to coarse crystalline(D2) dolomite matrix. Sample: DGR3-CR279-2-839.80 Gull River. (**B**) Photomicrograph (PPL) of fracture-filling coarse saddle dolomite cement (SD) postdating, pervasive replacive micro to fine crystalline dolomite matrix (D1). Sample: DGR3-CR279-2-839.80 Gull River. (**C**) Photomicrograph (XPL) of planar-e replacive medium to coarse crystalline dolomite matrix (D2) and fine crystalline dolomite matrix (D1). Sample: T001925-1-890 (1-5) Cobourg. (**D**) Photomicrograph (XPL) of fracture-filling coarse saddle dolomite cement (SD) postdating planar-s, pervasive, replacive medium to coarse (D2) and fine crystalline (D1) dolomite matrix. Sample: DGR4-CR265-817.12 Gull River. (**E**,**F**) Paired photomicrographs (PPL and CL, respectively) of anhydrite cement postdating zoned, dull (rims) to dark red (cores) luminescent, coarse saddle dolomite cement (SD). Sample: DGR3-CR279-2-839.80 Gull River.

#### Silurian:

Dolomite in the Silurian succession occurs mostly in the Guelph, Salina (A1-carbonates, A2-carbonates, Figure 3), and Bass Island formations. Two types of replacive matrix dolomite and one of dolomite cement are observed (Table 1; Figure 7). These are: D1 which is a pervasive replacive micro to fine crystalline matrix dolomite and D2, which consists of medium crystalline, subhedral, zoned dolomite (Figure 7B,C). It also occurs as a less abundant selective replacive medium (>50–100 μm) euhedral crystalline dolomite, only observed along dissolution seams in partially dolomitized limestones and usually

characterized by cloudy-dark cores and clear rims (Figure 7C). Ferroan saddle dolomite (SD; Figure 7D–F) occludes fractures and vugs predating late blocky calcite cement (Figure 7E,F).

**Figure 7.** (**A**) Photomicrograph (PPL) of planar-s replacive micro to fine crystalline dolomite matrix (D1). The yellow arrow indicates the presence of framboidal pyrite in the matrix. Sample: DGR3-CR61-194.02 (3-11) Salina G-Unit. (**B**) Photomicrograph (PPL) of planar-s fine crystalline dolomite matrix (D1) and planar-e to planar-s, replacive medium crystalline dolomite matrix (D2), showing increasing intercrystalline porosity. Sample: DGR4-CR100-327.98 (4-11) Salina A1-Unit (Carbonate). (**C**) Photomicrograph (PPL—after staining) of planar-s fine crystalline dolomite and planar-e to planar-s, pervasive replacive medium crystalline dolomite matrix (D2) showing increasing intercrystalline porosity and non-ferroan cores followed by later iron-rich rims. Sample: DGR3-CR128-388.45 (3-16) Guelph. (**D**) Photomicrograph (PPL—after staining) of fracture-filling blocky calcite cement (BKC) postdating fracture-filling coarse saddle dolomite cement (SD) and fine crystalline (D1) dolomite matrix. Sample: DGR8-CR139-382.40 (8-13) Guelph. (**E**) Photomicrograph (CL) of fracture-filling, bright luminescent blocky calcite cement (BKC) postdating dull to non-luminescent coarse saddle dolomite cement (SD) and dark red luminescent fine crystalline dolomite matrix (D1). Sample: DGR8-CR139-382.40 (8-13) Guelph. (**F**) Photomicrograph (XPL) of fracture-filling blocky calcite cement (BKC) postdating fracture-filling coarse saddle dolomite cement (SD), replacive medium to coarse crystalline dolomite matrix (D2) and fine crystalline dolomite matrix (D1). Sample: DGR8-CR116-313.31 (8-9) Salina A2-Unit (Carbonate).

#### Devonian:

Two types of dolomites occur in the Devonian formations (Figure 3). These are: a pervasive replacive euhedral to subhedral micro to fine crystalline matrix dolomite (D1) and a pervasive replacive euhedral to subhedral medium crystalline, zoned matrix dolomite (D2) that is also associated sometimes associated with dissolution seams (Figure 8E). They are primarily observed in the Lucas and Amherstburg formations (Table 1; Figure 8).

**Figure 8.** (**A**) Photomicrograph (XPL) of planar-s fine crystalline dolomite matrix (D1) and planar-e to planar-s, replacive medium crystalline dolomite matrix (D2), showing increasing intercrystalline porosity. DGR8-CR22-59.06 (8-4) Amherstburg. (**B**) Photomicrograph (XPL) of chert (Sil) and silicified skeletal fragment partially replaced by fine crystalline dolomite (D1). Sample: DGR4-CR26-105.11 (4-5) Bois Blanc. (**C**) Photomicrograph (PPL) of planar-e to planar-s, medium crystalline dolomite matrix (RD2) pervasively replacing precursor limestone and partially replacing allochems and skeletal fragments. Sample DGR1-CR2-3-27.75 (1-4) Lucas. (**D**) Photomicrograph (PPL) of micro to fine crystalline dolomite matrix (D1) and planar-e to planar-s, replacive medium crystalline dolomite matrix (D2) associated with dissolution seams showing dark cores and clear rims. Sample: DGR3-CR6-42.30 (3-4) Lucas. (**E**) Photomicrograph (PPL—after staining) of planar-e, replacive medium crystalline dolomite matrix (D2) associated with dissolution seams Characterized by cloudy cores and clear rims. The undolomitized limestone matrix (LM) shows a pink color after staining. Sample: DGR1-CR3-2-31.85 (1-6) Amherstburg. (**F**) Photomicrograph (CL) of fracture-filling authigenic quartz (Qz) postdating dull to dark red luminescent medium crystalline dolomite matrix (D2). Sample: DGR4-CR26-105.11 (4-5) Bois Blanc.

#### 4.1.2. Calcite Cementation

Early and late calcite type of cement are present in the investigated Paleozoic formations. These types of cement display various crystal habits and sizes and occlude both primary as well as secondary porosity, pre- and post-dolomitization (Figures 4 and 9). The types of calcite cement identified include isopachous, syntaxial overgrowth, equant, dog-tooth, and blocky cement (Table 1).

Within Cambrian strata, fractures are occluded by a coarse, equant/ blocky calcite cement with strong, bright CL (Figure 9). In contrast, calcite cementation in the Ordovician formations is represented by an early equant calcite cement surrounding fossil allochems, isopachous rims around fossils and occluding hairline fractures (Figure 9), and later blocky calcite in fractures and vugs. The blocky calcite (BKC) in fractures postdates saddle dolomite cement (Figure 9).

In the Silurian successions, isopachous calcite cement growing around coated grains and exhibits microcrystalline, bladed texture (Figure 9; Table 1). Syntaxial calcite overgrowth cement (SXC; Figure 9D) forms ferroan, bright-luminescent crystals. Void-filling equant, pore-lining drusy calcite cement is observed in intergranular and intraskeletal pores, molds, and fractures. Late fracture and void-filling blocky calcite cement (BKC) postdate saddle dolomite cement. In terms of paragenesis, early calcite cement predated D1 and D2 and predated fracture- and pore-filling saddle dolomite which is postdated by late blocky calcite (BKC; Figures 4, 7F and 9) and evaporite cement.

**Figure 9.** (**A**) Photomicrograph (CL) of fracture-filling, zoned, blocky calcite cement (BKC) characterized by dull to nonluminescent cores and bright luminescent rims. Sample: T006056- 57-396.1 (2-15) Coboconk. (**B**) Photomicrograph (CL) of fracture-filling, bright luminescent blocky calcite cement (BKC) postdating zone, dull to non-luminescent coarse saddle dolomite cement (SD) and dark red luminescent medium crystalline dolomite matrix (D2). Sample: DGR2-CR133-1-843.96 Cambrian. (**C**) photomicrograph (XPL) of bladed vug-lining isopachous calcite cement (ISC) predating vug-filling blocky calcite cement (BKC). Sample: DGR6-CR243-898.24 Gull River. (**D**) Photomicrograph (XPL) of syntaxial overgrowth calcite cement (SXC) on crinoidal and echinoid skeletal fragments in a bioclastic grainstone. Sample: DGR1-CR154-440.92 (1-33) Manitoulin. (**E**) Photomicrograph (XPL) of drusy calcite cement (DC) and micro to fine crystalline dolomite matrix (D1) in fenestral dolostone. Sample: DGR1-CR49-142.28 (1-17) Bass Island. (**F**) Photomicrograph (XPL) of bladed isopachous calcite cement (ISC) filling interparticle pores between coated grains in a partially dolomitized grainstone. DGR8-CR122-332.26 (8- 10) Salina A1-Unit (Carbonate). (**G**) Photomicrograph (XPL) of bladed vug-lining isopachous calcite cement (ISC) predating vug-filling blocky calcite cement (BKC) and gypsum (Gy). Sample: DGR1-CR119-337.32 (1-27) Salina A1-Unit (Carbonate). (**H**) Photomicrograph (CL) of fracture-filling, zoned, blocky calcite cement (BKC) postdating dark red luminescent fine crystalline dolomite matrix (D1). Sample: DGR3-CR10-51.60 (3-5) Lucas. (**I**) Photomicrograph (PPL—after staining) of non-ferroan, drusy calcite cement (DC) lining a fine crystalline dolostone. The coarser crystals show zonation characterized by iron-rich cores followed by non-ferroan rims. Sample: DGR3-CR335.35 (3-2) Lucas.

> In the Devonian formations, four types of calcite cement occlude the pore spaces represented by syntaxial overgrowth, dog-tooth, equant/drusy, and blocky. Paragenetically, early calcite cement (e.g., syntaxial overgrowth, dog-tooth; Figures 4 and 9) predate replacive dolomite matrix D1 and predate the late fracture- and pore-filling blocky calcite (BKC) cement.

#### 4.1.3. Other Diagenetic Phases

Together with the calcite and dolomite diagenetic phases reported above, silica-bearing minerals and sulfate minerals are also present in the studied successions (Figure 10).

**Figure 10.** (**A**) Photomicrograph (XPL) of anhydrite cement (ANH) postdating coarse saddle dolomite cement (SD) and planar-s medium crystalline dolomite matrix (D2). Sample: DGR3-CR279-840.20 Cambrian. (**B**) Photomicrograph (XPL) of anhydrite cement (ANH) postdating coarse saddle dolomite cement (SD) and planar-e to planar-s medium crystalline dolomite matrix (D2). Sample: DGR3-CR279-2-839.80 Gull River. (**C**) Photomicrograph (XPL) of gypsum cement (Gy) postdating coarse saddle dolomite cement (SD) and micro to fine crystalline dolomite matrix (D1). Sample: DGR8-CR116- 313.31 (8-9) Salina A2-Unit (Carbonate).

Quartz is common in the Cambrian rocks typically replacing skeletal and nonskeletal components such as corals in coral-bearing facies and coated grains, such as ooids (Figure 5G). Silicification pre-dates dolomitization but postdates early isopachous calcite cementations (Figure 4). In the Ordovician strata, silicification replaced part of the muddy calcitic matrix, with microcrystalline dolomite (D1) replacing the silicified micrite and fossils. Anhydrite and gypsum are present in Cambrian, Ordovician, and Silurian rocks as fracture-fill, which postdates replacive dolomite (Figure 10). The rare occurrence of authigenic K-feldspar is observed by SEM filling some pores in Ordovician and Silurian rock but predated D2. Also, very minor pyrite and fluorite are observed in Silurian dolostones.

#### *4.2. Geochemical and Fluid Inclusion Results*

#### 4.2.1. Oxygen, Carbon and Sr Isotopes

Table 2 summarizes the stable and radiogenic isotope results of dolomite and calcite cement types. Figures 11–15 show the relationships among these isotopes with respect to their petrographic characteristics.

**Figure 11.** Cross plot of δ13CVPDB vs. δ18OVPDB of different generations of dolomite. Note the overlap in many values and the shift from the postulated range for respective seawater composition.


**Table 2.** Summary of the isotopic composition of dolomite and calcite components in the studied successions.


n. 5 5 4 Avg. −1.58 −9.12 0.70803 Stdev. 0.43 0.85 0.00003 Max. −1.22 −8.01 0.70806 Min. −2.40 −10.54 0.70798

**Table 2.** *Cont.*

Devonian BKC

**Figure 12.** Cross plot of δ13CVPDB vs. δ18OVPDB of different generations of calcite cement of different ages. Note the overlap in many values and the negative shift of late calcite from the postulated range for respective seawater composition.

**Figure 13.** Cross plot of δ18OVPDB values vs. 87Sr/86Sr ratios of different generations of dolomite. Note the overlap in many values and the shift from the postulated values for respective seawater composition except for dolomite associated with dissolution seams.

**Figure 14.** Cross plot of δ18OVPDB values vs. 87Sr/86Sr ratios of different generations of calcite cement. Note the highly radiogenic ratios in Cambrian blocky calcite.

**Figure 15.** Comparison between 87Sr/86Sr secular variation [69] and 87Sr/86Sr ratios of the investigated dolomite.

Fine-crystalline matrix dolomite in the Cambrian carbonates (D1) has a wide range of <sup>δ</sup>18O values from −16.9 to −2.9‰, <sup>δ</sup>13C values from −7.6 to −0.3‰ and 87Sr/86 Sr ratios from 0.70928 to 0.71716, respectively. Medium crystalline dolomite (D2) δ18O values vary from −8.2 to −6.9 ‰, <sup>δ</sup>13C values vary from −2.9 to −1.5‰, and 87Sr/86 Sr ratios from 0.70927 to 0.71100. Saddle dolomite cement <sup>δ</sup>18O values range between −10.7 to −8.7‰, from −4.2 to −0.5‰ for <sup>δ</sup>13C, respectively. In all dolomite types there is a co-variant

relationship between both isotopes (Figure 11). Blocky calcite cement (BKC) occluding fractures have <sup>δ</sup>18O values ranging from −15.2 to −7.9‰ for <sup>δ</sup>18O and <sup>δ</sup>13C values −5.3 to −3.5‰. Its 87Sr/86 Sr ratios range from 0.71000 to 0.71031 (Figure 12).

The fine crystalline replacive matrix dolomite (D1) in the Ordovician carbonates has <sup>δ</sup>18O values vary from −12.4 to −5.3‰, and less negative <sup>δ</sup>13C values compared to the Cambrian D1 dolomite (ranging from −1.2 to 1.1‰; Figure 11). Its 87Sr/86 Sr ratios range from 0.70808 to 0.72346, respectively. The highly enriched ratios are samples that occur along dissolution seams (Figure 13). D2 overlaps with D1 with δ18O values ranging from −9.6 to −4.1‰, <sup>δ</sup>13C values from −0.4 to 2.2‰, and 87Sr/86 Sr ratios from 0.70814 to 0.71255, respectively. The enriched ratio also represents a sample occurring along dissolution seams (Figure 13). The <sup>δ</sup>18O values for saddle dolomite (SD) vary from −10.8 to −8.3‰, and <sup>δ</sup>13C values from −2.5 to 2.2‰, respectively. Its 87Sr/86 Sr ratios vary from 0.70844 to 0.71066. The <sup>δ</sup>18O values for the early calcite matrix ranging from −9.2 to −3.8‰, and from −3.9 to 2.1‰ for carbon isotopes. One sample 87Sr/86 Sr ratios for matrix calcite is 0.70809. For late blocky calcite cement (BKC) the <sup>δ</sup>18O values vary from −9.8 to −3.9 and from −1.8 to 1.4‰ for <sup>δ</sup>13C. The 87Sr/86 Sr ratios for this cement vary from 0.70797 to 0.70809.

There is a significant overlap in δ18O and δ13C values in D1, D2, and SD. Finecrystalline matrix dolomite in the Silurian carbonates (D1) has δ18O values ranging from −8.0 to −4.4‰, <sup>δ</sup>13C values from −2.8 to 3.9‰ and 87Sr/86 Sr ratios from 0.70871 to 0.71201, respectively. Medium crystalline dolomite (D2) <sup>δ</sup>18O values vary from −8.3 to −3.8, <sup>δ</sup>13C values vary from −0.3 to 4.0‰, and 87Sr/86 Sr ratios from 0.708447 to 0.70891. Saddle dolomite cement <sup>δ</sup>18O values range between −9.6 to −6.1‰, from 0.1‰ to 3.69 for δ13C, respectively. Its 87Sr/86 Sr ratios vary from 0.70852 to 0.70872. The δ18O values for the early calcite matrix ranging from −8.8 to −4.7‰, and from −0.8 to 3.5‰ for carbon isotopes. In one sample, 87Sr/86 Sr ratios for the matrix calcite is 0.70815. Blocky calcite cement occluding fractures has <sup>δ</sup>18O values ranging from −10.1 to −4.5 for <sup>δ</sup>18O and <sup>δ</sup>13C values −3.75 to 3.1‰. Its 87Sr/86 Sr ratios range from 0.70805 to 0.70838.

Most of the δ18O values for D1 in the Devonian carbonates are slightly more enriched than D2 and they vary between −7.0 to −4.4 for <sup>δ</sup>18O and from 0.5 to 4.2 for <sup>δ</sup>13C, respectively. Its 87Sr/86 Sr ratios range from 0.70788 to 0.70937. Medium crystalline dolomite (D2) <sup>δ</sup>18O values vary from −10.3 to −3.9, <sup>δ</sup>13C values vary from 0.6 to 4.6‰, and 87Sr/86 Sr ratios from 0.70809 to 0.70937. The δ18O values for the early calcite matrix ranging from −9.6 to −4.4‰, and from 2.4 to 4.1‰ for carbon isotopes. The one analysis of Sr isotopes of this calcite yields a ratio of 0.70805. Blocky calcite cement (BKC) occluding fractures and vugs has <sup>δ</sup>18O values ranging from −10.5 to −8.0 for <sup>δ</sup>18O and <sup>δ</sup>13C values −2.4 to −1.2‰. Its 87Sr/86 Sr ratios range from 0.70798 to 0.70806.

#### 4.2.2. Major, Minor, Trace and REE Elements Content

Analytical results of major, trace, and rare-earth (REE) element concentrations for the dolomites and calcites from various age groups are summarized in Table 3. All analyzed replacement and cement dolomite types are non-stoichiometric, consistently calcium-rich (average CaCO3 ranging from 54.3 to 61.8 mol% for D1 and D2, and from 57.7 to 63.3 mol% for SD, respectively; Figure 16). There are no significant differences in Mg content in zoned D2 and SD between the cores and rims, as evidenced from the semiquantitative EDX analysis. In general, Mn and Fe contents in the Cambrian and Ordovician dolomites are much higher than the Silurian and Devonian dolomites (Table 3). There is also a covariant trend between Mn and Fe for these dolomites (Figure 17). However, the average of Sr content in all dolomites does not show large variations among the investigated dolomite types (vary between 83 and 253 ppm, Table 3), but many samples from the Ordovician dolomites show high Sr content and covariant trends with Fe and Mn (Figure 17). Such a relationship is not apparent in the dolomites from other age groups.


**Table 3.** Summary of the major and trace element composition of the studied successions.

**Figure 16.** Histogram plots of CaCO3 (mol%) of: (**A**) D1; (**B**) D2, and (**C**) SD.

**Figure 17.** Cross plots of the investigated dolomites: (**A**) between Fe and Mn; (**B**) Sr vs. Mn concentrations, and (**C**) Sr vs. Fe.

Samples of dolomite and fracture- and vug-filling calcite cement from different age groups were analyzed for REE content (Table 4). The ∑REE of Cambrian dolomite types and calcite cement is significantly higher than those in Ordovician, Silurian, and Devonian strata (Table 4) compared with the average of warm water brachiopods shown by Azmy et al. [70]. ∑REE of Devonian samples record the lowest values among the rest.

Replacive dolomite D1 and D2 sampled from Cambrian formations show different trends (Figure 18) compared to the pattern of warm water brachiopods shown by Azmy et al. [70], with a dominant-negative La anomaly and positive Ce anomalies (Figure 19).

**Figure 18.** Post-Archean Australian Shale (PAAS) normalized rare-earth element (REE) pattern for average values of fine (D1) and medium (D2) crystalline dolomite, and saddle dolomite (SD) samples, compared with average PAAS normalized REE pattern of water brachiopods [70].

**Table 4.** Summary of the REE (ppm) composition of dolomite and calcite components in the studied successions.



**Table 4.** *Cont.*

**Figure 19.** Cross plot Ce (Ce/Ce\*)SN–La (Pr/Pr\*)SN anomaly of D1, D2, and SD samples of different age groups.

Both dolomite D1 and D2 (Figure 18) sampled from Ordovician formations show different trends compared to the pattern of warm water brachiopods shown by Azmy et al. [70], with a negative La anomaly and both cases of positive and negative Ce anomalies (Figure 19). Both D1 and D2 exhibit higher average ΣREE (11.63 ± 21.56 ppm, and 8.54 ± 7.36 ppm, respectively) compared with the average of warm water brachiopods shown by Azmy et al. [70]. Saddle dolomite cement (SD) patterns show lower average ΣREE (9.55 ± 4.36 ppm) than D1 but slightly higher than D2. It also shows, similarly to D1 and D2, slightly different patterns compared to warm water brachiopod's pattern with both negative La and Ce anomalies.

Replacive dolomite D1 and D2 sampled from Silurian formations show different trends compared to the pattern of warm water brachiopods shown by Azmy et al. [70] (Figure 18), with a minor negative La anomaly and both cases of positive and negative Ce anomalies (Figure 19). Brachiopods, D1, and D2 patterns with both negative La and Ce anomalies as well as a slight negative Eu anomaly (Figure 19).

Both fine and medium dolomite matrix samples (D1 and D2, respectively) from Devonian formations have shale-normalized patterns sub-parallel to those of modern warm water brachiopods from Azmy et al. [70]. They both show the typical characteristics of seawater patterns such as depletion of LREE over HREE and slightly positive La and Gd anomalies (Figure 19).

There is a subtle difference in Cambrian and Ordovician dolomite vs. Silurian and Devonian samples in terms of Sm/Yb vs. Eu/Sm (Figure 20A), La/Yb vs. La/Sm (Figure 20B), Y/Ho vs. Sm/Yb (Figure 20C) and Y/Ho vs. Eu/Sm (Figure 20D) ratios. For example, higher Sm/Yb ratios (Figure 20A), lower Sm/Yb, Eu/Sm (Figure 20C,D) in the Cambrian and Ordovician dolomites than Silurian and Cambrian dolomites.

**Figure 20.** (**A**) Plots of Sm/Yb vs. Eu/Sm; (**B**) La/Yb vs. La/SmN; (**C**) Y/Ho vs. Sm/Yb, and (**D**) Y/Ho vs. Eu/Sm of dolomite samples from different age groups.

#### 4.2.3. Fluid Inclusions Microthermometry

Microthermometric measurements of primary fluid inclusions from D2, SD, and BKC from each age group provide melting and homogenization temperatures (Th) and estimates of salinity (cf. [61,71] for each phase (Table 5 and Figure 21). Primary fluid inclusions within D1 crystals were too small to measure (cf. [71]). However, monophase fluid inclusions are present in this dolomite. For calcite, measurements were taken for coarse crystalline cements (BKC) in vugs and fractures. The measured Th provides an estimate of minimum entrapment temperatures. The low first melting temperatures and ice-melting temperatures suggest that the fluid within the fluid inclusions may approximate a H2O-NaCl-CaCl2 system. The examined inclusions are parallel to crystal facets and the vapor bubble size for the two-phase (liquid-vapor) fluid inclusions. There is a notable overlap in Th and salinity values among all dolomite types in Cambrian, Ordovician and Silurian samples (Figure 21), showing high Th and salinity. However, lower Th and salinity values are recorded for D2 in the Devonian samples. Fluid inclusions within D2 from Cambrian crystals yielded Th that range from 75◦ to 132 ◦C (102.9 ◦C, n = 19), with salinity estimates based on final ice-melting temperatures (Tmice) that correspond to salinities range from 23.2 to 27.2 wt.% (25.4 wt.%, n = 10). Fluid inclusions in SD crystals have Th ranging from 111◦ to 156 ◦C (131.6 ◦C, n = 5). Measurable 2-phase fluid inclusions in BKC giving Th ranging between 85◦ and 141 ◦C (108 ◦C, n = 11), and salinities from 22.1 to 23.6 wt.% (22.9 wt.%, n = 2).


**Table 5.** Summary of the fluid inclusion data of the studied successions.


**Table 5.** *Cont.*

Fluid inclusions within D2 from Ordovician samples yielded Th that range from 66◦ to 132 ◦C (96.3 ◦C, n = 41), with salinity, estimates that range from 21.6 to 30.5 wt.% (26.5 wt.%, n = 12; Table 5). Fluid inclusions in BKC giving Th ranging between 68◦ and 153 ◦C (109.7 ◦C, n = 27), and salinities from 24.8 to 30.3 wt.% (27.2 wt.%, n = 7).

In Silurian sample Th values for D2 are ranging between 49.7 ◦ and 134.1 ◦C (93.4 ◦C, n = 42), and salinities from 15.1 and 25.2 (23.1, n = 15). Fluid inclusions in SD crystals have Th ranging from 101.2◦ to 193.4 ◦C (123.5 ◦C, n = 53), and salinities from 25.2 to 32.6 (28.5 wt%, n = 29). Measurable 2-phase fluid inclusions in BKC giving Th ranging between 70.2◦ and 194.3 ◦C (129.1 ◦C, n = 53), and salinities from 21.4 to 32.3 wt.% (27.6 wt.%, n = 25).

Fluid inclusions within D2 from Devonian samples yielded Th that range from 69.9◦ to 102.3◦C (83.7 ◦C, n = 29), with salinity range from 18.9 to 21.8 wt.% (20.7 wt.%, n = 12). For BKC higher Th values than samples from other age groups are measured and vary between 146.2◦ and 255.1 (202.0 ◦C, n = 48), and lower salinity values vary between 15.1 and 23.9 wt.% (20.7 wt.%, n = 25).

**Figure 21.** Cross plots of (**A**) Homogenization temperature vs. salinity of dolomite types; (**B**) calcite phases; (**C**) Histogram plot of Th for dolomite samples; (**D**) Histogram plot of Th for calcite samples; (**E**) Histogram plot of salinity of late calcite; and (**F**) Histogram plot of salinity for dolomite types.

#### **5. Data Analysis and Discussion**

Petrographic investigations, stable-isotope, elemental analyses, and fluid-inclusion microthermometry provide important clues regarding dolomitization and other diagenetic processes, fluid chemistry, and the history of fluids migration in the sedimentary basin. The Paleozoic successions in the Huron Domain of Michigan Basin were affected by

early (marine to shallow eogenetic zone; [72]) and deeper (mesogenetic zone) diagenetic processes [72]. These include early calcite cementation, silicification, partial dissolution, early dolomitization and later dolomitization, recrystallization, fracturing, stylolitization, late calcite, and sulphate cementation. There are similarities in the diagenetic processes in Cambrian and Ordovician carbonates and Silurian and Devonian carbonates which allow us to consider and discuss these carbonates together (Figure 4).

#### *5.1. Constraints from Petrography-Paragenetic Sequence*

The paragenetic sequence of the main diagenetic events in the studied Paleozoic successions (Figure 4), which was constructed originally based on textural relationship and geochemical and fluid inclusion evidence, provides a useful platform for further discussion of the data obtained in this study.

As reviewed in Section 2 above, the depositional environments of the Cambrian and Ordovician carbonates encompass various lithofacies assemblages, including supratidal, lagoonal to subtidal environments (e.g., [47,55]). The diverse lithofacies changes observed in these formations are directly controlled by relative sea-level changes and sources of detritus. These were evident especially in the Cambrian strata where mixed siliciclasticcarbonate facies are common. Ultimately, these depositional processes also affected the diagenetic changes encountered in these sediments during shallow and deeper burial stages. In carbonate-rich units, the bioclastic grainstone facies has been influenced by early stabilization of calcitic/aragonitic components and dolomite replacement at shallow burial (D1; Figure 5B) to be succeeded by medium to coarse crystalline dolomite (D2; Figure 5B–D,I), possibly at an intermediate burial. Minor presence of saddle dolomite (SD) occluding vugs postdates D2 but predates BKC, which is occluding fractures; both of which probably formed later than D2 in the intermediate burial stage (Figure 4; Figure 5E). In the ooid grainstone lithofacies, which is extensively silicified (Figure 5G), the early isopachous calcite cement was silicified. However, the silicified ooids are well preserved and show evidence of very early diagenetic, the pre-compaction origin of these allochems. The siliciclastic-dominated facies were cemented earlier by calcite cement later to be compacted with the creation of some minor fractures occluded by equant calcite cement, which was partly dolomitized. The basal Cambrian sandstone is assumed to have been a major migratory path for diagenetic fluids [73,74]. Strata in the eastern region of the basin pinches out towards the margins. These thin fingers of clastic material which have large amounts of bioclastic material increase the likelihood for extensive dolomitization of a thin layer.

In the Ordovician formations (e.g., Coboconk Formation) diagenesis was initiated by early cementation that produced the nodular facies on the seafloor. Burial and compaction affected the formation of these nodules and may have created hairline fractures, later occluded by dolomite. Dolomitization (D1; Figure 6) occurred by selective early replacement of skeletal components such as bryozoans and echinoderms and the micritic matrix, as well as, pervasively replacing the intra-nodule matrix during shallow burial. These dolomites were crosscut by hairline fractures and stylolites during a burial at a later stage, suggesting the formation of D1 at a shallow burial. Silicification, which occurs in some facies of the Coboconk Formation occurred prior to dolomitization as evidenced by silica replacing the matrix and fossils that later were selectively dolomitized. In the Gull River Formation, a similar paragenetic sequence to the Coboconk Formation is observed (Figure 4). Initially, dolomitization proceeded with fine crystalline (D1) and medium to coarse crystalline matrix replacement dolomite (D2). The minor saddle dolomite (SD) cement postdates earlier dolomite formation and occludes vugs and fractures and is also postdated by late, blocky calcite cement (BKC), which is occluding vugs and fractures (Figure 6). Other compaction features such as stylolites postdate or sometimes appear to be contemporaneous with the medium to coarse crystalline dolomite (D2; Figure 6C). In both the Coboconk and Gull River formations, the intra and internodule dolomites typically present as petrographically distinct. These distinctions reflect differences in the underlying substrate that was originally dolomitized and relate to permeability, not fluid

chemistry, or timing (notwithstanding some evolution of brine chemistry over time). An early, shallow burial origin is suggested for the fine crystalline (D1) dolomite as evidenced by the selective replacement of this dolomite, its shape and size, and it predates stylolite formation, consistent with a shallow burial diagenetic environment. In contrast, the medium and coarse crystalline (D2) dolomite may have formed at medium burial stages as evidenced by its increasing crystal size, subhedral and anhedral shape (e.g., [26]), and close association sometimes with compaction features, such as stylolites. The stylolites observed in the samples all occur at ~90◦ to the fractures. This relationship may be explained by compressive stress, imposed by the Taconic Orogeny, which was occurring concurrently with deposition. These stylolites may have acted as conduits for Mg-rich fluids (e.g., [75]).

In contrast, the depositional environments of the Silurian and Devonian successions encompass various lithofacies assemblages, including shallow marine ramp and deeper water subtidal lithofacies dominated by limestone and dolostone lithofacies with intercalations of shales and evaporite deposits in both age groups (e.g., [47,76,77]). In the Silurian and Devonian carbonates early calcite cementation (e.g., ISC and SXC; Table 1) and micritization initiated on the seafloor to be followed during burial by occlusion or porosity by later calcite cementing vugs and fractures (e.g., DC, and BKC; [78]). BKC postdates saddle dolomite (SD; Figure 4). Evaporite minerals, such as gypsum and anhydrite represent the final stage of cementation in the Silurian rocks. Dolomitization commenced early (D1) formed replacing calcitic and aragonitic allochems and muddy calcitic matrix [3]. This dolomite may have formed at a shallow burial, in a reflux evaporative setting (e.g., [33,79,80]). This is supported by its fine crystalline crystal size [68] and the cross-cutting relationship with compactional features, such as dissolution seams and stylolites which crosscut this dolomite. This dolomitization phase is followed by the formation of medium to coarse crystalline replacive dolomite (D2) at the medium burial stage. The increase in crystal size of D2 may reflect dolomitization of coarser lithofacies (such as grainstones) and/or recrystallization of precursor dolomite matrix (D1; [23]). In the Silurian carbonates, only saddle dolomite (SD) is observed occluding fractures, postdating D1 and D2 but predating blocky calcite (BKC; Figure 7). This dolomite formed possibly in an intermediate burial setting.

Evaporite minerals (gypsum and anhydrite; Figure 10) represent the last diagenetic phases and form in a deeper burial realm. These minerals occlude fractures and vugs in Cambrian, Ordovician, and Silurian, carbonates.

#### *5.2. Controls of Dolomitization-Geochemical Evidence*

#### 5.2.1. Stable and Radiogenic Isotopes

Cambrian:

The stable isotopic composition of D1, D2, and SD in Cambrian rocks (Figure 11) are notably lighter than the postulated dolomite precipitated in equilibrium with marine waters (e.g., [69]). These values are also distinct from other age groups (Figure 11). There is little isotopic fractionation of 13C/12C with temperature, and hence the δ13C values of these dolomite reflect isotopic fractionation of C related to fluids that are enriched in lighter CO2 originated from oxidized carbon, which can be related to seawater-sourced, reduced, bacterial sulfate reduction process for D1 and hydrocarbon sources for D2 and SD [2].

However, the δ18O values of the dolomites reflect both the temperature and the composition of the parent fluids, and thus the departure from equilibrium values reflects precipitation of these dolomite under a different diagenetic regime [81]. The negative shift in δ18O of D1 and D2 (Figure 11) signifies the effect of temperature increase during burial, change in fluid composition, or both [81]. Recrystallization of D1 with increasing temperature/change in fluid chemistry can also result in a negative shift in D1 (e.g., [23]). This recrystallization is also observed petrographically, as shown by the increasing crystal size of D1 and the presence of overgrowth zonation in D2 (Figure 5). The co-variant trend between carbon and oxygen isotopes in these dolomites reflect varied water-rock interaction with increasing burial (cf. [82]). Water-rock interaction is also evidenced by the presence of more radiogenic 87Sr/86Sr ratios in D1 and D2 (Figures 13 and 15) in comparison

to the postulated range for Cambrian seawater. This can be related to the presence of more siliciclastics, rich in clays and feldspars, in the Cambrian rocks as well as can be sourced from the Precambrian shield [74,83]. The same trend and comparable negative isotopic values are evident in late calcite cement (BKC; Figure 12) and thus reflecting similar diagenetic conditions during later fluid flow events for these carbonates unique to Cambrian paleohydrologic flow system. Similarly, BKC also shows more radiogenic Sr isotope ratios than the Cambrian seawater values (Figure 14).

Ordovician:

Most of the oxygen and carbon isotopic values of the Ordovician replacive dolomites (D1 and D2) show overlapping signatures with a departure by at least 3‰ from their coeval postulated seawater equilibrium values with respect to their age (Figure 11), especially for oxygen isotopes. This can be related to the recrystallization of D1 to D2 at elevated temperatures during the intermediate burial stage (Figure 4). C isotope values are closer to their seawater equilibrium values reflecting the buffering of this isotope with original carbonate values (e.g., [26,29]). However, δ13C values of saddle dolomite show slight negative shifts compared to the postulated equilibrium values. This may reflect an oxidized organic C source during burial. The 87Sr/86Sr ratios of many samples of D1 and D2 show comparable and/or slight radiogenic signature to seawater values (Figure 13; [69]) with the exception of few dolomites that show highly radiogenic rations (Figures 13 and 14). These dolomite samples occur along dissolution seams and thus reflecting a radiogenic fluid source of Sr, perhaps from Cambrian siliciclastic fluid. In contrast, the isotopic composition of early calcite cement shows preservation of the original seawater signature (Figure 12), while late calcite cement shows some negative trend in both oxygen and carbon isotopes relative to the postulated values in Ordovician calcite. This reflects the diagenetic progression and precipitation of this cement from fluids modified at the burial. This is also shown by 87Sr/86Sr ratios of BKC which show seawater signature (Figure 14).

Silurian and Devonian:

The isotopic data of the overlying Silurian and Devonian carbonates show overlaps between δ13C and δ18O values (Figure 11). Dolomite samples from both age groups show similar δ13C values compared to those of carbonates deposited in equilibrium with seawater of respective age except for few D1 samples from the Silurian carbonates that show a slight negative shift, which can be related to a depleted carbon source related to bacterial sulfate reduction at the shallow burial [84]. This is supported by the presence of minute framboidal pyrite crystals embedded within RD1 (this study and Tortola et al. [3]). Also, few BKC samples from the Silurian carbonates show a negative shift from the equilibrium values (Figure 12). This can be related to an oxidized organic carbon source, possibly related to interaction with hydrocarbons (e.g., [2]). However, δ18O values of D1 show evidence of dolomite recrystallization as these values depart from their respective age equilibrium values. This negative shift in oxygen isotopes also reflects increasing temperature and/or changes in the diagenetic fluid composition with increasing burial (cf., [3,26]). 87Sr/86Sr ratios of many samples, both dolomite and calcite, from these age groups cluster close to their coeval seawater composition (Figure 13). However, a few dolomite samples show slight enrichment of radiogenic Sr ratios. This can be related to the presence of modified seawater fluids upon burial.

#### 5.2.2. Major, Minor, and REE Elements

The average values of CaCO3 in all analyzed replacement and cement dolomite types show that all dolomite types are non-stoichiometric, consistently calcium-rich (Table 4, average CaCO3 range from 54.3 to 61.8 mol% for D1 and D2, and from 57.7 to 63.3 mol% for SD, respectively; Figure 16). This suggests that the precursor mineralogy of these dolomites is of high-Mg calcites and aragonite and the stabilization of these minerals to dolomite occur in semi-closed diagenetic systems with less fluid mobility and limited supply of Mg ions (cf., [18,85]). Hence, dolomitizing fluids reflect a lowering Mg/Ca ratio in these fluids which encourages Ca-rich, non-stoichiometric dolomite to form. In general, Mn and Fe contents in the Cambrian and Ordovician dolomites, which are indicators of redox conditions during dolomitization, are much higher than the Silurian and Devonian dolomites (Table 3) reflecting a decrease in oxidizing conditions in the older successions. There is also a covariant trend between Mn and Fe for these dolomites (Figure 17) reflecting the geochemical affinity of these two elements [86]. However, the average of Sr content in all dolomites does not show large variations among the investigated dolomite types (vary between 83 and 253 ppm, Table 3) and are comparable to those of other ancient dolomites but lower than the estimated values for dolomite in equilibrium with seawater (~500–800 ppm, (e.g., [87]). However, many samples from the Ordovician dolomites show high Sr content and covariant trends with Fe and Mn (Figure 17), perhaps reflecting the aragonitic precursor mineralogy of these dolomites. Such a relationship is not apparent in dolomites from other age groups.

Rare earth elements (REEs) have been commonly utilized for the reconstruction of a paleoenvironmental history of ancient seawater [88–90], diagenetic alteration, and the temperature of fluids precipitating carbonate minerals [70]. REE can also provide information about the relative water depth and the proximity to hydrothermal sources [91]. Some workers indicate that carbonates might preserve their REE compositions and/or patterns during diagenetic processes [89,92]. However, many studies have revealed that consequent diagenetic alteration modifies REE concentrations and/or patterns [70,93]. Therefore, the redistribution of REE in carbonates is an indicator that can be utilized to ascertain the potential contributing factor to observed REE concentrations in diagenetic minerals, especially in high water/rock ratio diagenetic systems [70]. The rare earth elements generally exist in the stable trivalent (3+) oxidation state, whereas Eu (Eu2+, Eu3+) and Ce (Ce3+, Ce4+) exhibit distinct behavior due to the redox state of the natural system (e.g., [93]). Ce (Ce3+) can be oxidized to Ce (Ce4+) in marine systems under surface conditions, whereas high temperature and/or high reducing conditions cause Eu (III) to be reduced to Eu (II) [65,94]. These characteristics provide the utilization of REEs as tracers of chemical processes during the diagenetic history of carbonates [89,95–97].

The ∑REE of Cambrian dolomite types and calcite cement (Table 3) is significantly higher than those in Ordovician, Silurian, and Devonian strata when compared with the average of warm water brachiopods shown by Azmy et al. [70]. ∑REE of Devonian samples record the lowest values among those examined in this study. The higher concentration of ∑REE of Cambrian dolomites signify that the dolomitizing fluids contained higher concentrations of REE, likely related to the interaction of these fluids with basement rocks.

Replacive dolomite D1 and D2 sampled from Cambrian formations show different trends (Figure 18) compared to the pattern of warm water brachiopods shown by Azmy et al. [70], with a dominant-negative La anomaly and positive Ce anomalies (Figure 19). These dolomites may have formed under suboxic to anoxic conditions [70]) from early crustal fluids highly enriched in REE particularly MREE and HREE relative to dolomites from other age groups. Saddle dolomite from the Cambrian rocks show evidence of different fluids that produced D1 and D2 are likely related to hydrothermal fluids that originated from deeper sources and transported via fractures and faults [2]. The flat Ce\* anomaly, slightly positive Eu\*, flat LREE over HREE, and sight MREE enrichment with low La/Sm in saddle dolomite reflect higher temperature fluids of seawater origin under low oxic conditions modified by siliciclastic contamination.

Both D1 and D2 (Figure 18) sampled from Ordovician formations show different trends compared to the pattern of warm water brachiopods shown by Azmy et al. [70], with a flat to slight Eu anomaly, negative La anomaly, and both cases of positive and negative Ce anomalies (Figure 19). Both, D1 and D2 exhibit higher average ΣREE (11.59 ± 22.44 ppm, and 8.54 ± 7.36 ppm, respectively) compared with the average of warm water brachiopods shown by Azmy et al. [70]. Saddle dolomite cement (SD) patterns show lower average ΣREE (9.55 ± 4.36 ppm) than D1 but slightly higher than D2. It also shows, similarly to D1 and D2, slightly different patterns compared to Warm water Brachiopod's pattern with

both negative La and Ce anomalies. High Eu\*, high La/YbN\* and Y/Ho can indicate interaction with hydrothermal fluids [90].

Silurian replacive dolomites (D1 and D2) show different trends compared to the pattern of warm water brachiopods shown by Azmy et al. ([70]; Figure 18), with a minor negative La anomaly and both cases of positive and negative Ce anomalies (Figure 19). Brachiopods, D1 and D2 show patterns with both negative La and Ce anomalies as well as a slight negative Eu anomaly (Figure 18). These data suggest a seawater-dominated source of diagenetic fluids evolved during burial and increasing water/rock interaction [3]. The higher ΣREE and a more pronounced enrichment in LREE in SD and BK reflect a different source of REEs and/or involvement of saline fluids.

Both, fine and medium dolomite matrix samples (D1 and D2, respectively) from Devonian formations have shale-normalized patterns sub-parallel to those of modern warm water brachiopods from Azmy et al. [70]. They both show the typical characteristics of seawater patterns such as slight depletion of LREE over HREE and slightly positive La and Gd anomalies (Figure 19) potentially modified from the interaction of meteorically derived fluids and fluids of marine parentage. This can occur at near-surface, slightly reducing conditions. Although there is no evidence of karstification in the samples collected and analyzed in this study, a previous study [98] suggested that dissolution processes were mostly active in the shallow subsurface, usually <200 m depth in south-western Ontario. Thus, the possibility of meteoric water influence cannot be discounted completely.

There is a subtle difference in Cambrian and Ordovician dolomite vs. Silurian and Devonian samples in terms of Sm/Yb vs. Eu/Sm (Figure 20A), La/Yb vs. La/Sm (Figure 20B), Y/Ho vs. Sm/Yb (Figure 20C) and Y/Ho vs. Eu/Sm (Figure 20D) ratios. These differences may reflect changes in fluid composition, temperature, and/or redox conditions [99,100]. For example, Y/Ho vs. Sm/Yb ratios showing two possible fluids (marine and hydrothermal), and their interactions more evident in Cambrian and Ordovician dolomites. In contrast, the effect of seawater on the formation of Silurian and Devonian dolomites is more important. Other ratios document the presence of modified seawater as the primary fluid type in the formation of these dolomites [100].

#### 5.2.3. Fluid Inclusions

Previous studies on fluid inclusions in Cambrian to Devonian dolomite and calcite in both central and marginal parts of the Michigan basin demonstrated fairly a wide range of Th and salinity values (e.g., [2,3,32,33,101–103]). These investigations reported high Th and salinity values in the central part of the basin compared to the margin of the basin, which some of them attributed to hydrothermal fluids which originated from the central part of the basin. Haeri-Ardakani et al. [33] further suggest a thermal anomaly related to the mid-continent rift under the Michigan Basin during late Devonian-Mississippian is the major source of heat for hydrothermal fluids.

The presence of single-phase fluid inclusions in D1 suggests very low temperatures, up to 50 ◦C [71]. However, the salinity-homogenization temperature plots (Figure 21A,B) for both dolomite (D2 and SD) and calcite (BKC) demonstrated a clear and subtle difference between age groups. For instance, higher Th and salinity values are recorded in SD and some D2 in Ordovician and Silurian samples compared to Devonian dolomite. In Cambrian, D2 overlaps with Ordovician samples. In contrast, BKC in Devonian samples show the highest Th and lowest salinity values compared to the rest (Figure 21B). High salinity and homogenization temperatures values are shown by Silurian dolomite and calcite. This reflects the effect of fluids sourced from saline brines related to evaporite dissolution in the Silurian rocks (e.g., [3,32]). The differences in homogenization temperature and salinity values in BKC (Figure 21) between Silurian and Devonian may reflect divergent fluid systems, such as a fluid that is characterized by high temperature and high salinity in the Silurian rocks and much higher temperature and lower salinity fluids in the Devonian rocks. These phenomena can also reflect the effect of Alleghenian Orogeny in Carboniferous time (Figure 1) as a factor that affected the fluid composition and sources.

In all analyzed D2 and SD dolomite and late BKC in the investigated age groups the high Th values are significantly higher than the maximum burial temperatures suggested by other studies (e.g., [2,32,33,101,102]). This may suggest an involvement of hydrothermal fluids during the formation/recrystallization of dolomite and precipitation of late calcite cement. Increasing temperature of formation/recrystallization, as well as water/rock interaction, can explain the negative departures from equilibrium δ18O values of dolomite.

#### *5.3. Dolomite Recrystallization*

Dolomite crystal shape, its boundaries, and its characteristics are directly related to the crystallization environment in which they formed [104]. Recrystallization of earlier formed dolomite (e.g., [23,105–107]) can involve change in texture, such as coarsening of crystal size or increase in nonplanar boundaries, change in structure (progressive ordering and strain), change in composition, which include isotopes, trace elements, stoichiometry, fluid inclusions, and zonation. It can also include a change in paleomagnetic properties [11].

Recrystallization of D1and D2 dolomites in the studied rocks has been invoked because of petrographic and geochemical data. Alteration of early formed dolomite is a common diagenetic process during progressive burial and/or increasing water/rock interactions (e.g., [8,108,109]). SEM and petrographic evidence for recrystallization of D1 include: (1) increase of crystal size (Figure 22); (2) lack of zonation (e.g., overgrowth) under CL, especially in D2; (3) lack of stoichiometry in D1 and D2 can be related to recrystallization of D1 to D2 by Ca-rich, warmer basinal fluids, perhaps under a semi-isolated diagenetic system (e.g., [3,13,110]). Geochemical and isotopic evidence for recrystallization with increasing burial include: (1) a negative shift of δ18O from postulated values of marine dolomite of respective ages (Figure 11), (2) a slight enrichment of Sr isotopic ratios (Figure 13), (3) high Th values in D2 (Figure 21), and (4) trace elements fractionation in D1 and D2 (Figure 17).

**Figure 22.** SEM and BSE images of D1 (**A**), D2 (**B**,**C**) and SD (**D**).

#### *5.4. Evolution and Origin of the Diagenetic Fluids*

Fluid flow in sedimentary basins can be linked to the tectonic evolution of the basin which can involve tectonic thrusting, sedimentary loading, uplift, and compression [4]. These fluids can transfer dolomitizing fluids as well as metals and hydrocarbons. Recent studies on diagenetic fluid flow in southern Ontario provided models on the fluid composition and its source(s). Coniglio and Williams-Jones [111], Middleton et al. [55] and Coniglio et al. [112] advocated that lateral and vertical movement of Mg-bearing fluids in Michigan Basin resulting from the transformation of clay minerals, originating from younger, older, or correlative argillaceous sediments from deeper in the Michigan and Appalachian basins towards the basin margin was considered as a mechanism for migration of dolomitizing fluids. Yoo et al. [113] also considered compaction driven flow to explain the mass and heat from the Michigan Basin. Haeri-Ardakani et al. [32,33] suggested that hot brines in the central part of the basin migrated through the basal Cambrian sandstone and ascended through a fracture network precipitating dolomite and late-stage calcite while carrying hydrocarbons.

Hobbs et al. [114] suggested two geochemical systems present on a regional scale; a shallow system (<200 m) below ground surface containing fresh through brackish waters and an intermediate to deep (>200 m) system containing brines associated with hydrocarbon reservoirs. Petts et al. [115] using microthermometric and isotopic data of secondary minerals in Cambrian and Ordovician strata suggested the presence of hydrothermal brines mixed with connate water. Al-Aasm and Crowe [2] in their investigations of dolomitization in the Cambrian and Ordovician carbonates in the Huron Domain identified two somewhat isolated digenetic fluid systems; an earlier Cambrian system and a later Ordovician system, both displaying unique geochemical attributes. Bouchard et al. [116] suggested that the fluids in the Ordovician low permeability sequences in the eastern part of Michigan Basin have originated as Ordovician seawater that was chemically overprinted by evaporated seawater that infiltered during the late Silurian and mixed with pre-existing brines from a deeper part of the hydrostratigraphic section. More recently, Tortola et al. [3] suggested the presence of a diagenetic fluid system that affected Silurian carbonates and was altered by salt dissolution post-Silurian and a Devonian fluid system that affected by the Alleghanian orogeny and resulted in a different fluid source and composition.

Petrographic relationships and isotope data can be used to infer changes in the composition of pore fluids in the Paleozoic successions during progressive diagenesis. Using the mean Th values and the range of δ18O for D2 and SD from the investigated time frames, the δ18O composition of the precipitating fluids was determined. The relationships among mineral δ18O values, pore water δ18O values, and temperature for diagenetic phases from the investigated carbonates are shown in Figure 23A,B. The suggested pathways for pore water evolution reflect changes in temperature, changes in fluid chemistry, or both. For dolomite types/generations (Figure 23A), there is a subtle difference in temperature and hence fluid composition between D2 and SD in Cambrian strata. The dolomitizing fluids for D2 reflect marine parentage, while for SD it is more enriched in δ18O values suggesting the modified brines formed under higher temperatures. In contrast, there is a significant overlap in the δ18O composition of dolomitizing fluids for Ordovician D2, SD, and Devonian D2 and a comparable range of formation/recrystallization temperatures. However, there is a measurable difference between D2 and SD in Silurian samples, as D2 formed under lower temperatures and comparable seawater δ18O values but higher temperatures and more enriched values for SD formation. This can be related to an increase in salinity with burial, which can be linked to the Salina salt (dissolution) in the Upper Silurian, and the incursion of hydrothermal fluids during the formation of SD (e.g., [117]).

**Figure 23.** (**A**) Oxygen isotope values of dolomite the investigated formations plotted on the temperature-dependent, dolomite-water oxygen fractionation curve [118]. (**B**) Oxygen isotope values of calcite cement from the investigated formations plotted on the temperature-dependent, calcite-water oxygen fractionation curve [119].

The formation of blocky calcite (BKC) later in the diagenetic sequence reflects variations in formation temperatures and fluid composition in the investigated age groups (Figure 23B). All calcites formed at higher temperatures and deviate from the seawater composition of corresponding ages. However, similarities in δ18O fluid composition between BKC and D2 in the Cambrian samples suggest the fluids could have originated from the Precambrian rocks below and migrated to the overlying porous Cambrian strata (cf. [74]). In contrast, the composition and temperatures of fluids that formed the late calcite in other age groups changed, becoming more enriched and thus reflecting evolved basinal brines [2]. The pathways of these fluids are shown by increasing both δ18O values and temperature with time (from older to younger systems). Variations in pore fluid isotopic composition in the Paleozoic age groups suggest compartmentalization of these diagenetic fluids resulting in the formation and/or recrystallization of dolomite and calcite cementation. This phenomenon is also supported by radiogenic Sr isotopic ratios of these diagenetic phases as well as REE distribution of these carbonates (Figures 14 and 23).

Fluid compartmentalization in carbonate strata in the Huron Domain is also illustrated in Figure 24, which shows some of the important geochemical parameters used in this study vs. their stratigraphic position. There is a subtle difference in stable isotopic composition between different age groups, especially in the Cambrian samples. In all dolomite types there is a significant departure in δ18O values from the postulated range of dolomite that was precipitated in equilibrium with seawater of respective ages. This can be related to increasing temperature during burial as well as modification of pore fluid isotopic composition, in SD this negative shift can be related to the effects of hydrothermal fluid interaction with host rocks; perhaps limited in extent and distribution. Salinity and Th measurements in the Silurian also show a distinctive pattern with higher and salinity temperatures of formation comparing to other rock units, which suggest the effect of evaporite dissolution and migration of fluids in these rocks.

**Figure 24.** Stratigraphy vs. isotopic and microthermometric data for the investigate age groups.

The conceptual model for fluid flow mechanism in the Huron Domain shown in Figure 25 illustrates the concept of fluid compartmentalization in different stratigraphic rock units suggested in this study. The current burial depth is shown also in the model based on various studies [80,111,120]. As mentioned earlier this region of the Michigan Basin has experienced multiple stages of tectonic evolution, both passive and active tectonic events which affected fluid composition and migration pathways (e.g., [2,5,32,37]). Migration of warm, basinal brines from the center of Michigan Basin towards the margin [32] transported dolomitizing fluids and other solutes. The source of heat of these fluids could be related to the reactivation of the buried mid-continental rifts during late Devonian to Mississippian time [121]. A limited influx of hydrothermal fluids moved upwards from the Precambrian basement via fractures and faults [2,47] results in minor precipitation of saddle dolomite in fractures and vugs and altered earlier formation of matrix dolomite. This likely occurred during Taconian and Acadian Orogenies. Interaction of connate fluids of marine parentage with warm, basinal brines controlled the formation of dolomite and calcite cement at a shallow to intermediate burial depth. We believe that diagenesis is horizontally (stratabound) uniform across the Huron domain, displaying few signs of significant vertical connectivity beyond formation tops, except in areas with local heterogeneity (faults; [122]).

**Figure 25.** Conceptual model for fluid flow in the study area illustrating fluid compartmentalization in different stratigraphic rock units and migration pathways. These fluid pathways vary with time and effect of tectonic setting (stratigraphic setting is modified from [123]).

> Downward migration of saline fluids that originated from the dissolution of Silurian evaporites and possible interaction with basinal brines moving up-dip for the center of the basin was responsible for the modification of diagenetic phases [33]. In the Devonian hydrogeochemical system, the interaction of basinal brines with possible infiltration of younger fluids, perhaps of meteoric water origin, resulted in the formation of late dolomite and calcite cement (e.g., [55,111,112]). This probably occurred during the Alleghenian Orogeny (e.g., [3]). Sutcliff et al. [124] in their study of late calcite cement recovered from sub-horizontal fractures near the base of the Silurian sequence dated show U–Pb ages of 318 ± 10 Ma determined by LA-ICP-MS and 313 ± 1 Ma by ID-TIMS. Hence, these authors also suggested that the fluid flow of hydrothermal brines was influenced by the Alleghanian Orogeny which caused the migration of hydrothermal fluids westwards from high-level Alleghanian plutons through the basal Cambrian rocks.

#### **6. Conclusions**

Based on petrographic, geochemical, and fluid inclusion data presented, the following conclusions have been drawn with regards to the diagenetic evolution of the Huron Domain of southern Ontario:


**Author Contributions:** Conceptualization, I.S.A.-A.; writing—original draft preparation, I.S.A.-A.; supervision, I.S.A.-A.; writing, review and editing, R.C.; writing, review and editing, M.T. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research is funded by Nuclear Waste Management Organization (NWMO) contract #816150 to ISA.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** The data will be available with the authors.

**Acknowledgments:** I.S.A.-A. and M.T. acknowledge the help from M. Price during isotope analysis. The critical review by M. Hobbs is greatly appreciated. Thanks to two anonymous reviewers for their constructive comments.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


## *Article* **Origin of Drusy Dolomite Cement in Permo-Triassic Dolostones, Northern United Arab Emirates**

**Howri Mansurbeg 1,2,\*, Mohammad Alsuwaidi 3, Shijun Dong 3, Salahadin Shahrokhi <sup>4</sup> and Sadoon Morad <sup>3</sup>**


**Abstract:** While the characteristics and origin of drusy calcite cement in carbonate deposits is well constrained in the literature, little attention is paid to drusy dolomite cement. Petrographic observations, stable isotopes, and fluid-inclusion microthermometry suggest that drusy dolomite cement in Permo-Triassic conglomerate/breccia dolostone beds in northern United Arab Emirates has precipitated as cement and not by dolomitization of drusy calcite cement. The low δ18OVPDB (−9.4‰ to −6.2‰) and high homogenization temperatures of fluid inclusions in drusy dolomite (Th = 73–233 ◦C) suggest that dolomitization was caused by hot basinal brines (salinity = 23.4 wt% NaCl eq.). The δ13CVPDB values (+0.18‰ to +1.6‰) and 87Sr/86Sr ratio (0.708106 to 0.708147) indicate that carbon and strontium were derived from the host marine Permo-Triassic carbonates. Following this dolomitization event, blocky calcite (Th = 148 ◦C; salinity = 20.8 wt% NaCl eq.) precipitated from the hot basinal brines. Unravelling the origin of drusy dolomite cement has important implications for accurate construction of paragenetic sequences in carbonate rocks and decipher the origin and chemistry of diagenetic waters in sedimentary basins.

**Keywords:** drusy dolomite; dolomitization; cementation; hot basinal fluids

#### **1. Introduction**

Calcite cement with a drusy mosaic texture, displaying crystal size growth from pore wall to pore center, is widely distributed in limestones [1,2]. Such cement commonly fills moldic pore and hence has been interpreted to be sourced by dissolution of aragonitic allochems [3–7]. The stable isotopic composition of drusy mosaic calcite cement coupled with the presence of moldic porosity have been interpreted to indicate a meteoric origin [8–10]. Dolomite cement in carbonate rocks commonly develops rhombic crystals that do not show systematic size variation within the pores [11]. Nevertheless, rare drusy mosaic dolomite cement has been reported in carbonate deposits too [12–15]. The origin of such dolomite is not fully explored, particularly regarding whether it is a primary precipitate or formed by dolomitization of precursor drusy calcite cement. The fluid flow and related diagenetic evolution, particularly dolomitization of carbonate succession in northern United Arab Emirates (UAE), is interpreted to be controlled by the tectonic evolution of the region [16,17].

In this paper, we constrain the origin of drusy mosaic dolomite cement in conglomerates/breccia of the Bih Formation (Permo-Triassic), Ras Al Khaima, northern UAE (Figure 1), using petrographic, stable isotopic, and fluid-inclusion microthermometric analyses. This study provides important insights into the conditions of formation of drusy dolomite cement and its associated diagenetic fluids. The presented data and discussions

**Citation:** Mansurbeg, H.; Alsuwaidi, M.; Dong, S.; Shahrokhi, S.; Morad, S. Origin of Drusy Dolomite Cement in Permo-Triassic Dolostones, Northern United Arab Emirates. *Water* **2021**, *13*, 1908. https://doi.org/10.3390/ w13141908

Academic Editor: Kristine Walraevens

Received: 4 June 2021 Accepted: 6 July 2021 Published: 9 July 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

<sup>1</sup> General Directorate of Scientific Research Center, Salahaddin University-Erbil, Erbil 44001, The Kurdistan Region, Iraq

help in constructing more accurate paragenetic sequences in Permo-Triassic dolostones, northern UAE, and in other carbonate successions with similar depositional and diagenetic settings.

**Figure 1.** Location of the study area, northwest of the UAE (marked in green). Modified from [18–20].

#### **2. Geological Setting**

The Bih Formation (Permo-Triassic) crops out in the Musandam Peninsula, Ras Al Khaimah, northern United Arab Emirates (Figure 1). The formation is suggested to be age-equivalent to the Khuff Formation [19] which is an important hydrocarbon reservoir in the Arabian Gulf region. Deposition of the Bih Formation occurred in subtidal to intertidal (shoals, lagoon, and mudflats) environments on a passive margin of the Arabian plate [20].

The formation is divided into two units separated by a marker interval, which is called mid-Bih conglomerates/breccias [21]. The origin of this interval is uncertain, being interpreted to have resulted from the collapse of evaporite beds [20] tectonic deformation [20] or are lag deposits formed owing to marine transgression [22]. Strohmenger et al. [20] proposed that the Bih Formation corresponds to sequence KS4 through KS7 of the Khuff Formation (Figure 2). Two main tectonic events have affected the study area [23] including: (i) The obduction of Oman ophiolite on the Arabian platform during the Late Cretaceous. During this tectonic phase, the passive northeastern margin of the Tethys ocean became

compressive and resulted in the close of the Tethys ocean [24]; (ii) the Zagros orogeny (Paleocene to Middle Miocene) [25–28] which caused large scale folding and thrusting [26].

**Figure 2.** (**A**) Chronostratigraphic framework of the Khuff Formation in the subsurface compared to the stratigraphic range of the Bih Formation [21]. (**B**) Logging graph showing the thickness of conglomerate/breccia beds in the field.

#### **3. Samples and Methods**

Fifty-six samples were collected from the conglomerate/breccia beds and associated dolostone beds. Thin sections were prepared for all the samples after impregnation with blue epoxy. Representative polished uncovered thin sections were prepared and examined by a Technosyn Cold cathodoluminescence microscope (CL), backscattered electron imaging (BSEI), and energy dispersive X-ray analyses (EDS) attached to a Quanta scanning electron microscope. Fluid-inclusion microthermometry was conducted on two representative double polished 50 μm thick sections using a pre-calibrated LINKAM THMS600G stage fitted onto the Nikon E600 microscope.

Synthetic pure H2O and CO2 inclusions standard were used to calibrate the thermocouple. Both homogenization (Th) and final ice melting temperatures (Tmice) were measured whenever the size of the inclusions permitted. The former is the minimum entrapment temperature of the inclusion (i.e., the temperature of mineral precipitation) [29]. The latter is the temperature at which the last ice melts in an aqueous fluid inclusion and reflects the trapped diagenetic fluid salinity [29]. The fluid salinity (wt% NaCl eq.) was inferred from Tmice using the Bodnar equation [30]. In order to avoid overheating, low heating rates were used for measuring homogenization temperatures (Th).

Stable carbon and oxygen isotopes were obtained for the 36 micro-drilled drusy dolomite cement with a GVI IsoPrime continuous flow isotope ratio mass spectrometer system (CF-IRMS). The δ13C and δ18O isotopic values are determined on carbon dioxide (CO2) from carbonate minerals by the reaction with 100% phosphoric acid (H3PO4) at 90 ◦C in a vacuum using standard procedures based on the principles reported in [31]. Samples and standards are weighed into Exetainer™ septum vials, sealed and flushed with pure helium. In addition, 100% phosphoric acid is injected into the vial, which is then placed in an aluminum tray maintained at 90 ◦C for more than 1 h. The CO2 is extracted automatically with a double-hole needle Gilson auto-sampler connected to a GVI MultiFlow. The MultiFlow contains a 500-uL sample loop, and a GC column that separates CO2, N2, O2, and H2O. Helium carrier gas then transports the purified CO2 into the GVI IsoPrime continuous flow isotope ratio mass spectrometer system (CF-IRMS) that measures the isotope ratios. Four representative micro-drilled drusy dolomite cement samples were analyzed for Sr isotopes. One mg of each sample is collected to react with suprapure HCl for 24 h. Purified Sr is obtained by chromatographic separation with 2.5 mL of AGW 50 × 8 (Biorad) cation exchange resin [32]. The analysis equipment is a Finnigan MAT 262 7-collector solid-source mass spectrometer with a single Re filament applying 1 μL of ionization enhancing solution [32].

#### **4. Results**

#### *4.1. Field Observations*

The base of the Bih Formation does not crop out, and the lowest stratigraphic part consists of a medium-to-thick-bedded brown dolomitized packstone–grainstones with locally well-developed coarse-crystalline dolostone and vuggy textures. About 180 m above the base of the exposed Bih Formation, the mid-Bih conglomerate/breccia occurs (approx. 20 m thick; Figures 2B and 3A). The conglomerate/breccia beds of Wadi Bih, which are interbedded with dolograinstones, are massive overall but locally show a low-angle, through cross stratification (Figure 3B). The beds are cross cut by a series of normal faults (Figure 3A). The basal contact between conglomerate/breccia and underlying dolograinstone beds is sharp and erosional (Figure 3C). The lower contact between dolograinstone beds with the conglomerate/breccia beds is marked by the presence of rip-up dolomudstone intraclasts (Figure 3D).

Stylolites occur within as well as along boundaries between the conglomerate/dolostone beds. Above the mid-Bih breccia, there are cycles of interbedded dolostones and dolomitic limestones, with algal laminations, bird-eye structures, and dissolution molds. The succession continues with 90 m cliff-forming, dark gray coarse-crystalline dolostones.

**Figure 3.** Field photographs showing: (**A**) The conglomerate/breccia beds are cut across by a series of normal faults, the conglomerate/breccia beds are shown by white lines, the faults are shown by orange lines. (**B**) Parallel lamination occurs at 14 cm above the basal contact in the conglomerate/breccia beds. (**C**) The contact (yellow line) between conglomerate/breccia and underlying massive dolostone bed. (**D**) The conglomerate/breccia beds have a gradational contact with overlying dolograinstone bed (arrows).

#### *4.2. Petrography of the Drusy Dolomite*

The conglomerate/breccia dolostone beds are cemented by drusy mosaic dolomite with a crystal size increasing progressively from 20 μm along the pore wall to 2 mm in the pore center (Figure 4A–D). In some cases, the centermost parts of the pore are occupied by blocky calcite (500–700 μm across). Dolomite crystals adjacent to this blocky calcite are partly calcitized (Figure 5A,B). The drusy dolomite crystals display sector CL zoning with alternating non-luminescent and bright orange luminescent zones (Figure 5C,D). The microcrystalline equant to bladed dolomite crystals (20–100 μm across) along the pore walls are dull luminescent. The blocky calcite crystals are dull luminescent with a local thin orange luminescent zone.

Dull luminescent calcite occurs along stylolites and fills fractures varying in width from 5 to 60 μm (Figure 5E,F). The stylolites (amplitudes up to 2 cm) and dissolution seams cut across the dolostone pebbles and have bed-parallel and bed-oblique orientations. Some of the stylolites are cutting across both pebbles and dolomite between pebbles (Figure 6A–D). The drusy mosaic dolomite has Fe and Mn contents below the EDS detection limits. Fe and Mn contents as low as few ppm can still generate zoning in dolomite [33–35].

**Figure 4.** Optical photomicrographs (PPL) showing: (**A**) The fine dolomite crystals (arrows) in mosaic dolomite between pebbles are considered to be relics of fine-crystalline dolostone pebbles. The distribution of the fine dolomite crystals shows the rounded shape of pebble (dashed line). (**B**) The fine-crystalline dolostone pebble is partly replaced by the medium-crystalline dolomite (40 μm). Parts (arrows) of the pebble are not subjected to replacement. (**C**) Fracture cutting across the mosaic dolomite cement between pebbles was filled by the mosaic calcite cement (5 to 60 μm). The average width of this fracture is 70 μm. (**D**) Fracture-filling calcite along a pebble.

**Figure 5.** *Cont.*

**Figure 5.** Optical (**A**; PPL) and companion CL (**B**) photomicrographs showing several generations of dolomite cement, which filled intergranular pores between pebbles are characterized by a distinguishable luminescent degree. Optical (**C**; PPL) and companion CL (**D**) photomicrographs showing: Coarse mosaic dolomite cement and coarse calcite cement, which is reddish stained in intergranular pores between pebbles; concentric growth pattern of planar-e to planar-s mosaic dolomite and non-luminescent calcite cements precipitated in pore space, where precursor dolomite was dissolved. Four generations of dolomite cement have been identified in the companion CL (**D**) photomicrograph. Optical (**E**; PPL) and companion CL photomicrographs (**F**) showing the calcite precipitated in intercrystalline pores between pebbles (white arrow) and within stylolite (yellow arrow) is characterized by a non-luminescent zone.

**Figure 6.** BSE images showing: (**A**) The presence of dolomite along a dissolution seam within a pebble. The core of rhombic dolomite crystals is dissolved while the rim is preserved. (**B**) Calcite cement (Ca) along the dissolution seam within the pebble, which precipitates after dissolution of dolomite (shown by an irregular boundary of dolomite crystal, yellow arrows). (**C**,**D**) stylolite/seams (S1 and S2) cutting across each other and the clay mineral and TiO2 occur along the stylolite/solution seam within a pebble.

#### *4.3. C, O, and Sr Isotopic Composition of Dolomite Cement*

The δ13CVPDB values of drusy dolomite vary from +0.2‰ to +1.6‰ and δ18OVPDB from −9.4‰ to −6.2‰. The <sup>δ</sup>13CVPDB values of fracture-filling calcite cement vary between −3.8‰ and +1.9‰ and <sup>δ</sup>18OVPDB from −7.6‰ to −4.0‰. The 87Sr/86Sr ratios of two drusy and mosaic dolomite samples are 0.708106 and 0.708147, respectively (Figures 7 and 8A,B).

**Figure 7.** Cross plot of stable δ13CVPDB versus δ18OVPDB values of samples from the conglomerate/breccia beds and the dolograinstone beds in the conglomerate/breccia interval, as well as isotope values of Permian to Triassic marine carbonates. Stable oxygen isotope values of drusy dolomite are lower compared to other dolomites and slightly lower than sea coeval marine carbonates.

#### *4.4. Fluid-Inclusion Microthermometry of Carbonate Cement*

Fluid inclusions (1–20 μm across) in the drusy mosaic dolomite cement display twophases (liquid and vapor). The crystal sizes of drusy dolomite are divided into fine (50–100 μm), medium (>100–500 μm), and coarse crystals (>500 μm). Results can be summarized as follows: (1) Homogenization temperatures (Th) of fine-sized dolomite crystals (50–100 μm), range from 93 to 143 ◦C (av. 118 ◦C); (2) Th of medium-sized dolomite crystals (100–500 μm) varies from 73 to 188 ◦C (av. 130.5 ◦C). The salinity of these fluid inclusions could not be measured owing to the small size of the inclusions (<2–3 μm); (3) Th values of coarse crystalline dolomite (>500 μm) filling the centermost part of the pores vary from 118 to 233 ◦C (av. 175.5 ◦C), while salinity ranges from 21.8 to 25.1 wt% NaCl eq.; and (4) the calcite cement that postdates dolomite has a Th value of 148 ◦C and salinity of 20.8 wt% NaCl eq. The distribution of Th in calcite and drusy dolomite is shown in histograms (Figure 9).

**Figure 8.** (**A**) Cross-plot of 87Sr/86Sr against δ18OVPDB showing two samples from the pebble have a slightly higher Sr isotope value than coeval marine carbonates. Two samples from drusy and mosaic dolomites are within the isotopic range of Permian to Triassic marine carbonates. The sample from drusy dolomite showing a lower δ18OVPDB value than coeval marine carbonates. (**B**) Measured 87Sr/86Sr ratios of different types of dolomite cements in the Bih Formation (Permo-Triassic) plotted on the sea water 87Sr/86Sr curve [31].

**Figure 9.** Histograms of fluid-inclusion homogenization temperatures (Th) in fine-crystalline dolomite (**A**), mediumcrystalline dolomite (**B**), coarse-crystalline dolomite (**C**), and calcite (**D**).

#### **5. Discussion**

Constraining the origin of the drusy dolomite in terms of whether it has been precipitated as cement, which is a common cement in limestones [13,14] or formed by the replacement of precursor drusy calcite [15] has implications for unravelling the role of diagenetic paleofluids, as well as for improvements in the definition of paragenetic sequences in carbonate rocks and related reservoir-quality features. There are two possible interpretations of the formation of drusy dolomite cement: (1) Precipitation as cement, and (2) fabric-preserving dolomitization of precursor drusy mosaic calcite. Evidence supporting precipitation as drusy dolomite includes the sector cathodoluminescence zoning and its close association with the latter precipitated blocky calcite.

Fluid-inclusion microthermometry of drusy mosaic dolomite and the subsequently precipitated calcite cement indicate formation from hot brines during compressional tectonic events in conjunction with the obduction of Oman ophiolites [16,36]. However, it is not immediately clear whether the drusy mosaic dolomite has been precipitated as cement [13] or formed by dolomitization of drusy mosaic calcite by the hot basinal brines [15,37]. Drusy calcite is a common cement in limestones [1,5–7,38]. Thus, the common presence of drusy mosaic dolomite in totally dolomitized limestones may suggest formation by fabric-preserving dolomitization [15]. The latter process suggests similar dissolution rates of drusy calcite and precipitation rate of dolomite during the dolomitization process. A similar fabric-preserving dolomite of Bih Formation has been reported in earlier studies of carbonate reservoirs in the world [39–41].

The high Th (73 to 233 ◦C) and salinity (20.8 to 25.1 wt% NaCl eq.) of drusy mosaic dolomite and subsequently precipitated cement indicates that they have precipitated either during deep-burial diagenesis or by upward flux of hydrothermal/hot basinal fluids [42]. The increase in Th with the increase in dolomite crystal size (av. 118 to 175.5 ◦C) is interpreted to indicate that the drusy dolomite precipitated with the increase in temperature during progressive burial. The maximum burial temperatures reached by the Bih Formation is approximately 190 to 200 ◦C [22]. The wide range of Th of different carbonate phases indicates that some fluid inclusions have been subjected to re-equilibration during burial processes [29]. The relatively systematic Th zoning from pore walls inwards (i.e., increase in crystal size) supports the formation of drusy dolomite by precipitation from pore fluids than by dolomitization of precursor drusy calcite cement. Apparently, the chemistry of the dolomitizing fluids has changed with time and reflected on the variation of Fe-rich and Fe-poor zoning of CL pattern of the drusy dolomite cement (cf., [2,36]).

The similar Th (148 ◦C) and salinity values (20.8 wt% NaCl eq.) of late blocky calcite and fracture-filling medium to coarse calcite cement is attributed to the circulation of hot brines [17]. The considerably lower <sup>δ</sup>13CVPDB values (−3.8‰ and −0.9‰) of fracture-filling calcite cement than drusy dolomite samples might indicate the deviation of dissolved carbon from degradation of hydrocarbons or organic matter [42]. Using the δ18OVPDB values (−9.4‰ to −6.2‰) and Th (73 to 233 ◦C) of drusy dolomite, the fractionation equation of Land [43] it is inferred that precipitation has occurred from evolved brines (δ18OVSMOW = +0.29‰ to +10.96‰). The wide extent of variations in the degree of geochemical evolution of the brines could be related to variable fluid sources and variable degrees of interaction with the host rocks in the basin. The positive correlation between C and O isotopes of drusy mosaic dolomite (Figure 7) suggests increasing input of 12C into the formation waters with the increasing temperature [44]. The small variations in Sr isotopic compositions between the various types of carbonate cements suggest a related origin (Figure 8A,B). The 87Sr/86Sr ratios in one of the pebbles (0.708202) is within the isotopic range of Permian to Triassic seawater values (0.7069 to 0.7083) [31]. The slightly higher ratio in another pebble sample (0.708424) than the range for Permo-Triassic seawater (Figure 8B) suggests that the hot basinal fluids have interacted with siliciclastic rocks [45,46]. The similar carbon isotope and 87Sr/86Sr ratio of drusy dolomite and coeval marine carbonates ([31]; 0.708106 and 0.708147, respectively) suggests that dissolved carbon and Sr were derived from the dissolution of marine carbonates and/or marine pore waters. The slightly lower oxygen isotope value of drusy dolomite than coeval marine carbonates (Figure 8A) is attributed to precipitation at elevated temperatures.

The development of stylolites and dissolution seams along the rim of pebbles or cutting across both pebbles and dolomite between pebbles indicate that stylolitization postdates the drusy dolomite cement. This complex relationship between fracturing, stylolitization, and dolomite cementation may be caused by multiple phases of tectonic movements. The cross-cutting of low amplitude seams by high amplitude stylolites in conglomerate/breccia samples could be a result of increasing compressive stress by tectonic movements generated during the obduction of Oman ophiolites [47,48]. The presence of dolomite cement (Th up to 190 ◦C; salinity up to 20.8 wt% NaCl eq.) along the stylolites suggests the flux of hot basinal brines and may indicate that the dolomite was precipitated by dolomitization of the host limestones prior to stylolitization. Dolomite cement apparently has been concentrated in the vicinity of the stylolitization surface as an insoluble residue (cf., [49]). Similar observations and interpretations have been presented for carbonate successions elsewhere [49–51]. The presence of clay minerals (locally pigmented by TiO2) along stylolites suggests that these minerals promoted the pressure dissolution of carbonates [52,53].

Constraining the depths at which the flux of hot dolomitizing basinal brines occurred is difficult without radiometric dating of dolomite [54,55]. However, the cross-cutting relationship between the stylolites and dolomite cement suggests that dolomitization took place at shallow depths [36].

#### **6. Conclusions**

The puzzling origin of pore-filling drusy mosaic dolomite cement in Permo-Triassic conglomerate/breccia beds outcropping in northern United Arab Emirates is constrained using the integrated field, petrographic, isotopic, and fluid-inclusion microthermometric analyses. There are two possible interpretations of the formation of drusy dolomite cement: (1) Precipitation as cement, and (2) fabric-preserving dolomitization of precursor drusy mosaic calcite. Fluid-inclusion microthermometry indicates the formation of this dolomite from hot basinal brines (Th = 73 to 233 ◦C; salinity = 21.8% to 25.1% wt% NaCl eq.). Using the average homogenization temperature and oxygen isotope of dolomite δ18OVPDB (−9.4‰ to −6.2‰) the brines are inferred to be geochemically evolved with <sup>δ</sup>18OVSMOW (+0.29‰ to +10.96‰). Lines of evidence supporting precipitation as drusy dolomite from pore waters include the sector cathodoluminescence zoning, the relatively systematic Th zoning, and with the presence of the latter precipitated blocky calcite, which has also precipitated from hot basinal brines (Th = 148 ◦C salinity = 20.8 wt% NaCl eq.). The construction of accurate paragenetic sequences is important for establishing exploration and production models for carbonate reservoirs. To achieve this goal, diagenesis should be linked to fluid history and tectonic evolution of the basin.

**Author Contributions:** H.M.: Conceptualization, Writing—original draft. M.A.: Writing—review & editing. S.D.: Formal analysis, Methodology. S.S.: Writing—review & editing, Visualization. S.M.: Writing—review & editing, Supervision, Resources, Investigation. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Not applicable.

**Acknowledgments:** The fluid inclusion analyses were performed by Daniel Morad at the French Petroleum Institute, Paris, France, the carbon and oxygen isotope analyses at the Environmental Isotope laboratory, University of Waterloo, Canada, and the strontium isotope analyses at Bochum University, Germany.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**

