Organic-rich shale samples from a lacustrine sedimentary sequence of the Newark Basin (New Jersey, USA) are investigated by combining Broad Ion Beam polishing with Scanning Electron Microscopy (BIB-SEM). We model permeability from this 2D data and compare our results with measured petrophysical properties. Three samples with total organic carbon (TOC) contents ranging from 0.7% to 2.9% and permeabilities ranging from 4 to 160 nD are selected. Pore space is imaged at high resolution (at 20,000x magnification) and segmented from representative BIB-SEM maps. Modeled permeabilities, derived using the capillary tube model (CTM) on segmented pores, range from 2.3 nD to 310 nD and are relatively close to measured intrinsic permeabilities. SEM-visible porosities range from 0.1% to 1.8% increasing with TOC, in agreement with our measurements. The CTM predicts permeability correctly within one order of magnitude. The results of this work demonstrate the potential of 2D BIB-SEM for calculating transport properties of heterogeneous shales.

1. Introduction

The pore network is the main control on transport processes in low-porous media, such as gas shales [1, 2]. Thus, it is of interest to develop improved methods for understanding controls and to upscale pore structural information from nanometer-scale to larger rock volumes [3]. Argon Broad Ion Beam (BIB) and Focused Ion Beam (FIB) milling are techniques that offer new insight in microstructures through high-resolution imaging. While BIB milling is commonly used to analyze mm-sized 2D sections, FIB is limited to μm-sized volumes [4, 5]. By using these imaging methods, it is possible to receive a solid idea of a shale’s internal (micro) structure and pore geometry (e.g., [611]). BIB milling in combination with SEM imaging delivers direct information about pore sizes, their morphology and spatial distribution, rock microstructure, and potential anisotropy as well as core damage due to sample handling (e.g., [3, 1225]). Additionally, in combination with energy dispersive X-ray spectroscopy (EDX), BIB-SEM enables the identification of individual grains [26, 27]. However, subjective selection of representative study areas on the cross-section as well as the resolution—a major part of the pore space can exist below resolution—restrict these techniques. Due to the small-scale nature of pores in shales other imaging techniques, for example, MicroCT, are not suited for pore imaging. So far SEM is the most practical tool to image these pores properly [28]. Nevertheless, it remains challenging to upscale pore system characteristics from the nanometer-scale of SEM investigations to the centimeter scale of petrophysical experiments to the meter scale of reservoirs.

As the properties of a rock’s pore system directly govern permeability, image-derived pore-scale models (also known under the term digital rock physics) are used to study porous media, especially for conventional rocks, such as sandstone or carbonate [29, 30]. However, building a pore network model for 3D flow simulations in heterogeneous fine-grained rocks is an unresolved issue. Some approaches extract 3D pore data directly from FIB-SEM [31] or use the Markov chain Monte Carlo (MCMC) reconstruction method on 2D images [32, 33]. Whether these complex pore network models are representative for the sample is another uncertainty [17, 31, 34]. Kelly et al. 2015 pointed out that FIB-SEM investigation volumes of shales are not suited to resolve local heterogeneities, and thus, larger areas of investigation are necessary. Moreover, extracting a percolating pore network in shales with FIB-SEM remains a challenge since most pore throats are below the pixel resolution [15, 20, 35, 36].

In contrast, this work aims to determine pore system properties and permeabilities by BIB-SEM imaging in 2D from three different shale samples. Permeability predictions are based on a simple model [37, 38] that relates pore microstructure with permeability. The BIB-SEM derived permeability results are then validated by experiments on the exact same sample plugs.

2. Materials and Methods

2.1. Sample Material

Three samples (Table 1) with a decreasing porosity, permeability, and TOC trend according to their labeling (NJ-001, NJ-019 and NJ-023) were selected from cores drilled during the Newark Basin Coring Project (NBCP) from 1990 to 1993. These laminated, organic-rich lacustrine shales from the Triassic Passaic and Lockatong formation (~214 to 222 Ma) are of good quality featuring no macroscopically visible fractures or secondary mineralization (e.g., [3942]). From each sample location, a cylindrical sample plug was drilled parallel to the bedding, which was used for porosity and gas permeability measurements as well as BIB-SEM imaging.

2.2. Bulk Measurements
2.2.1. XRD, TOC, and Analyses

