Porosity characterization using medical computed tomography scans of two deeply buried paleokarst aquifers in Silurian carbonates, southern Ontario, Canada

Paleokarst and paleokarst aquifers are not as well-documented as shallow karst, and studies of pores and pore networks are rare. Paleokarst lacks the open conduits typical of karst due to compaction, infilling, and diagenetic recrystallization by deep burial, so groundwater movement is through matrix porosity. In this case study, porosity networks in two saline aquifers in carbonate paleokarst of the Silurian Guelph Formation and Salina A-1 Carbonate Unit in southern Ontario (Canada) have been studied in drill cores at microscopic to macroscopic scales, utilizing medical computed tomography (CT) scans, optical petrography, macroscopic core examination, and scanning electron microscopy. The CT scans provided nondestructive three-dimensional visualization and quantification of pore distribution, size and volume, pore connections, and estimates of total porosity, similar to gas porosimetry. In the A-1 Carbonate, 50% of the pore volume comprises layers of macropores associated with algal laminations, with good horizontal connections. In the Guelph Formation, vuggy macropores contribute most of the pore volume. They are connected through subvertical fractures and rubble-filled karst conduits and by abundant small pores and micropores with poor horizontal connections. The Guelph paleokarst represents a longer period of subaerial exposure as evidenced by its greater thickness and geographic extent, former karst conduits, and near-total destruction of primary sedimentary fabrics. The findings provide real-world examples for construction of laboratory scale models of paleokarst aquifers, demonstrate the value of a multi-scale approach to porosity studies, and showcase the value of medical CT scans with its unique ability to visualize pore connections within drill cores.


Introduction
As defined by Simms (2014), paleokarst is "karst that has been buried by younger rocks", and in the context of this study, is karst that formed in the geologic past during an extended period of subaerial exposure of carbonate (and evaporite) bedrock with subsequent burial and modification in the deep subsurface.Karst and paleokarst both enhance the porosity and permeability of the affected rocks and their ability to store and transport groundwater.Karst and karst aquifers are widespread and well-documented in sedimentary basins (e.g.Bakalowicz 2005; Ford and Williams 2007; Chen et al. 2017).Paleokarst and paleokarst aquifers are not as well documented and likely have often not been recognized.Recognition is complicated by compaction, infill of karst conduits by rubble, and deep burial diagenetic recrystallization.As noted by Simms (2014), any sedimentary basin with significant accumulations of carbonate rocks is likely to have paleokarst in the intermediate to deep subsurface, and paleokarst aquifers.
Regional unconformities representing extended periods of subaerial exposure are common in the Paleozoic bedrock of southern Ontario, Canada, and where they occur in carbonate rocks, porous and permeable paleokarst horizons are often present, forming regional confined aquifers (Carter et al. 2021a).In the intermediate to deep subsurface (>100 m), they contain brackish to saline sulphur water or brine with up-dip connections to the shallow freshwater regime.They are utilized for disposal of oilfield brine, as a source of brine for ice-control on highways and are being considered for potential injection of carbon dioxide for carbon capture and storage (Carter et al. 2021a).In the geologic past, these paleokarst horizons provided pathways for regional migration of hydrocarbons into subsurface hydrocarbon traps.In the present day, hydrochemical and isotopic zonation of groundwater in these aquifers are important considerations in studies to establish a safety case for siting of deep geologic repositories for long-term isolation of nuclear wastes (Hobbs et al. 2011).Despite their importance, porosity development and the formation of connected pore systems has not been studied in detail in any of these paleokarst aquifers.
The subject of this study is the Guelph Aquifer in Ontario, a confined brine aquifer (Carter et al. 2021a) in a regional paleokarst horizon in dolomitized carbonates of the Guelph Formation, and the recently documented subregional A-1 Carbonate Aquifer developed in massive to algal-laminated lime mudstone (Figs. 1 and 2).Paleokarst has been documented in the Guelph Formation in several studies (Smith 1990;Carter et al. 1994;Brunton and Brintnell 2020;Carter et al. 2021a).These same studies documented paleokarst in the overlying limestones and dolostones of the A-1 Carbonate Unit localized in paleotopographic highs above pinnacle structures.The subregional paleokarst of the A-1 Carbonate was only recently identified by Carter et al. (2021a) at the location of the present study.This study focuses on the Bruce nuclear site on the eastern shore of Lake Huron (Fig. 1).
Microscopic studies of matrix porosity in karstic carbonates are rarely completed as water movement in modern shallow karst is principally controlled by large-scale (cm to m) conduits of solution-widened fractures and joints (e.g.Bakalowicz 2005; Worthington 1999).The Guelph and A-1 Carbonate paleokarst have been thoroughly recrystallized during burial diagenesis in the deep subsurface with mechanical compaction and diagenetic infill of open conduits and now have hydrogeologic characteristics of a porous media, presenting a unique opportunity to characterize the matrix porosity of an ancient karst interval(s).This study was made possible by the availability of a full-diameter drill core from closely spaced deep research boreholes completed at the Bruce nuclear site to assess its suitability as a deep geological repository for nuclear wastes.The present work employs a multiscale approach using a combination of medical computed tomography (med-CT) scan technology, macroscopic core examination, gas porosimetry, optical petrography, and scanning electron microscopy (SEM) analysis.These direct observation methods are aimed at bridging the gap between well log observations and core plug studies and attempts to optimize the representativeness of geological sampling in subsurface studies.The pore fabric of the two aquifers is also compared to assess differences in pore type, sizes and connectivity related to (1) the different carbonate substrates, and (2) maturity of karst development.Mohammedi et al. (2021) review modern approaches to constructing laboratory scale models of karst aquifers that are representative of actual field conditions.The present study addresses the general lack of real-world examples for construction of physical models for paleokarst aquifers.The methods employed here could have widespread application to deep saline aquifers in paleokarst throughout southern Ontario and within other sedimentary basins.

Geological setting
Southern Ontario is underlain by up to 1,500 m of Paleozoic marine sedimentary rocks comprised principally of carbonates and shales, with thick evaporites in the Michigan Basin and sandstones in the Appalachian Basin (Armstrong and Carter 2010;Carter et al. 2021b).Bedrock formations straddle a regional structural arch, dipping shallowly into the Michigan and Appalachian basins and forming northwesttrending subcrop belts (Fig. 1).These rocks were subject to many periods of subaerial exposure and erosion as evidenced by the presence of numerous disconformities, resulting in the formation of several paleokarst horizons in the carbonate rocks.The formations of interest to this study are the Lower Silurian Guelph Formation and the Upper Silurian A-1 Carbonate Unit of the Salina Group (Fig. 1).The two units underlie most of southern Ontario west of the Niagara Escarpment, and at the Bruce nuclear site occur between 326 and 390 m below the ground surface (bgs) (Intera 2011).The upper contact of the Guelph Formation is a regional unconformity (Brunton and Brintnell 2020;Carter et al. 2021a).
Regionally, the Guelph Formation is comprised of dolostones and limestones deposited on an eastward deepening, shallowing-upwards carbonate ramp (Brunton andBrintnell 2011, 2020).The Bruce site is located within the inter-pinnacle lithofacies belt of the Guelph Formation (Fig. 2).Within the inter-pinnacle lithofacies the Guelph paleokarst aquifer is composed of 2-8 m of a distinctive porous and permeable breccia identified as a paleokarst rubble/paleosol breccia (Smith 1984(Smith , 1990;;Smith et al. 1993;Carter et al. 1994), resulting from two or more periods of subaerial exposure over an extended period of time in the geologic past.These exposure events resulted in thinning of the Guelph Formation and obliteration of primary depositional features and fabrics, extending downward into the argillaceous lime wackestones of the uppermost Goat Island Formation.Subsequent burial beneath the evaporitic Salina Group preserved these porous and permeable carbonates in a regional paleokarst horizon (Fig. 2), forming the regional Guelph Aquifer (Carter et al. 2021a) underlying a large part of southern Ontario.Two water samples from the Bruce nuclear site contained total dissolved solids (TDS) concentrations of 365,604 and 375,468 mg/L, respectively (Intera 2011).The average measured horizontal hydraulic conductivity in the unit is 3 × 10 -8 m/s with an up-dip hydraulic gradient (Intera 2011).
At the study site, the Guelph paleokarst is gradationally overlain by a series of salinification-upwards carbonate-anhydrite-halite cycles comprising the Salina Group (Carter et al. 1994;Armstrong and Carter 2010).Thick beds of halite occur in the central Michigan Basin, thinning easterly into Ontario with lateral facies change to anhydrite, dolomite and increasing amounts of shale (Armstrong and Carter 2010) indicating transition to a shallow sabkha depositional environment.The A-1 Carbonate is the basal unit of the second carbonate-evaporite cycle of the Salina Group (Carter et al. 1994).The Salina Group, including the A-1 Carbonate, forms a thick regional aquitard beneath much of southern Ontario west of the Niagara Escarpment (Carter et al. 2021a).At the study site the lower (nonkarstic) A-1 Carbonate has horizontal hydraulic conductivity of only 9 × 10 -12 m/s (Intera 2011).
Through most of southern Ontario the A-1 Carbonate is a regional aquitard but at the Bruce site an anomalous water-bearing porous and permeable paleokarst interval occurs within the uppermost 3-4 m of the formation, at a depth of approximately 340 m bgs (Fig. 1).The subsurface extent of the paleokarst aquifer is delineated by water intervals encountered during the drilling of petroleum wells (Fig. 2b) up to 40 km downdip from the subcrop belt.Horizontal hydraulic conductivity of the paleokarst in deep boreholes at the Bruce site averaged 2 × 10 -7 m/s (Intera 2011).Two water samples contained total dissolved solids of 26,760 and 30,455 mg/L, interpreted to represent mixing of the original dense brine, typical of the Guelph Formation and the overlying Salina Group formations in the deep subsurface, with glacial meltwater flowing downdip during glaciation (Intera 2011;Carter et al. 2021a).
The Guelph Formation paleokarst represents a considerable period of exposure at a regional unconformity, whereas the A-1 Carbonate, with features typical of an epikarst (Klimchouk 2004), is hypothesized to represent a much shorter period of exposure on the shallow edges of the Michigan Basin.

Materials and methods
The study materials comprise 15 full-diameter samples of the diamond drill core of the Guelph Formation and the Salina A-1 Carbonate Unit, obtained from four boreholes: DGR-1, DGR-3, DGR-4, and DGR-8.The boreholes were drilled at the Bruce nuclear site between 2007 and 2011 (Fig. 3) to assess the feasibility of long-term disposal of low to intermediate level nuclear waste in a deep geological repository (Raven et al. 2010; Nuclear Waste Management Organization 2011).The boreholes were continuously cored from the bedrock surface, reaching total depths varying from 465.1 to 871.3 m.Depth intervals range from 326 to 339 m for the A-1 Carbonate samples and 380-390 m for the Guelph Formation samples.Core diameters are 75 mm with individual sample lengths from 10 to 60 cm.Study methods include macroscopic core examination, geophysical log interpretation, med-CT scan imagery, full-diameter core analyses of porosity and permeability, and petrographic and SEM studies (Fig. 4).High-quality hydraulic conductivity data is also available from other studies (Intera 2011).

Drill core logging
Lithological core logging included the complete cored interval of the Goat Island and Guelph formations and the A-0 Carbonate, A-1 Evaporite and A-1 Carbonate units of the Salina Group.Detailed logging focussed on carbonate rock types (Dunham 1962), porosity classification (Choquette and Pray 1970) and diagenetic features including dolomitization, dissolution, cementation and karstification.

Petrography and SEM
Five thin sections have been prepared for carbonate microfacies analysis using criteria listed in Flugel (2010).The samples were cut parallel to core axis and impregnated with blue epoxy resin to highlight pore spaces.Sample impregnation and thin section preparation were completed by Precision Petrographics Ltd. in Langley, BC.Petrographic analysis was based on examination under transmitted light using Leica Macroscope Z16 APO in the Surface Science Western Laboratory (SSWL) at University of Western Ontario.Thin section examination focused on grain and matrix types, porosity types and sizes, and diagenetic fabrics.Analysis of energy-dispersive X-ray (EDX) and porosity types of four polished thin sections was completed using the Hitachi SU3500 scanning electron microscope (SEM) with an energy-dispersive X-ray (EDX) detector at SSWL.The specimens were vacuum dried for 20 min and coated with gold using an Osmium Plasma Coater, which increases the conductivity of the samples and reduces damage from electron beam when being imaged.Images were made in backscattered electron (BSE) mode (Goldstein et al. 2018).

Conventional core analysis
Four samples representative of the paleokarst intervals were collected by the authors from archived drill core and sent for whole core porosity and permeability measurements at AGAT Laboratories, Calgary, Canada.Samples were dried in a convection oven for 96 h at 108 °C before measurement.Helium gas porosity was measured relying on Boyle's law.Grain volume (gv) was measured by the change in helium pressure in vacuum.Bulk volume (bv) was measured by calipering on right-cylindrical samples that were assumed to have a perfect cylindrical geometry.Porosity was calculated by the difference of grain volume and bulk volume divided by bulk volume: (bv-gv)/bv.Permeability, including horizontal maximum permeability (Kmax), horizontal minimum permeability (K90) and vertical permeability (Kv) were measured by the flow rate of nitrogen in steady state in all four fulldiameter samples.

Medical CT-scanning
Medical CT-scanners measure X-ray attenuation inside a sample using a rotating X-ray source and detector.The result is a three-dimensional (3D) volume of the sample expressed in Hounsfield units (HU).X-ray attenuation, hence, HU values, is a function of density, material (atomic number for elements or the effective atomic number for compounds) and size.While it is difficult to calculate a precise density using CT-scanning with only one energy level, changes in density are routinely observed.Density can be assessed using dual-energy CT scanning (DECT) (Wellington and Vinegar 1987;Siddiqui and Khamees, 2004;Al-Owihan et al. 2014) but, although promising, it requires higher CT-scanning times and computational power (therefore it is more expensive), and, still presents several challenges related to the choice of the calibration method, the beam hardening corrections to apply or the complexity of managing image noise and resolution for very heterogeneous samples (Victor et al. 2017;Martini et al. 2021).

Acquisition and data processing
CT-scan analyses were completed on 15 full-diameter core samples from the 4 boreholes (Fig. 4) at the Institute national de la recherche scientifique, Eau Terre Environment (INRS-ETE), Quebec, with a Siemens SOMATOM Definition AS+ 128.The increment between each CT slice was set to 0.1 mm with a 0.6 mm thickness.Images were recorded in DICOM format compatible with the open-source Fiji (Schindelin et al. 2012) or commercial 3D visualization software such as AVIZO or Dragonfly.Protocol and data analysis were adopted from Larmagnat et al. (2019), the main acquisition and reconstruction parameters are summarized in Table 1.

Quantitative analyses
Porosity was derived from CT-scan images using two methods: (1) a dry-state, and (2) two-state fully saturated (Larmagnat et al. 2019).In the dry-state analysis, every voxel is considered to be a mixel, i.e. a mixture of porosity (air) and solid phase material.The value of one mixel, in Hounsfield units (HU) is therefore an average value of its content.The porosity evaluation for the conventional mixel approach of dry state images has to be taken with caution because it assumes a constant and known mineralogy throughout the sample.With consideration of this limitation, the porosity estimate provides a proxy value of the total sample porosity.Subsequently, the second two-state protocol involved saturation of the core with water and displacement of air from connected pores.The connected porosity calculation assumes each voxel is a mixel and it is not influenced by mineralogy heterogeneity within the rock samples anymore, simply the change in density introduced by the water saturation, as saturated and dry volumes were subtracted.Thus, only the variation in density from the two states is left, and porosity is calculated by using that density variation.In addition to whole core porosity estimation, considering various two-dimensional (2D) and 3D views of each sample individually provides qualitative information such as macroporosity distribution vertically and horizontally, at the centimetric scale.Based on the scanning parameters (including voxel size) and the nature of the rock samples analysed, the spatial resolution limit for this study is considered to be 0.3 mm, meaning that the smallest pore diameter that can be fully resolved is 0.3 mm (at least nine voxels needed in the 2D slice, or 27 voxels in the 3D).
To gain more information on the macropores, porosity analysis was performed using the Aviso software from Thermo Fisher Scientific.It begins by applying a threshold on the porosity volume obtained from the subtraction of saturated and dry volumes.A threshold creates a binary volume by determining for every voxel whether its value is higher or lower than the threshold value, resulting in a 3D mask of the selected pore's shape-for example, applying a threshold of 5% to a porosity volume will result in a binary volume in which every voxel with a porosity value higher than 5% will now have a value of one, and zero for the rest.Morphological operations are then applied to the binary volume in order to smooth and clean the mask, with the most important one being the "opening operation".Such morphological noise-reduction operations are routinely performed on binarized data (e.g.Wildenschild et al. 2002).This operation along with the "remove small spots" operation will have the effect of removing very small structures (i.e.objects smaller than 0.92 mm 3 ), the reason being that they are indistinguishable from noise in the image.In most cases, multiple threshold values were attempted on the same volume as well as one or two opening values (Fig. 9).A label analysis is then applied to the mask in order to get information on the pores such as their size, position and average porosity value.Figures detailing pore-size class contribution to porosity are done using this data.Among them, an estimate of the pore diameter distribution with depth or the pore distribution histogram can be produced.Lastly, a 3D color-coded rendering of pores and a 3D color-scaled rendering of volume porosity are placed side by side to better understand the macropores' structure and porosity profile.The color-coded image is useful to see which pores are connected to one another as they will have the same color.The color-scaled render will instead give visual information on variations in porosity.The thresholding was applied to capture the macropore network, as only pores larger than 1 mm 3 remain in the segmented dataset.This documents the macropore network without giving unrealistic weight to microporous areas that are detected by the saturation processing but not fully resolved at the med-CT scale (Russ 2016;Kaestner et al. 2008).

Comparison of multi-scale methods
The different methods utilized here provide data on different aspects of the core samples examined, including pore sizes and connections and quantitative measurements of porosity or permeability (Table 2).The multi-scale capability ranges from ~100-0.001mm and ensures that the full range of pore sizes represented in the two paleokarst aquifers is captured and analyzed and the presence/absence of solution conduits can be recognized.Petrographic microscopy and SEM provide millimetric to near-nanometric scale, identification of grains, mineralogy, and porosity-occluding cements.CT scanning uniquely provides nondestructive 3D visualization of connected pores within the rock and the complexity of the pore networks they form (Van Geet et al. 2003;Mees et al. 2003) This is directly relevant to permeability.Conventional gas porosimetry is considered to be the industry standard for reservoir property studies (Anovitz et al. 2015;Njiekak et al. 2018).It provides a bulk volumetric estimate of connected porosity and permeability.Due to differences in resolution between med-CT and gas porosimetry the porosity estimation may be dissimilar; however, where macroporosity is significant and/or there is limited micropore connectivity, the difference will be minimal.These methods are complementary and together provide a complete understanding of the physical porosity network, how it was formed and its postburial modification.

Guelph Formation
The Guelph Formation unconformably overlies the Goat Island with an irregular karstic contact with infills of brown dolostone rubble (Fig. 5a).At the Bruce nuclear site, it is 3.54-5.35m thick.It is a dark brown to greyish brown breccia of pebble-sized karstic rubble dolostone clasts in a dolo-mudstone matrix.A medium grey, argillaceous dolo-mudstone 2.3-3.52-m-thickforms the lower Guelph Formation.Fossils are rare, whereas moldic porosity of gastropods and brachiopods are common in subvertical brownish karstic conduit infills in a medium grey dolo-mudstone matrix.Within the karstic conduit infills, dissolution-enhanced vugs and cavities are very common, partly plugged by dolomite cements.
The upper Guelph Formation comprises 1.24-1.83m of very porous, vuggy dolostone of recrystallized rubble clasts with a medium brown, dolo-mudstone matrix.Karstic conduits and cavities are common, partly plugged by dolomitic cements.There are no preserved fossils.Contacts of the rubble clasts are marked by irregular, dark grey residues with high organic content (e.g.Fig. 5a).The uppermost 30-50 cm is more karstic, with very common irregular vugs ranging from 1-3 cm in size, partly plugged by dolomite and halite.Convoluted stromatolites with dissolutionenhanced fenestral vugs occur as a 0.8-m-thick capping unit in well DGR-8.The Guelph Formation is sharply overlain by thinly laminated bituminous limestone of the A-0 Carbonate Unit.

A-1 Carbonate Unit
The A-1 Carbonate has a sharp contact with nodular to enterolithic anhydrite and minor dolostone of the underlying A-1 Evaporite.It consists of tan-grey to black, finely crystalline, variably organic-rich limestone and is 40.8-42.7 m thick with a lower member that is 0.8-15.2m thick, medium brown to black, oil-stained, lime mudstone with isolated vugs and cavities.The middle unit comprises 22.05-38 m of cyclic algal laminites, evaporitic algal laminites, enterolithic anhydrite and anhydrite nodules.Lenses or clusters of gypsum laths are very common.
The uppermost 1.73-3.58m comprises a very porous paleokarst of light brown to dark grey lime mudstone and crinkly laminated microbialites (Fig. 5b) with abundant fenestral vugs, moldic and intercrystalline pores.Stromatactic laminites occur with very abundant irregular, 1-3 cm vugs concentrated along and parallel with laminite layers.Dissolution-enhanced fenestral vugs 1-2.5 cm in diameter are very common, as well as subvertical fractures and leaching features.Oil stains are locally present.

Guelph Formation
The Guelph Formation was confirmed by energy dispersive x-ray spectroscopy (EDS) to be dolomitic limestone to dolostone.The uppermost Guelph Formation consists of karst rubble (D2) in a dolo-mudstone matrix (D1), cemented by finely crystalline dolomite crystals and euhedral dolomite cements.Oil stains accumulate along irregular seams at the contact of D1 and D2 (Fig. 6a-b).D1 dolomite is finely crystalline (<50 um), anhedral to subhedral, intergrown rhombic crystals with partially preserved grains and precursor textures.It forms porous sucrosic mosaics with planar crystal boundaries (Fig. 6a) or nonporous mosaics with irregular crystal boundaries (Fig. 6b).D2 dolomite is medium to coarsely crystalline (>50 um), loosely packed, sucrosic anhedral to euhedral rhombic crystals with either cloudy cores with clear rims (e.g.Fig. 6a) or tightly interlocking mosaics of anhedral nonplanar crystals (Fig. 6b).This type of dolomite is nonmimetic of mineral precursors.D2 has very common vugs and intercrystalline porosity.
Dolomite cement usually occurs as limpid dolomite and saddle dolomite crystals.Limpid cements are white, euhedral, finely to medium crystalline with rhombic terminations and uniform extinction.They constitute minor amounts of the dolomite cements lining molds and vugs.Coarse-crystalline euhedral saddle dolomite is the common lining or filling of irregular vugs (Fig. 6b).

A-1 Carbonate
The uppermost A-1 Carbonate Unit was confirmed by EDS analysis to be limestone.Both samples from the uppermost A-1 Carbonate consist of peloids and algal laminations cemented by fibrous calcites as rims of the grains (Fig. 6c).The dominant carbonate constituents are brownish grey, circular to horizontally elongate peloids, ranging in size from 0.2 to 1 mm oriented parallel to algal laminations or bedding planes (Fig. 6e).Pseudo-cortex of peloid grains resemble algae-originated ooids; however, these pseudo-cortices are formed by fibrous diagenetic calcite cement (Fig. 6f).Peloids are organic-rich with common oil stains (Fig. 6e,f).Fractures may occur in the peloid grains that form voids or are partly infilled by lime muds or calcite cements.EDS analysis indicates the fractures within peloids are also infilled with pyrite and fluorite with fenestral vugs showing oil staining.Irregular vugs are partially infilled with euhedral calcite crystals, white framboidal pyrite and have oil staining.Small crystals of halite infill the tiny pores in the limestone matrix.
Med-CT-images showed good contrast between pores, grains/matrix and dense mineral phases.This allowed the generation of 3D images at a macroscopic scale documenting pore geometry, connectivity and distribution that are not possible from petrographic studies inherently limited to the size and representability of thin sections (Fig. 7).The method delivers information either as total (dry state) porosity (conventional mixel approach) or connected porosity ('two states' saturation approach) that cannot be obtained using image analyses on thin sections or conventional gas porosimeter analyses.

Guelph Formation
Intercrystalline pores controlled by dissolution of calcite precursors and dolomite constitute >30% of the overall porosity and are well-developed in D1 and D2.In D1, intercrystalline porosity is <5 um, with larger interconnected pores (100 um) in D2.Some of the intercrystalline pores were solution-enhanced due to dissolution of dolomite (Fig. 6a-d).Intercrystalline porosity comprises smaller, irregular pores among calcite crystals that are variably plugged by halite.In the DRG-4 core at 375.08 m the D1 and D2 dolomites have an irregular contact, with intercrystalline pores that are enhanced by karstification and variably plugged by halite with an accumulation of oil stains (Fig. 6a,b).In the matrix of D2, intercrystalline pores are reduced by euhedral dolomite cements (Fig. 6b).Pin-point pores on the surface of the coarsely crystalline D2 are locally plugged by halite (Fig. 6d).
Moldic pores formed by selective dissolution of fossil fragments or crystals comprise 5-20% of total porosity.Pores are commonly <1mm in diameter (Fig. 6b).Moldic porosity is often reduced by saddle dolomite cementation that hinders observation of the dissolved crystal/fossil precursors.
Karstic vugs contribute >50% of the macro-porosity in the Guelph Formation and occur within both the karstic conduits and within karst rubble (Fig. 5a).Vugs are millimetre to centimetre-sized voids with irregular shapes and variable sizes.Many vugs could be solution-enlarged moldic pores (Fig. 6b).Med-CT imaging provides clear illustrations, at the pluri-centimetre scale, of karstic conduits that have oblique to subvertical contacts with the karst rubble (Fig. 7a-b).

A-1 Carbonate Unit
Channel porosity is the principal porosity type in the uppermost A-1 Carbonate Unit and forms horizontally wellconnected pore networks after dissolution of fibrous calcite cements above algal laminations.The channel porosity systems are 10s mm long, and 0.5-1 mm thick (Fig. 6e).CT imaging suggests the channel pores are coalesced intercrystalline pores and vugs (Fig. 7).Subvertical connections of the channel pore systems are located along degassing pathways or desiccation cracks.Fracture porosity identified by CT imagery is volumetrically minor, with partial plugging by calcite sparite, but may have played a major role in connecting pore systems.
Intra-particle porosity commonly occurs within the peloid grains (6f).Pores are usually <1 um, poorly connected and partly occluded by fibrous calcite cements or lime mud sediments.Intra-particle porosity may also develop at the contact of pseudo-cortex (fibrous cements) and the peloid grains (Fig. 6g,h).
Intercrystalline porosity is very common among calcite cements (Fig. 6e,f) and crystalline calcite (Fig. 6g).These pores are created by dissolution of carbonate crystals, grains and cements.Intercrystalline micro-porosity (Fig. 6f) forms in the limestone matrix of finely crystalline calcite with dispersive halite crystals as micro-pore infills.Pin-point pores 1-3 um in size are variably developed on the surface of calcite crystals (Fig. 6g).
Moldic porosity formed by selective dissolution of crystals within peloid grains.The void spaces left by dissolution of peloids are partly plugged with internal lime mud sediments (Fig. 6f).Fenestral vugs (0.2-0.5 mm in size) are commonly developed above algal laminations (Fig. 6e), partly filled by calcite cements and internal lime mud sediments.Fenestral vugs may have resulted from degassing associated with degradation of organic matter.Karstic dissolution has enlarged and created horizontal connections of fenestral vugs that are well documented on med-CT images (Fig. 7).Coarsely crystalline euhedral calcite cement fills vugs and reshapes the original pores.

Guelph Formation
In samples of Guelph Formation paleokarst rubble, the average matrix porosity (connected porosity) is ~5% (Fig. 8c), whereas the total porosity (assuming uniform mineralogy throughout the core) is around 12% (Fig. 8c).Large vuggy macropores locally increase porosity to 14%.The largest macropores appear to be connected through a network of fractures and small pores and micropores unresolved at the med-CT scale (see for example Fig. 8h).The CT scan cannot totally resolve small pore spaces (<0.03 mm 3 ), but the mixel values give a representation of partially filled/plugged pore spaces (e.g.grey gradient zones in Fig. 9), which are consistent with petrographic observations in thin sections (Fig. 6f).
Overall, the connected porosity in Guelph dolostone varies from 8.47 to 11.62%, with an average value of 9.23% and a standard deviation of 1.08%.By comparison, conventional core analysis of two samples (DGR-8 and DGR-4) reported average connected porosity of 14.5 and 12.6%, respectively.
The Guelph dolostones have highly porous intervals along depth, with large vuggy macropores responsible for rapid variations (Fig 8b,d).CT images and 3D renderings indicate the connection between macropores is through an irregular network of fractures and karst conduits.There is also an important contribution of small pores that cannot befully resolved at the med-CT scale but are visible as gray gradient on both dry state and two-state images.

A-1 carbonate unit
A-1 Carbonate Unit samples have an average matrix porosity along depth of up to 4%, but locally porosity per slice could be up to 18%.Connected porosity is formed predominantly of rounded to elliptic macropores (Fig. 8f).The largest macropores are not homogeneously distributed within the sample but rather concentrated in centimetric thick layers with connected porosity per slice up to 15%.
For the A-1 Carbonate samples, the average connected porosity is considerably more variable than those of the Guelph samples, ranging from 1.48 to 14.46% with an average value of 8.06% and a standard deviation of 4.82%.Conventional core    (c and g) Porosity profile calculated using the dry state (solid blue line) data where solid phase density is considered to be calcite density throughout the sample and two-state saturation protocol (dashed orange line).(d and h) CT-image after 'two states' saturation.The areas in white correspond to areas where water replaced air after saturation.The greyscale represents the amount of water filling the porous rock material.In other words, the greyscale indicates the porosity value in every voxel, from 0% (black) to 100% (white).A subvertical karst conduit is clearly visible in the Guelph Formation (d), and a network of fractures with elevated porosity in the A-1 Carbonate (h) is interpreted as healed desiccation cracks small pores that cannot be fully resolved by the med-CT scale imaging technique.

Macroporosity estimation following segmentation
Three-dimensional renderings of segmented pore networks provide an interactive tool to visualize pore distribution and pore connectivity, and can also inform on the concentration of pores along depth.In Guelph Formation dolostones samples, the number of pores can vary significantly along depth, and their volume as well, from a few cubic millimeters up to several hundreds of cubic millimeters.Very small pores (<10 mm 3 ) are very abundant (Fig. 10), but very large pores (500 mm 3 and higher) are the major contributor to porosity, representing in some instances 96% of the pore volume.For A-1 Carbonate samples, small pores are very abundant, but large pores, 80 mm 3 and higher, often represent 50% of the total pore volume and are often concentrated in specific layers (e.g.Figs.8f,i and 9e-h).To complete the analysis, the threshold sensitivity was checked on six samples (Fig. 12).The macroporosity varies linearly with the value of the threshold, nevertheless the samples fell into two groups.For five samples, a 1% change in threshold value changed the estimated microporosity by ~0.4% and for one sample it was ~0.1% (Fig. 12).These two groups have different macroporosity percentages and distinctly different initial threshold values selected.For samples with greater macroporosity (3-10%), the initial threshold value was 10%.For the other sample (DGR1_22.3),the Fig. 10 Statistical analyses for two illustrative segmentation analyses.For the Guelph Formation a threshold of 13% and an opening of 2 pixels were selected (a-d), whereas for the A-1 Carbonate, a threshold of 10% and an opening of 2 pixels was used (e-h).a and e Frequency plots of porosity distribution with core depth (b and f) -Pore size distribution with depth, with pore volume in mm 3 .c and g 3D volume representation of the pore network with the color bar on the right-hand side representing the amount of water that was replacing air space.This information can be seen as the 'porosity intensity'.The interconnected network of porous fractures in the bottom half of the core (g) are interpreted to be healed desiccation cracks.On the left-hand side, 3D renderings of a segmented pore network where the different colors represent one connected pore space in 3D.(d and h): histogram of the pore contribution to the total volume of rock.Note: pores bigger than 500 mm 3 have been set to 500 mm 3 for graphic purposes only.The actual volumes were used for porosity calculations Fig. 11 Summary chart comparing averaged core porosity for samples of Guelph (a) and A-1 Carbonate (b) using three different methods.For the Guelph Formation, the averaged core value obtained for DGR8_22.6 is negative, due to a spatial registration issue, prior to the 3D volumes subtraction.Only four full diameter core samples were analyzed using a conventional gas porosimeter Fig. 12 Sensitivity analysis of the threshold selection on the macropores value (in %).For a selection of six samples, three threshold values were selected based on visual analyses on 2D sections (Fig. 9).This sample selection illustrates that the samples with the coarsest pore network (e.g.large fracture; DGR3-21.6;DGR1-22.2) are the most sensitive to the choice of threshold value porosity was <3% and initial threshold value was 5%.This highlights the variability amongst samples, with the greatest sensitivity for samples with higher porosity.The choice of threshold values is a delicate one that was done on a sampleby-sample basis based on the visual analysis of the macropores and knowledge of the samples.

Bulk porosity method comparison
The whole core porosity was obtained by three different methods-gas porosimetry, dry-state porosity and twostate porosity (Fig. 11).Conventional gas porosimetry is considered to be the industry standard for reservoir properties studies (Anovitz et al. 2015;Njiekak et al. 2018), but is inherently different than the proposed two-state porosity estimation since the fluids involved are different (gas versus water).Additionally, the med-CT data was segmented to provide a macroporosity analysis and porosity values are therefore expected to be higher than porosity obtained using conventional gas porosimetry.
The whole core CT analysis provides an opportunity to assess porosity along depth over individual lengths of cm and in composite sections (e.g.Figs.9).The connected porosity values obtained using the CT two-state datasets are lower than the ones obtained using gas porosimetry and this is expected since the CT dataset, with its coarse resolution does not fully account for the contribution of micropores.Therefore, the two-state porosity values can be interpreted as the macropores contribution to total porosity.The connected porosity values obtained using the two-state datasets show how dramatically the porosity value can change along the depth.Porosity profiles along sample length give an average porosity value for each horizontal slice every 100 μ (Fig. 9d,i).

Discussion
The association of aquifers with shallow karst in carbonate bedrock in sedimentary basins is well documented (e.g.Bakalowicz 2005;Ford and Williams 2007;Chen et al. 2017).Paleokarst and paleokarst aquifers are not as well documented but any sedimentary basin with significant accumulations of carbonate rocks is likely to have paleokarst in the intermediate to deep subsurface (Simms 2014).Although open conduits are not preserved, paleokarst intervals have enhanced porosity and permeability conducive to groundwater flow.Paleokarst horizons are the principal geological control on confined carbonate bedrock aquifers in southern Ontario (Carter et al. 2021a).This is the first study that attempts to characterize the porosity systems in two of these paleokarst aquifers and interpret the physical processes that contributed to formation of connected pore systems in these rocks.The study methods used here may have wider applications.
The interpinnacle Guelph Formation was first recognized as a regional paleokarst horizon by Smith (1990) at a site over 160 km to the south, where it was described as a paleosol breccia.The nearly complete destruction of primary depositional features and development of vertical karst conduits, now filled by rubble, and regional distribution beneath most of southern Ontario, are indicative of mature karst developed over an extended period of subaerial exposure.This paleokarst forms the regional Guelph Aquifer (Carter et al. 2021a) and is sealed by the overlying Salina Aquitard.Three-dimensional hydrostratigraphic modelling by Carter et al. (2022) illustrates the regional extent of the aquifer and hydrochemical zonation of the confined groundwater from shallow freshwater, to intermediate saline to brackish sulphur water, to deep brine within the aquifer.Interpolated static level surfaces show up-dip flow directions in the deep brine regime.Flow directions are down-dip in the shallow to intermediate regimes with isotopic evidence of penetration of glacial meltwater into the intermediate regime when Ontario was covered by continental ice sheets during the Pleistocene (Carter et al. 2021a).
The A-1 Carbonate paleokarst is interpreted to be an epikarst with only a subregional extent and preservation of primary sedimentary structures such as peloids and algal mat laminations.Sediment-filled desiccation cracks are clearly visible in CT imaging.These cracks are interpreted to have allowed penetration of meteoric water into the laminated microbialites in the upper A-1 Carbonate when they were exposed on the sabkha surface, with lateral dissolution controlled by algal mat laminations.The slightly higher porosity of the cracks is hypothesized to be the result of crack infill by wind-blown silt as observed in modern supratidal sabkha environments (e.g.Shinn 1983; Kendall 2010).Petrographic characteristics are not known as none of the cracks were transected by a thin section in this study.
The paleokarst aquifers have horizontal hydraulic conductivities that are several orders of magnitude higher than their carbonate rock precursors (Carter et al. 2021a).Hydraulic conductivity measurements were obtained by Intera (2011) from the same boreholes at the Bruce site.There is no nonkarstic Guelph Formation for comparison, but the lithologically similar underlying Goat Island Formation has hydraulic conductivity of 9 × 10 -12 m/s vs 3 × 10 -8 m/s for the Guelph paleokarst.The nonkarstic lower A-1 Carbonate at the Bruce nuclear site has hydraulic conductivity of only 9 × 10 -12 m/s vs 2 × 10 -7 m/s for the A-1 Carbonate paleokarst.
Hydraulic conductivity of the two paleokarst aquifers is also several orders of magnitude less than a karst limestone in a shallow groundwater environment (Freeze and Cherry 1979).This is interpreted to be the result of mechanical compaction and recrystallization during burial diagenesis in the deep subsurface over very long periods of time, dolomitization of the Guelph Formation, and partial infill of karstic pores by carbonate cements and sediments.Former karstic conduits in the Guelph Formation are infilled with rubble, late dolomite cements and dolomudstone, which has created a complex and heterogeneous pore network where groundwater storage and movement are controlled by interconnected pores and not by karstic conduits.There are abundant micropores in both aquifers but large macropores contribute to the largest proportion of the pore volume in both.
The complexity of porosity types and distribution in carbonates is well known and has been the subject of intensive studies and discussions (e.g.Choquette and Pray 1970;Wilson 1980;Ahr et al. 2005;Lønøy 2006).This study corroborates the advantages of combining methods and scale of observations when assessing porosity and particularly in complex systems such as paleokarst (Anovitz et al. 2015;Niiekak et al. 2018).Large-scale features such as vugs, fractures and primary bedding were best observed visually by core logging.The CT scanning technology provided 3D visualization of pores and pore connections at macroscopic to microscopic scales, and their heterogeneity and complexity, and enabled quantitative assessments of the heterogeneity of porosity distribution.Healed desiccation cracks in the A-1 Carbonate which were readily observed with the CT images were unrecognized in macroscopic examination and were not intersected by any thin sections.Conversely, dispersed spherical and ovoid macropores in the A-1 Carbonate were identified as moldic dissolution of peloids using optical thin sections petrography, which could not be determined using only CT imaging.SEM provided images of pores and cements and elemental analyses at the nanometric scale.The biggest disadvantage/difficulty in applying the combination of multiscale methods used here may be the scarcity of fulldiameter drill core in aquifers.These are readily available in studies of oil and gas reservoirs but are rarely acquired in aquifers, particularly in deep paleokarst where groundwater salinity is greatly elevated.
Compared to gas porosimetry, CT scanning provides a continuous porosity profile, is nondestructive, and produces similar quantitative results in terms of porosity evaluated on core plugs (Fig. 11).By providing a continuous record of porosity along several tens of centimeters, the med-CT technique solves the issue of core plug representativeness in reservoir analyses.This is crucial in karst studies when plug-derived porosities would be strongly dependent on the location of core plugs because of the presence of large vugs (e.g.Figs.7 and 9c).As for its limitations, med-CT analysis and data processing is still time consuming (and thus expensive), particularly when computing power is limited, and is therefore not yet widely available at a commercial scale.There are also limitations in spatial resolution and the need for human-based expert decisions in choosing porosity masks to produce statistical analyses.That being said, in this case, the goal was to characterize pores geometries and their distribution in 3D.Regardless of the choice of threshold, the pore shape does not vary significantly and neither does the pore distribution.This is common practice in the digital rock characterization field.Compared to the increasing number of pore-scale studies done in micro-CT (e.g.Blunt et al. 2013;Noiriel et al. 2021), the med-CT requires less computing power and has the potential to bridge the scale gap between indirect borehole measurements and direct core-plug analyses.It has been used in the oil and gas sector for decades (Wellington and Vinegar 1987;Cromwell et al. 1984).

Summary and conclusions
This study is a rare examination and characterization of porosity development in two saline paleokarst groundwater aquifers from the pluri-centimeter (core logging, med-CT) to the microscopic (optical microscopy, SEM) scale utilizing medical CT scan technology and a two-stage saturation protocol.Med-CT scanning for pore analysis is more commonly used in petroleum reservoir studies (e.g.Withjack et al. 2003) as full-diameter drill core is usually readily available for reservoir studies due to the economic value of the hydrocarbons, but is rarely acquired in saline aquifers.
Rubble-filled subvertical karstic conduits in the Guelph Aquifer demonstrate the former presence of a groundwater conduit flow system.Mechanical compaction and diagenetic recrystallization in a deep burial environment has plugged the conduit network.This has created a complex and heterogeneous pore network where groundwater storage and movement is controlled by interconnected pores and not by primary karstic conduit geometry.The Guelph paleokarst represents a longer period of subaerial exposure than the A-1 Carbonate as evidenced by its greater thickness and geographic extent, the nearly complete destruction of primary sedimentary fabrics, and the presence of rubble-filled subvertical karst conduits.Large, irregular, vuggy macropores formed by karstic dissolution form the greatest percentage of the porosity in the Guelph Formation, with dolomite infilling of pores and common occurrence of crude oil lining pore spaces.Med-CT images and 3D renderings indicate the connection between macropores is made through an irregular network of fractures and subvertical karst rubble conduits with poor horizontal connections.There is also an important contribution of very abundant small pores (<10 mm 3 ).
In the A-1 Carbonate preservation of primary features, including algal laminations, peloids and desiccation cracks, has been beneficial for horizontal pore connections.Rounded to elliptic macropores are concentrated in layers associated with algal laminations, with good horizontal connectivity and porosity of up to 15%.Despite the greater maturity of the Guelph paleokarst the A-1 Carbonate has better hydraulic conductivity of 2 × 10 -7 m/s vs 3 × 10 -8 m/s for the Guelph, which is attributed to the horizontal interconnections of macropores along algal laminations.
Medical CT scans provided 3D visualizations of pore networks within the core samples that greatly enhanced the understanding of pore distribution and pore connections in the two paleokarst aquifers.Porosity calculations using med-CT scans have similar results to conventional full-diameter gas porosimetry.The technology has excellent capability for nondestructive characterization and analysis of porosity, including: • Ability to visualize the heterogeneity of porosity distribution within the entire core volume in three dimensions • Quantification of pore size range, distribution of pore sizes and pore volume attributable to each pore size, and relate these parameters to a location within the core • Visualization and quantification of connected pore space, a necessary requirement of permeability • Total porosity estimates similar to conventional gas porosimetry Saline paleokarst aquifers are probably widespread in other sedimentary basins but remain poorly documented.The combination of methods and approach developed here have the potential to be applied in other sedimentary basins and pore systems in paleokarst and also other aquifer types, where understanding of the multiscale complex porosities is difficult to capture.The present study also addresses the general lack of real-world examples for construction of physical models for paleokarst aquifers (e.g.Mohammedi et al. 2021) Acknowledgements An internal review was completed at the GSC by Alexandre Desbarats.Assistance of INRS CT-scan laboratory staff, including Mathieu Des Roches, was critical to the completion of the study.The NWMO supported access to the DGR drill core at the Bruce nuclear site.The authors would like to thank journal reviewers Mohamed Soufiane Jouini, two anonymous reviewers, and the journal editor Maria-Theresia Schafmeister whose comments have significantly improved this manuscript.
Funding Open Access funding provided by Natural Resources Canada.Funding for this study was provided by the Natural Resources Canada (NRCan), Geological Survey of Canada (GSC) and NWMO.This work is a contribution of the GSC Groundwater Geoscience Program Archetypal Aquifer project.This is NRCan contribution 20220439.Conflict of interest On behalf of all authors, the corresponding author states that there is no conflict of interest.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made.The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material.If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.To view a copy of this licence, visit http:// creat iveco mmons.org/ licen ses/ by/4.0/.

Fig. 1 a
Fig. 1 a Location map showing study area in red.b Regional bedrock geology of southern Ontario (Canada) showing northwesterly trending subcrop belts of Paleozoic bedrock formations, subcrop boundaries of the Guelph Formation and A-1 Carbonate Unit, location of Bruce nuclear site, and county boundaries.Bedrock formations dip

Fig. 2 a
Fig. 2 a Facies belts of the Guelph Formation in the subsurface of southern Ontario showing extent of the Guelph interpinnacle paleokarst and the Bruce study site location.County boundaries shown for geographic reference.Modified from Carter et al. 2021a.b Inferred

Fig. 3 Fig. 4
Fig. 3 Location of sampled deep boreholes at the Bruce nuclear site

Fig. 7
Fig. 7 Med-CT images in a dry state illustrating macroscopic pore spaces.a Guelph sample DGR4-21.3 at 390.61 m displays a complex network of dissolution-enhanced vugs forming an interconnected macropore network.b Guelph) sample DGR8-21.8 at 382.13 m displays solution-enlarged fractures, and cavernous vugs.c A-1 Carbonate sample DGR1-22.1 at 325.95 m displays rounded to elliptic