Bulk mineralogical compositions were derived from X-ray diffraction (XRD) patterns of randomly oriented powder preparations. The measurement was done on a Bruker AXS D8 Advance diffractometer using a CuKα-radiation produced at 40 kV and 40 mA (analysis performed by Rietveld refinement). Details of the sample preparation and evaluation are described in [43].

Total organic carbon (TOC) and total inorganic carbon (TIC) were measured on powdered samples with a liquiTOC II analyzer. The instrument analyzes the released CO2 during one temperature ramp without previous acidification.

Vitrinite reflectance () measurements were performed on polished blocks under oil immersion ( = 1.518) using a Zeiss Axio Imager microscope. Details on sample preparation, the analytical procedure, and instrumentation are described in [44, 45].

2.2.2. He-Porosity and Gas Permeability

He-porosities were calculated by combining skeletal densities from helium expansion (pycnometry) and bulk densities from cylindrical sample plug dimensions ( = ). Details on the setup and measuring procedure were previously described in [46, 47].

Gas permeability coefficients were measured on the same cylindrical sample plugs drilled parallel to bedding with helium gas as permeate at 25°C (298 K) in dry condition. They were installed into triaxial flow cells and then loaded to a confining pressure level of 40 MPa. After installation, the system was flushed with helium and leak-tested. Permeability measurements were conducted at confining pressure levels of 40, 30, 20, and 10 MPa during unloading. At each confining pressure level, nonsteady state flow tests (pressure pulse decay) were then performed at various pore pressures from 0.5 to 5 MPa. From the pressure incline/decline versus time series, apparent permeability coefficients were calculated. Details on the experimental setup, pressure pulse decay tests, and gas permeability calculation procedure are described in [46, 48].

Due to gas slippage, measured apparent gas permeability coefficients are higher than “intrinsic” permeability coefficients. Gas slippage effects were corrected by applying the Klinkenberg-correction on the apparent permeability data of a given confining pressure level [46, 47, 49, 50]. The Klinkenberg-corrected permeability versus stress couples was then extrapolated by an exponential expression to “unstressed” conditions to be able to compare them to the results of “unstressed” BIB-SEM measurements [4951].

2.3. Microstructural Investigation
2.3.1. Sample Preparation

The end of all plugs used for the permeability measurement were cut off and used for BIB-SEM investigation. Subsamples were cut dry into rectangular blocks of 3 × 5 × 5 mm using a low-speed microdiamond saw. Subsampling was based on macroscopic investigations, that is, checking for milling locations perpendicular to the bedding plane featuring visibly different kinds of layers, for example, darker (more clayey) and brighter (siltier) layers. These locations were then BIB-polished by a JEOL SM-09010 polisher to produce planar, Gaussian-shaped cross-sections of approx. 2 mm2. The samples were subsequently coated to prevent charging of the sample during imaging. For details about the protocol see [21].

The end trim of the plugs used for the permeability measurement of samples NJ-001 and NJ-019 was additionally polished perpendicular to the bedding using sand paper (SiC) in preparation of cm2 scale energy dispersive X-ray spectroscopy (EDX) analysis.

2.3.2. BIB-SEM Analysis

The Field Emission SEM used for image acquisition in this work is a Zeiss SUPRA 55 equipped with a backscattered electron (BSE), secondary electron (SE2), and EDX detector.

Pores can be identified on BIB-polished surfaces using a SE2 detector at low acceleration voltages of 3 to 5 kV [21]. The BSE detector was used to record the density contrast at the same area on the cross-section applying acceleration voltages of 15 to 20 kV. Element maps, showing the intensity of an element, were produced simultaneously using the EDX detector. Subsequently, the recorded images were stitched together automatically forming large mosaics by using the software Oxford Instruments© AZtecEnergy. All BIB-SEM images were scanned perpendicular to the bedding.

Several image mosaics were scanned at different magnifications to gain microstructural information at multiple scales. For each sample, at least one BIB cross-section was prepared, and typical layers were selected based on their mineralogical composition and mapped at high resolution with the SE2 detector, and when relevant, also with the BSE and EDS detector: () plug analysis on complete plug surfaces at 125x magnification (BSE and EDX with a pixel resolution of 2.4 μm); () subsample overviews of 250x magnification for identifying cross-section locations (BSE and EDX with pixel sizes of 1.2 μm and 2.4 μm, resp.); () detailed recording of the BIB-SEM cross-sections at 2.500x magnification (BSE and EDX with pixel sizes of 120 nm and 240 nm, resp.); plus (4) high-resolution mosaic maps of 10,000x to 20,000x magnification for in-depth mineral phase and porosity identification (SE2, BSE, and EDX with pixel sizes of 15 nm, 30 nm, and 60 nm resp.). High-resolution mapping comprises three mosaic maps (SE2, BSE, and EDX) for sample NJ-001 and NJ-019, each, plus one map for sample NJ-023.

Assuming that the pore space is strongly linked to the mineral phases, the representative elementary area (REA) of mineralogy is an approximation for the representative area of porosity [19]. A box counting method (as described in detail in [21]) was employed to determine decreasing fluctuations in the mineralogical distribution (±5%) at which an area is considered to be representative.

2.3.3. Digital Image Analysis

The segmentation of mineral phases included thresholding and automatic image treatment based on BSE and EDX maps using different toolsets within ESRI© ArcGIS. The mineralogical composition was quantified by assessing the element intensities of the EDX maps, where high Si counts were interpreted as quartz, high K, and intermediate Si as K-feldspar, high Na as Na-feldspar (albite), and high Ca counts as carbonates (dolomite/ankerite). Pyrite and OM were segmented from BSE micrographs and the remaining area was considered as clay. Besides an identification of mineral phases, mineral grain sizes, shapes, and spatial distributions were qualitatively assessed.

Porosity was segmented from SE2 images utilizing a “seed and grow” algorithm [52]. Appropriate morphological indicators to quantify and distinguish distinctive characteristics between pores or cracks are elongation, derived through a pore’s axial ratio (see Appendix A, (A.1)), and roundness, obtained by the circularity of a single pore’s boundary (see Appendix A, (A.2)). For the flow model, segmented pores were converted into polygons and evaluated by using the capillary tube model (CTM).

2.3.4. Upscaling to Larger Rock Volumes (Plugs)

Macroscale EDX analysis of the surface of the plug’s end trim with 125x magnification allowed distinguishing layers with typical mineralogy that can be used for upscaling to a respective plug’s mineral content. To compare permeability predictions of a single layer to the measured permeabilities of respective plugs (both parallel to the bedding), typical layers of sample NJ-001 and NJ-019 were assessed via EDX intensity measurements. The total plug area was then separated into typical layers to obtain the area proportions of each representative layer.

To identify the properties of distinct layers, each layer is characterized by a BSE/EDX mosaic map (Figure 8). The quantitative properties of these maps, in particular the dominant mineralogical composition, are then used to distinguish three main types of layers within the plug. This procedure allows extrapolating representatively to the plug scale (cm). Finally, CTM estimated permeabilities of each single mosaic map are normalized to the 2D areal fraction of the plug’s respective layer surfaces to receive realistic (upscaled) permeabilities (Table 5).

2.4. Capillary Tube Modeling

We used the CTM similar to [51]. It assumes idealized capillary bundles [37, 38] and represents the permeability by a combination of Darcy’s law with the Hagen–Poiseuille equation, describing laminar flow through a cylindrical pipe (see Appendix B, (B.1)–(B.10)). With the tortuosity (ratio of the true length of a pore path compared to the straight-line distance between those points ) and the hydraulic pore radius (to consider the complexity of the pore shape, pore area divided by its perimeter) of a capillary system, the individual flow through a single pore is described by the permeability (m2):where stands for the porosity. However, since there is no uniform pore radius in these shales and BIB-SEM enables quantification of a range of pore sizes, the equation above is rewritten to derive permeability for a single pore (m2) by taking into account its effect on porosity (a pore’s area (m2) divided by the total mosaic area (m2)) and sum these up for the total permeability (m2):

The tortuosity characterizes flow pathways in shales and describes the connectivity of the pore network. If tortuosity equals 1, flow may be fully described by the Hagen–Poiseuille equation [53]. Values of 10 imply a traversed flow path of ten times the actual (theoretically shortest) path length. In accordance with the work of [51] a tortuosity of = 10 is applied in this study.

Besides, modeling of the cumulative permeability coefficients can be optimized by correcting the visible porosities for isolated IntraP-pores and induced microcracks.

3. Results

3.1. Bulk Properties
3.1.1. Mineralogy and Maturity