Fig. 8
Fig. 8 Compilation of med-CT results for (a-d) Guelph dolostone sample DGR 3-21.4 (at 390.61 m) of and (e-h) A-1 Carbonate sample DGR1-22.3 (at 328.2 m). a and e Density profiles, in Hounsfield unit (HU), measured along the length of sample.b and f Coronal CT-scan image, scanned in a dry state, where black areas correspond to open spaces whereas white areas correspond to the densest materials.(c and g) Porosity profile calculated using the dry state (solid blue line) data where solid phase density is considered to be calcite density throughout the sample and two-state saturation protocol (dashed orange line).(d and h) CT-image after 'two states' saturation.The areas in white correspond to areas where water replaced air after saturation.The greyscale represents the amount of water filling the porous rock material.In other words, the greyscale indicates the porosity value in every voxel, from 0% (black) to 100% (white).A subvertical karst conduit is clearly visible in the Guelph Formation (d), and a network of fractures with elevated porosity in the A-1 Carbonate (h) is interpreted as healed desiccation cracks

Fig. 9
Fig.9Threshold selection for sample selected to be representative of Guelph dolostone (a-d), and sample DGR1-22.3 (at TVD 328.2 m) for the A-1 Carbonate (e-h).a Coronal section through the connected porosity volume of Guelph sample, i.e. 3D volume obtained after 'two states' saturation.b-d Three different segmentation scenarios based on threshold and opening selections.Subvertical karst conduit is clearly imaged.e Coronal section through the connected porosity volume for A-1 carbonate sample.f-h Three different threshold segmentation based on different threshold and opening selections.Scenarios g and h seem to be best suited to illustrate the concentration of pores in centimetric layers

Table 1
Medical CT-scanner parameter values for acquisition and reconstruction stages.mAs milliampere-second

Table 2
Comparison of multi-scale investigation methods used in this study

Table 3
Summary of porosity types in the A-1 Carbonate and Guelph paleokarst aquifers