The samples vary in TOC, with the highest organic content in sample NJ-001 and lowest organic content in sample NJ-023 (Table 2). The maturity of all samples ranges from ~2.1 to 2.7% . Their mineral composition is different from organic-rich marine shales with respect to () absent or low quartz content from <1 to 13%; () significant higher proportion of feldspar(s) with ~36 to 47%; and () a similarly high amount of carbonates with ~20 to 38% throughout all samples (Table 3).

3.1.2. He-Porosity and Gas Permeability

Bulk densities vary between 2.54 g/cm3 (NJ-001) and 2.68 g/cm3 (NJ-023) and grain densities between 2.68 and 2.72 g/cm3. As a result, calculated He-porosity values of sample NJ-001 are highest with approx. 5.1% followed by NJ-019 with 4.1% and NJ-023 with 1.4% (Table 2).

Similarly, Klinkenberg-corrected permeabilities (extrapolated to zero stress) decrease from sample NJ-001 to NJ-023. Sample NJ-001 shows highest permeability with  nD (approx.  m2) followed by NJ-019 with  nD (ca.  m2) and sample NJ-023 shows lowest permeability with  nD (approx.  m2) (Table 2).

3.2. Qualitative Description of Mineralogy, Microstructure, and Pore Morphology

All samples are characterized by a very heterogeneous fabric with a variety of different clasts embedded in a fine-grained matrix framework. Grains occur in a wide range of sizes and mostly touch each other. Layering lies within submillimeter to centimeter scale (NJ-001 slightly more distinct than NJ-019 or NJ-023). In between clasts, dispersed OM, intercalated by clay minerals, is evident and qualitatively analyzed to be gradually less from sample NJ-001 to NJ-023 according to the chosen TOC sequence. The given samples feature three typical submillimeter layers with an abundance of large albite (Figure 1(a)) or carbonate (dolomite/ankerite) clasts (Figures 1(b) and 4). Additionally, distinct sulfur-rich layers featuring a variety of euhedral pyrites are visible in sample NJ-001 (Figure 1(c)). As clasts in all samples are mostly in contact to one another, the fabric is classified as a grain supported matrix. The matrix itself (gray) features dispersed OM (black) intercalated by other types of minerals (more or less elongated) possibly belonging to the phyllosilicates (e.g., mica).

The feldspar-rich laminae are responsible for the high porosities in the investigated samples (Figure 1 versus Figures 4 and 5), especially as NJ-001 features significantly more silty layers on the plug scale (Figure 8). Similar laminae, responsible for regulating flow, were also found by Lei et al. 2015 [54] in Chinese organic-rich lacustrine shales.

Mineralogy and Microstructure. Carbonate grains (blue, Figure 1) represent the biggest and most distinct clasts of typical rhombohedral habits, ranging in size from <1 μm up to approx. 50 μm in diameter in sample NJ-001. In sample NJ-019 and NJ-023 Ca-clasts (blue, Figures 4 and 5) are less common and typically not larger than 20 μm (commonly <5 μm) in size. Feldspars (pink, Figure 1) also feature some relatively large clasts of up to about 25 μm, but their size is mainly below 10 μm. In contrast to sample NJ-001, the samples NJ-019 and NJ-023 exhibit two different types of feldspar, sodium-containing feldspar, here albite (lower gray values on the BSE) as well as potassium-containing feldspar (higher gray values on the BSE). The morphology of albite clasts (pink, Figures 4 and 5) is relatively distinct and their size represents the biggest clasts within these samples with less than 5 μm up to approx. 100 μm in diameter. The morphology of K-feldspar clasts (green) is similar, though their sizes are of only a few μm up to about 25 μm in diameter. In sample NJ-001, quartz grains (green, Figure 1) are usually small (<10 μm), though they exist also in sizes of more than 20 μm while quartz is almost absent in sample NJ-019 and NJ-023. Pyritic grains (yellow, Figure 1) of sizes up to 10 μm in sample NJ-001 are common while pyrite is very rare in sample NJ-019 and NJ-023.

Comparing the qualitative results to typical features of marine shales, lacustrine shales are often differentiated by () thin recurring (seasonal) laminae [55, 56]; () siltier layers, particularly rich in albite and carbonate (dolomite/ankerite); () no/little presence of calcareous fossils [18, 21, 57] or framboidal pyrite [25, 58]; () overall relatively low visible porosity; and () relatively low visible OM porosity for these high maturities.

Pore Morphologies. Porosity in all samples is predominantly distributed between clasts and OM interfaces (see also [59]) and thus allocated to InterP porosity according to the pore space classification of Loucks et al. [60]. However, if certain types of grains contain pore space, its typical morphology and distribution can be assessed as well (Figure 2).(i)Few albite grains embody large pores (IntraP) of up to ca. 10 μm (in direction of the longest axis), with angular morphologies, typically occurring at the clast’s rims progressing into the center (see also [61, 62]). Completely isolated and centered pores are considered as fluid inclusions (see also [63]) and excluded in the CTM as described in chapter 0 (cf. Figures 2(a) and 2(b)).(ii)Some carbonates feature cluttered pores (InterP) along their edges (cf. Figure 2(d)), plus very few grains show large pores (IntraP) of up to ca. 5 μm (cf. Figure 2(h)).(iii)Quartz and pyrite show almost no associated porosity (cf. Figures 2(b), 2(f), 2(k), and 2(l)).(iv)Relatively small pores are located within the dispersed occurring OM (cf. Figures 2(j) and 2(l)).(v)All other types of intergranular porosity (InterP) can be considered to be either very small pores or very narrow elongated slits along slightly dilatant grain boundaries interfaces (cf. Figures 2(e) and 2(f)).

The visible porosity of sample NJ-019 and NJ-023 is significantly less but occurs similar to the cavities described in sample NJ-001 (Figure 6). Similarly, pores are found predominantly in between the interfaces of OM and clasts (cf. Figures 6(b), 6(e), 6(f), 6(g), and 6(h)) and some large pores are situated within feldspar or carbonate grains (cf. Figures 6(a) and 6(d)). In contrast, clear clay platelets seem to feature elongated pores (cf. Figures 6(i) and 6(j)).

High magnified images revealed porosity (IntraP) that is omnipresent within the dispersed OM in sample NJ-001 (Figure 3(a)), but not within the dark, comparatively large individual particles of the OM (Figure 4(b)). Almost no pores are associated with the OM in sample NJ-019 and NJ-023 (Figure 6) which is affirmed by additional SE2 scans under very high magnifications (Figure 3(b)). Similar to the Triassic black shales investigated by Loucks et al. 2017 [64], most of the OM pores are located in migrated solid bitumen.

This indicates the presence of at least two different OM types, interpreted as(i)dispersed bitumen that migrated into cavities between grains;(ii)primary terrestrial vitrinite and inertinite particles.

Furthermore, the investigated lacustrine shale samples on average exhibit less visible porosity (~0.8%) than their marine organic-rich counterparts at similar maturity. In contrast to other studies on overmature organic-rich marine shales (e.g., [15, 22, 24, 26]), pores within the given lacustrine samples are rarely associated with the OM, even though the dispersed OM and primary terrestrial vitrinite or inertinite particles could be identified as OM Types A and C found in marine shales [22]. The dispersed OM identified in all samples only exhibits pronounced OM porosity in the high-resolution images of NJ-001 (Figures 2(j), 2(l), and 3(a)). Besides, the number of OM pores does not increase with further maturity as stated in numerous studies of lacustrine shales [65, 66], mostly since samples NJ-019 and NJ-023 contain less OM and little to none OM porosity. Considering the high vitrinite reflectance values () and the lack of significant visible large OM pores, the hydrocarbon potential can be classified as poor. The relatively small amount of visible OM porosity could also be related to the fact that lacustrine shales are dominated by kerogen of terrestrial source generally hosting less OM porosity [20, 67, 68].

3.3. Quantitative Description of Mineralogy and Pore Space
3.3.1. Mineralogical Composition

The mineralogical proportions of the investigated samples vary when comparing the microstructure of sample NJ-001 to NJ-019 and NJ-023, respectively (Figure 7). Qualitative trends described before, such as distinct layering in sample NJ-001, are recognizable quantitatively as well.

EDX results show that the mineral content varies most with up to approx. 27%, 16%, and 8% for carbonate, feldspar (albite), and pyrite, respectively, in the maps of sample NJ-001. Variability of the mineralogical proportions between maps of sample NJ-019 is up to approx. 36%, 28%, and 21% for clays, quartz/silicates (correctly classified as K-feldspar by XRD analysis), and albite, respectively (Figure 9). Bulk XRD results only exhibit profound variations in comparison to the EDX results of the mineralogy of sample NJ-023 with variations of clays (ca. 30%) and carbonates (ca. 33%). All other distinguished mineral phases in every sample vary by less than 7% compared to the EDX findings.

3.3.2. Porosity and Pore Morphology

All pores below 18 pixels (equal to 72 and 144 nm in diameter at 20,000x or 10,000x magnification, resp.) in size must be considered with caution since pores of these sizes are below the practical pore resolution (PPR). Microcracks (narrow elongated pores) that may originate from drilling, core recovery operations, drying, or sample preparation are identified and removed.

Total visible porosities between all mosaic maps range from approx. 2.0% (NJ-001) to almost no visible porosity (NJ-023), while considering only pores above the PPR lowers these values slightly (up to ca. 1.7%) (Table 4). Visible porosities differ between 1.4% in the carbonate-dominating layer (Map #2), 1.8% in feldspar/carbonate mixed layer (Map #3), and 2.0% in the feldspar-rich layer (Map #1) for sample NJ-001 (only Map #1 of sample NJ-001 was imaged with a pixel size of 30 nm instead of 15 nm). The highest visible porosities in sample NJ-019 (0.9%) are also located in the feldspar-rich layers (Map #2), while the clay and silicate-rich layers (Map #1 and Map #3) feature both about 0.5% (Figure 9). In general, the pore orientation follows the bedding.

Initial classifications of pore types according to Desbois et al. 2009 [4] and Heath et al. 2011 [69] distinguished pores not only based on the location of their occurrence (as proposed by Loucks et al. [60]) but also in consideration of their actual shapes. Similarities are observed in the given samples as well, featuring several different types of pores as previously described.

3.3.3. Proportion of Organic Matter Porosity

Very high magnification images (40,000 to 80,000x) revealed omnipresent porosity within the OM in sample NJ-001, while samples NJ-019 and NJ-023 show no OM porosity (Figure 3). Based on simple thresholding the porosity of the organic phase in sample NJ-001 was estimated to be ca. 6 ± 2%. Given the high OM proportion of about 8.4% in this sample, the OM porosity may be a significant fraction of the total porosity.

3.3.4. Pore Size Distributions

Pore frequency histograms with power law based bin sizes are given in Figure 10. The normalized distributions suggest a power law distribution over about four orders of magnitude, except a few outliers, for instance, some single large pores, which do not line up with the log-log best fit (Figure 11). The power law exponent varies between approx. 1.9 and 2.4 with and 1.8 and 2.3 without cracks. In sample NJ-001, small pores seem to contribute more to the total porosity due to their higher power law exponents () compared to NJ-019 ().

3.4. Permeability Predictions

Cumulated permeability graphs of each individual mosaic are presented in Figure 12. Modeling of the cumulative permeability coefficients was done after correcting the visible porosities for(1)isolated, relatively large IntraP-pores within albite (Figure 13);(2)pores below the PPR;(3)visible microcracks (all pores with an AR threshold value of ≤0.2 and circularity of ≤0.3) which are assumed to be artefacts from sample handling.

Through the investigation of characteristic plug layering, three representative layers of each formation were identified (Maps #1 to #3 of samples NJ-001 and NJ-019, resp.). For upscaling, the predicted permeability values of each mosaic are normalized (by the weighted average) to the total area proportions of each plug. Hence, the weighted average permeability of the three samples changes from 665.8, 323.1. and 2.3 nD to 310.1, 43.1, and 2.3 nD, respectively, after taking the mineralogical composition of the plugs NJ-001, NJ-019, and NJ-023 into account (Table 5).

The sum of modeled permeability coefficients of microcracks ranges from < 1 nD (NJ-001 Map #3) up to about 12 nD (NJ-001 Map #1) within all maps, depending on the frequency of these cracks. Pores below the PPR do not add up to relevant calculated values (predicted values of ca. < 0.1 to 2 nD between all maps).

4. Discussion

4.1. Comparison of Bulk Properties and BIB-SEM Results

The mineralogical compositions of EDX and XRD are similar for samples NJ-001 and NJ-019 (Figure 7), even though XRD results are much more accurate than measuring the bulk composition of an area via EDX. Accordingly, significant differences in carbonate and clay content of EDX versus XRD findings occur in sample NJ-023. The difference is related to sampling heterogeneity between the plugs and the cut-offs used for XRD due to rapid changes in the mineralogical composition within a few centimeters.

Except for the fraction of OM, no other mineral phase shows a clear correlation between its fraction and the porosity trend (Figure 14), implying that the OM is the controlling factor of total and visible porosity in the investigated samples, particularly in NJ.001. Although the recognition of OM porosity in sample NJ-001 lies mostly below the PPR with the given pixel resolution, a proportion of the OM porosity is visible at 20,000x magnification (Figures 2(j) and 2(l)). Also, a considerably higher power law exponent in NJ-001 compared to the other samples corresponds to the porous nature of the OM (Figure 11 and Table 4), as exponents above 2 are an indicator of a dominating contribution of small pores to the total porosity [18, 67]. On the contrary, the lower total porosity and permeability of sample NJ-019 and NJ-023 feature lower values. Additionally, differences in porosity along the investigated samples are most likely linked to the occurrence of silty layers. The biggest pores are located around feldspars and carbonate grains (as shown on the local porosity maps in Figures 1 and 4) implying that porosity and permeability are mostly bound to these two mineral phases.

Moreover, solid trends between mineralogy, microstructure, and pore properties throughout all investigated scales imply representative sampling. The bulk measurements exhibit clear correlations between petrophysical properties such as the trend of decreasing He-porosities as well as permeability coefficients along with decreasing TOC content (Table 2). Visible porosities (Table 4) as well as experimentally measured porosities (BIB-SEM versus MIP versus He-pycnometry) follow this trend as well.

4.2. Validating BIB-SEM Derived Permeability Predictions

Unlike sandstones, shales do not feature clear poro-perm relationships [51, 70]. As indicated by Busch and Amann-Hildenbrand [71], caution is advised when using permeability prediction models for shales, especially when the predictions are of empirical origin. In this study, we applied the theoretical CTM based on real pore geometries to predict permeability with an assumption on the tortuosity of the pore network.

We used the pore system characteristics of typical layers that were resolved by BIB-SEM for upscaling to plug scale (Figure 8) to allow a precise and meaningful comparison with measured permeabilities. We assume that the microstructural investigations are representative for certain dominant layers. Since NJ-023 features no clear layering and a very low porosity, we investigated macroscale layering perpendicular to the bedding only for sample NJ-001 and NJ-019 by EDX mapping.

For the permeability calculation, several adjustments on the pore data were conducted as they were expected to influence the modeled permeability coefficient. Very narrow, elongated pore shapes, identified as microcracks, do not influence the predicted permeabilities significantly (the predicted values of all maps are < 12 nD) (Figure 12). For the prediction of permeability using the CTM, pores below the PPR were excluded, even though they do not add up to relevant calculated values (the predicted values of all maps are < 2 nD). This is reasonable as their individual permeability coefficients are up to five orders of magnitude below the total permeabilities of the investigated micrographs, although pores below the PPR make up to ~39% of total pore frequency in some maps (Table 4). Considering the hydraulic radius rather than the geometric radius explains the minor difference between the permeability prediction of narrow elongated pores (i.e., microcracks) and tiny round pores (i.e., pores below the PPR). Additionally, in three of the total seven maps with considerable amounts of albite (NJ-001 Map #1, Map #2, and NJ-019 Map #2), large IntraP-pores within these albite clasts were excluded from the permeability calculations as they are assumed to be either completely isolated from the pore network (fluid inclusions) or “dead-end” pores. However, most pores typically occur at the rims of albite progressing towards the clast’s center, that is why they were not excluded in the permeability calculations.

Albite and its distinct porosity seem to control permeability since the albite rich maps show the highest calculated permeability (Table 5), even though isolated pores were excluded from the calculations. In general, porosity alterations are common in feldspars (e.g., [72, 73]). Carbonate clasts exhibit similar but slightly less porosity compared to albite, featuring the most pore space around grain boundaries (Figure 2(d)). However, this accounts only for carbonates in sample NJ-001 as their fabric is denser compared to the other samples (NJ-019 and NJ-023). The quantitative permeability predictions correspond well to these qualitative observations with the highest predicted permeabilities within the Na-rich maps (ca. 801 and 1,208 nD) and generally high permeability coefficients of the Ca-rich map in sample NJ-001 (ca. 538 nD) (Table 5). This suggests that the coarse-grained material is most important in controlling fluid flow due to less compaction resulting in larger pores. Furthermore, the porous nature of OM in NJ-001 suggests that the pores in the OM also add to transport pathways and add up to the calculated permeabilities. However, their small pore sizes probably result in small absolute permeabilities.

Limitations of stereological modeling apply to the CTM calculations as well. Obtaining adequate tortuosity values, especially from 2D data, is problematic. A 3D investigation could deliver different tortuosity factors. By choosing a tortuosity of = 10 a good fit between intrinsic and calculated permeability could be established. Values of around 5 to 10 are commonly used for permeability predictions in very tight rock samples that are strongly lithified and feature very low-porous rock fabrics (e.g., [71, 7477]), though one should be careful as different notations of the CTM exist in the literature (regarding versus ).

For validating the BIB-SEM based permeabilities, the results are compared to the experimentally derived permeabilities measured on the exact same sample plugs. Assessing the mineralogy and microstructure is key to allow upscaling to greater sample volumes for meaningful comparison. This was established by normalizing the results of the CTM to the plug scale in this study. The final results of the corrected and upscaled BIB-SEM permeabilities are 310 nD, 43 nD, and 2 nD for sample NJ-001, NJ-019, and NJ-023, respectively. These values are at maximum twice or half, respectively, when compared to the experimentally derived permeabilities (160.3, 43.1, and 3.9 nD). A good correlation of permeabilities modeled by the CTM compared to measured permeabilities by gas permeation was also found by Philipp et al. 2017 [51].

5. Conclusions

(i)BIB-SEM is a powerful tool that can deliver meaningful permeability values parallel to bedding from 2D pore areas based on simple capillary tube models as shown in this study for the heterogeneous Newark Shale.(ii)Upscaling pore system characteristics and permeability predictions from μm-sized mosaic maps to cm-sized plugs is achieved by assessing the mineralogy and associated pore properties of laminae with distinct characteristics individually via BIB-SEM and EDX.(iii)BIB-SEM derived permeabilities of the heterogeneous Newark Shale on the plug scale (upscaled) are representative and close to experimentally measured permeabilities.(iv)Silty laminae and OM porosity control porosity and permeability of the lacustrine Newark Shale.


A. Pore Parameters for the Statistical Evaluation

Pore elongation is based on a pore’s axial ratio () after assessing the “minimum bounding geometry” with the parameters width () and length ():The circularity describes the pore area in relation to its circumference and is a geometrical indicator for the overall roundness. It is expressed as a function of pore area () and pore perimeter ():

B. Detailed Derivation of the CTM

Fluid flow through porous media is described by Darcy’s law for incompressible media:

Here, is the volume flow rate (m3 s−1), the cross-sectional area (m2) of the porous medium, the permeability (m2), the dynamic viscosity (Pa s), and the pore pressure gradient (Pa m−1) in the direction of fluid flow.

The Hagen–Poiseuille equation describes laminar flow through a long cylindrical pipe of constant cross-section (a capillary tube):Here, is the hydraulic radius of the capillary, the dynamic viscosity (Pa s) of the permeating fluid, and the pressure difference (Pa) over the length .

Combining Darcy’s law (see (B.1)) with Poiseuille’s law (see (B.2)) delivers

With the assumption of a sample cube of the length and the cross-sectional area across which the pressure gradient is applied, the 3D-porosity of a single capillary tube with the length is defined as

Solving (B.4) for the area and combing it with (B.3) enable us to express the permeability asIncluding the tortuosity ( = ) in (B.5) as well yields a permeability coefficient of a single tortuose capillary pathway:

Porosity of a sample cube with a single capillary can also be derived from the 2D BIB-SEM data by using areas instead of volumes and thereby neglecting the 3rd dimension:Here, is the pore area (m2) and the total area of the mosaic map (m2). Inserting (B.7) into (B.6) yields

The hydraulic radius can also be calculated from the BIB-SEM data byThe total permeability (m2) of a specific area, for example, a complete SE2 mosaic map , is then calculated by adding up all single permeability values of each pore:

Conflicts of Interest

The authors declare that there are no conflicts of interest regarding the publication of this paper.


The authors thank P. Bertier of the Clay and Interface Mineralogy (CIM) research group at the RWTH Aachen University for providing them with the thorough XRD analysis. P. E. Olsen is acknowledged for his allowance to sample the Newark Basin Coring Project cores and his help in the IODP core store at Rutgers University.