Quantifying Centimeter- to Microscale Heterogeneities in Sedimentary, Compositional, and Geomechanical Properties of Siltstone Deposits in the Lower Triassic Montney Formation, Northeastern British Columbia, Canada

Although often described as homogeneous, marine mudstone and siltstone deposits in the ancient sedimentary record exhibit compositional, mechanical, and sedimentological fabric heterogeneities from the millimeter to the kilometer scale [1]. However, our current understanding of the heterogeneities in sedimentary fabric and compositional properties at the micro- (<mm)- to centimeter (cm)-scale and how they affect the mechanical properties in unconventional reservoirs remains poorly understood (e.g., [2–4]). Studies that have examined at fine-grained deposits in the subsurface and outcrop have recognized the complexity of centimeter-scale facies and their respective microscale characteristics that often comprise unconventional fine-grained reservoirs (e.g., [5–8]). One of the key parameters that is utilized from these facies is their mineralogical characteristics, which are linked to the mechanical properties of the rock, namely, the elastic moduli which determines the brittleness of the rock [9]. For example, studies in well-known unconventional reservoirs, such as the Barnett shale, Eagle Ford, Woodford, and Duvernay formations, have established the use of X-ray diffraction (XRD) to define a brittleness index [10–12]. This index uses mineralogical abundances of quartz, carbonate, clay, and in some equations total organic content (TOC) as a proxy to understand the influence that composition has on brittleness. However, it is important to note that each of the above formations comprises distinct geomechanical units due to contrasting compositional make up (e.g., carbonate-rich, chalk, organic-rich, and siliciclastic-rich facies) [10, 13, 14]. On the other hand, laboratory studies that focus on destructive methodologies to measure mechanical properties by using core plugs to run point load strength, triaxial, and brazilian tests provide insight on maximum principal effective stress, but significantly limit the resolution of samples to <2–3 centimeters (e.g., [11, 15–18]). Lastly, the most common, cost-effective routine for core analysis is sampling for composition using a handheld X-ray fluorescence (XRF) gun and a handheld hardness (HLD) piccolo tool, respectively. These methodologies sample along set intervals of every few centimeters to meters along entire cores (10’s to 100’s of meter), which offers insight into meter-scale mechanical and compositional heterogeneities [14, 19–21], but do not address how sedimentary facies at the centimeter and microscale might affect mechanical properties.

To overcome some of the limitations of the above methodologies and studies, namely, the characterization of small-scale heterogeneities of fine-grained rocks, this study addresses the influence that centimeter- to microscale heterogeneities have on mechanical properties and natural fractures in one of Canada’s most proliferous unconventional reservoirs, the siltstone reservoirs of the Lower Triassic MontneyFormation (northeastern British Columbia, Canada). The objectives of this study are (I) to characterize the sedimentary fabric, composition, and mechanical properties through a centimeter-scale gridding across a range of fine-grained sedimentary facies, (II) qualitatively document the microfabric and elemental distribution using a scanning electron microscope (SEM), (III) characterize observed natural calcite-filled fractures in the core, and (IV) explore the control that sedimentary fabric and composition have on mechanical and natural fracture properties. Collectively, this dataset aims to improve our understanding on the relationship between sedimentary fabric, composition, and geomechanical properties bellow well log resolution and a novel approach of integrating multiple datasets to understand the characteristic heterogeneities of fine-grained deposits in unconventional reservoirs.

The Lower Triassic Montney Formation of Western Canada is a westward thickening (up to 350 m thick) mixed siliciclastic-carbonate wedge, which extends from Alberta to British Columbia (Figure 1(a)). The formation comprises three formal subdivisions, the Lower, Middle, and Upper Montney members, which are bounded at the base by the Belloy Formation and the Sunset Prairie/Doig formations at the top ([22]; Figure 1(b)). Triassic strata were deposited along the western margin of the Pangea supercontinent during the Early Triassic (251.9–247.1 million years ago) [23–26]. During this time, the western margin of Pangea consisted of arid coastlines based on aeolian dunes, extensive anhydrite beds, and solution collapse breccia resulting from evaporate dissolution in upper Middle to Upper Triassic marginal marine and nonmarine successions [27–29]. Moreover, the paucity of detrital clay and presence of detrital feldspars have been attributed to a lack of chemical weathering on the continent, providing further evidence for an arid source area [23]. The nearly uniform grain size, ranging from silt to very-fine grained sand, suggests aeolian transport by suspension and saltation played a major role in transporting the sediments to the marine realm [22, 23], adding further support to an arid climate.

Triassic strata of the Western Canadian Sedimentary Basin were originally thought to have been deposited in a passive continental margin setting, with sediment being supplied exclusively from the North American craton [27]. However, recent regional tectonic reconstruction suggests instead that the North American craton was most likely tectonically active during this time with allochthonous terranes inferred to be an island arc basinward of the western margin of Pangea. This suggests that the Montney Formation may have likely been deposited in a back-arc or collisions retro-foreland basin setting with some sediment input coming from the western margin [22, 30].

The Montney Formation composition consists of siliciclastic-, bioclastic-, dolomitic-, and phosphatic-rich units, which are interpreted to be deposited primarily along a low-relief, mixed siliciclastic-carbonate ramp from shoreface to offshore [22, 23, 31, 32]. However, depositional environments associated with the Montney Formation range from turbidite succession (channel and fan complexes), storm- and wave-dominated shorelines, deltaic, and estuarine successions to biostromes (e.g., [22, 23, 32–38]).

The study area is within the “Northern Montney Field” ([39, 40], which is located in northeastern British Columbia, Canada (Figure 1). Two subsurface drilled cores were studied and totaled 145 meters in length. Note that core locations and names cannot be disclosed; however, they will be referred as well A and well B, respectively. The approximate locations of the wells are shown in Figure 1. Cored intervals and analyzed samples are from the Upper and Middle Montney members (Figure 2).

A sedimentological, compositional, and geomechanical analysis was conducted on drill core samples to describe rock heterogeneity. Data was collected from 10 samples to quantify sedimentary fabric (using HU values from CT scans), elemental composition (using XRF), and hardness (using Equotip Piccolo Bambino handheld microhardness). Each sample was characterized for the sedimentary fabric, composition, and mechanical properties using a 1 cm-by-1 cm grid (cf. [3]). See supplementary material 1 for workflow on micro- and centimeter-scale heterogeneity characterization used in this study and supplementary material 2 for equipment used to quantify elemental and geomechanical data.

At the centimeter-scale, samples were described based on the lithology, grain size, and sedimentary structures. Additionally, high-resolution (0.5 mm), nondestructive, X-ray computed tomography (CT) for imaging was used to characterize the fabric and structures in the representative core samples, with CT scans providing qualitative values to aid with centimeter-scale descriptions [41, 42]. Computer tomography scans rely on the variation of X-ray attenuation within the sediment cores which are expressed as Hounsfield units (HU; [43]). Hounsfield units are a function of bulk density, which is reflective of grain size and mineralogy [42]. Therefore, HU can be used as a proxy to numerically characterize sedimentary fabric; lighter, coarse-grained, more cemented layers have higher HU values, whereas darker, finer-grained, less cemented, and more mica-rich layers have lower HU values. Using ImageJ, a total of five Hounsfield units (HU) were average from a circle area of 1 centimeter to minimize human error. Maximum and minimum values were removed from the average calculation.

To characterize chemical composition, a Brucker Tracer handheld X-ray fluorescence (XRF) gun with a spot size of 8 mm was used on each sample. X-ray fluorescence analysis yielded a spectrum of elemental data for 877 data points. Calcium (Ca), magnesium (Mg), strontium (Sr), and manganese (Mn) have been used as proxy for carbonate, whereas aluminum (Al) was used as clay/mica volume. Iron (Fe) and sulphur (S) denote pyrite occurrence, silicon (Si) is used as a proxy for quartz, and the ratio of Si/Al is used as proxy for biogenic silica (cf. [44–46]). Elemental values were then normalized by not considering the light elements output from the XRF analysis due to variations in light elements between samples making up 57–67% of the bulk composition. This allowed for the direct comparison of elemental data points across the different facies [3]. For quality control, XRF values that diverged from the surrounding XRF points were ran 2-3 times to assure that XRF values were anomalous rather than associated with data acquisition error. Thin section samples were used for elemental mapping using energy dispersive spectroscopy (EDS) on the SEM to document the distribution of bulk composition (e.g., Ca, Al, and Si), qualitatively characterize the mineral phases, and describe digenetic properties of the samples using backscatter electron scanning for each microfacies.

An Equotip Piccolo Bambino handheld microhardness tool with a 3 mm diameter impact area was used to collect hardness data form the drill core samples. Data were collected using the same $1cm×1cm$ gridding system, and a total of 4760 microhardness tests were obtained. The Piccolo handheld microhardness tool measures the rebound of a small metal ball off the surface of the rock and records impact velocity and rebound velocity and gives an output value of HLD (handheld hardness). Six hardness readings were taken per point in which the first four were averaged and considered to represent the single grid hardness point. The first four points were averaged as after the fourth measurement HLD values remain consistent and reach a plateau [47]. Previous studies have interpreted that hardness correlates to unconfined compressive strength (UCS) in the Montney [19, 48] and, therefore, provides a cost-effective, nondestructive proxy to quantify geomechanical properties in the Montney Formation.

Fractures were characterized from 145 meters of unoriented slabbed cores and were subdivided based on geometry. Fracture geometries included brecciated, bed-parallel, and subvertical to vertical (Figure 2). Additionally, nature of termination, thickness, fracture aperture, and bounding fabric was recorded.

To explore and understand the relationship between sedimentary fabric, elemental composition, and microhardness, a discriminant analysis (DA) was used. Discriminant analysis maximizes the discrimination between predefined categorical groups (i.e., microfacies) based on their linear combinations of multiple numerical variables (i.e., elemental composition, fabric, and hardness) [49]. This classification technique can be used as an explanatory framework (i.e., to find trends) and as a predictive framework, in which the relationship between the input data is used to predict unknown datasets [50]. In this study, only the explanatory framework of a discriminant analysis was used. A discriminant analysis identifies linear combination of multiple variables, which best separates each categorical group. The data was run using the linear discriminant analysis package on XLSTAT (Addinsoft, New York, USA) and transformed into a single discriminant score in which the graphical representation is optimal, when more than 80% of the dataset can be accounted for. The resultant plot, therefore, explains >80% of the data and maximizes the component axes for greatest class-separation (i.e., what makes each microfacies unique).

To reduce the dimensionality of the dataset in XLSTAT, a total of eight computational steps take place: represent dataset as a matrix (eq. (1)), compute the mean for all features (i.e., elemental composition, fabric, and hardness) for each class (eq. (2)), compute the total mean of all data (eq. (3)), calculate the between-class matrix (eq. (4)), compute within-class matrix (eq. (5)), calculate a new transformation matrix using between- and within-class variance (eq. (6)), calculate the eigenvector and eigenvalues, and transform the original data using eq. (7). Note that the software calculates all these parameters automatically from the input dataset. Here, however, an overview is provided of the equations and description for each parameter that are in the software algorithm for the DA. For more detail on the mathematical procedure refer to Balakrishnama and Ganapathiraju [51] and Tharwat et al. [52], and references therein. Figure 3 provides a visualization example of all calculated parameters described below.

In this equation, $xij$ represents the $i$^th sample in the $j$^th class. For example, a hardness measurement (⁠$i$^th) within a microfacies and compared to the hardness mean of that specific microfacies (⁠$j$^th).

Using the transformation matrix, eigenvalues $λ$ and eigenvectors (⁠$V$⁠) are then solved to form the lower dimensional space in which the original dataset is projected onto. Notably, eigenvalues and eigenvectors are the parameters that provides the information regarding the DA space [52]. The eigenvectors represent the directions where the data is projected onto, whereas eigenvalues represent the scaling factors, length, and magnitude of the eigenvalues (i.e., how well that specific axis allows for the discrimination between microfacies). Therefore, a higher eigenvalue corresponds to a more robust ability of the model to differentiate between the classes (microfacies).

A total of 10 representative core samples were analyzed and consist of three primary facies: massive siltstone (⁠$F1$⁠; $N=2$⁠), stratified (siliciclastic) siltstone (⁠$F2$⁠; $N=4$⁠), and bioturbated siltstone (⁠$F3$⁠; $N=4$⁠), which represent the majority (91%) of centimeter-scale heterogeneity within the cored intervals (Figures 2 and 4). The other 9% consists of stratified (mixed siliciclastic-carbonate) siltstone and concretions. Only siliciclastic-rich intervals were analyzed to compare between similar lithological and textural characteristics. $F2$ and $F3$ are the main reservoirs of interest in the Montney Formation, and therefore, more samples were analyzed from these two facies. Each facies can consist of multiple microfacies (MF), which have been simplified based on the dominant microfabric, and include massive siltstone (MF1), pinstriped laminated siltstone (MF2), planar- to cross-stratified siltstone (MF3), and bioturbated siltstone (MF4) (Figure 4). All obtained data is summarized in Table 1 and displayed as histograms in Figure 5. Note that only one gridding per facies is illustrated (Figures 6–8); however, the data obtained for each grid can be found in the supplementary data.

Intervals of $F1$ comprises massive, structureless, siltstone with rare, and isolated coarser laminas (⁠$MF1$⁠; Figures 4 and 6(a)). Quantified sedimentary fabric using Hounsfield units (HU) from the CT scans suggests a homogeneous rock fabric with a narrow HU measurement range of 2692–2850 (Figures 5(a) and 6(a) and 6(b) and Table 1). Intervals assigned to $F2$ comprise of laterally discontinuous, pinstriped laminae siltstone (MF2) stratified with 1–3 cm thick planar- to cross-laminated siltstone (MF3, Figures 5(a) and 7(a) and 7(b)). Quantified sedimentary fabric using HU is subdivided into two ranges—MF2 exhibits a range of 2664–2971, similar to MF1, whereas MF3 has a HU range of 2676–3110 (Figures 5(a) and 5(b) and Table 1). Intervals assigned to $F3$ are characterized by bioturbated siltstone (Figures 5(a) and 7(a) and 7(b)). Bioturbation index (BI) is 2-3 (cf. [53]), and fabric is moderately mottled and disorganized, except for rare undisturbed planar- to cross-laminae (MF3 and MF4; Figures 8(a) and 8(b)). The quantified sedimentary fabric HU measurement ranges from 2612–3110.

Irrespective of the sedimentary fabric, the most common compositional elements are Si, Ca, and Al, which make up >85% of the total elemental volume. The massive siltstone of $F1$ exhibits a Si range of 57.2–60.8% (avg. 59.3%), Ca concentrations ranging 10–13.8% (avg. 11.6%), and Al ranging 9.8–12.1% (avg. 11%). The stratified siltstone of $F2$ is characterized by a Si range of 32–63% (avg. 56%), Ca ranging 10–23% (avg. 15%), and Al ranging 7–12% (avg. 10%). Lastly, the bioturbated siltstone of $F3$ exhibits Si concentrations of 32-63% (avg. 51%), Ca ranging 14–51% (avg. 20%), and Al ranging 6–11% (avg. 8%), see Figures 5(b), 5(c), 5(d), 6(c), 7(c), and 8(c) and Table 1 for more detail.

Like compositional characteristics, mechanical properties also exhibit intra- and intersample variations in hardness values (Table 1 and Figures 5(e), 6(d), 7(d), and 8(d)). The massive siltstone of $F1$ has HLD values ranging from 619–804. The stratified siltstone of $F2$ has HLD values ranging from 572–849. Microfacies that make up $F2$ (MF2 and MF3) consist of a bimodal, broad hardness distribution. Lastly, the bioturbated siltstone of $F3$ exhibits HLD values ranging from 621–855 (avg. 780). Both $F1$ and $F3$ exhibit a unimodal, narrow hardness distribution.

Based on elemental mapping of microfacies and their respective microfabric, all microfacies comprise quartz, dolomite, and plagioclase (K- and Na-rich) and variation in authigenic calcite cement and pyrite. Energy dispersive spectroscopy also allows for the indirect relationship between the mineralogical composition and the elemental proxies used in this study.

Microfacies 1 is the most abundant in massive siltstone (⁠$F1$⁠) and less commonly in stratified siltstone (⁠$F2$⁠) and bioturbated siltstone (⁠$F3$⁠). Energy dispersive spectroscopy under the SEM exhibits detrital and overgrowth cement of quartz (avg. 35%), dolomite (avg. 30%), and feldspar (K- and Na-rich; avg. 10 and 12%, respectively), whereas the matrix component (18%) comprises clay-sized dolomite, quartz, and feldspars in addition to abundant micas (primarily muscovite). Both euhedral and framboidal pyrite occur as authigenic minerals (avg. 5%) and occur in open pore space, as inclusions in calcite cement, or more commonly between muscovite grains. Additionally, authigenic pore-filling calcite cement is uncommon (avg. 8%; Figure 9(a)).

Microfacies 2 and 3 consist of siltstone with planar laminae and are differentiated based on the grainsize of the dominant, thicker lamina at the microfacies scale. Microfacies 2 comprises of thicker, dark, finer-grained laminae (lamina thickness 200–700 μm), and thinner, light, coarser-grained laminae (120–250 μm). Conversely, MF3 comprises of thick, light, coarser-grained laminae (300–400 μm), whereas the dark, finer-grained laminae are thinner (100–200 μm). Irrespective of the thickness, in both microfacies, a light lamina comprises silt-sized quartz, dolomite, and feldspar (K- and Na-rich). Light laminae are well-cemented and comprises intergranular calcite cement (avg. 15%) in addition to quartz (avg. 35%), K- and Na-rich feldspar (avg. 15% and 15%, respectively), and dolomite with Fe-rich or Fe-poor cement overgrowth (avg. 15%). Note that light laminae lack appreciable matrix component (<fine-grained silt; 2%) of micas, in addition to pyrite (<1%; Figures 9(b) and 9(c)). Conversely, dark laminae comprise silt-size quartz (30%), dolomite (10%), and feldspar (K- and Na-rich, 12% and 12%, respectively) and less calcite cemented (⁠$avg.<3%$⁠; Figures 9(b)–9(d)). Dark laminae are distinctively matrix-rich (avg. 18%) comprising clay-sized quartz and dolomite in addition to micas (primarily muscovite). Both euhedral and framboidal pyrite are also common (avg. 5%).

Microfacies 4 dominates $F3$ and is intercalated with MF1 and MF4. Like the previous microfacies, the main components are silt-sized quartz (avg. 40%), dolomite (avg. 22%), and K- and Na-rich feldspars (avg. 8%). In MF4, clear contrasting laminae (matrix-rich intercalated with matrix-poor) are not observed. Instead, localized, well-cemented, pore-filling calcite cement (avg. 12%) in addition to quartz and dolomite cement overgrowth gradationally grade both vertically and laterally to matrix-rich (avg. 10%) zones that collectively form a mottled or disorganized textural fabric in MF4 (Figures 9(e) and 9(f)).

A total three eigenvalues were determined by XLSTAT and discriminated between the four microfacies. Computed discriminant analysis explains 95% of the data, LDA axis 1 accounts for 79% $λ=3.458$ of separability between microfacies, LDA axis 2 accounts for 16% $λ=0.657$ of the separability, whereas LDA axis 3 only accounts for 5% $λ=0.258$ of discrimination. Therefore, $F1$ and $F2$ were used to form the lower dimensionality plot as they corresponded to the highest eigenvalues (Figure 10(a)). Based on the discriminant axes (Figure 9(a)), numerical data used for this analysis exhibit two eigenvectors (arrows) directions, towards the northwest and southwest quadrants. Eigenvectors that have a small angle between them exhibit a positive trend, whereas vectors with larger angles between them exhibit negative trends. Therefore, normalized Ca%, Si/Al ratio, Hounsfield units (HU), and HLD all exhibit positive trends, whereas negative trends are observed with Al%, Fe%, and S%. The eigenvector for normalized Si% is towards the northeast (i.e., neutral position), which likely reflects the similarity in quartz abundance irrespective of the microfacies (Figure 5(b) and Table 1) and, therefore, does not provide a good elemental proxy to understand the hardness distribution between facies at least within these samples in the Montney Formation. Based on microfacies (Figure 10(b)), three clustering can be observed between the four microfacies (see Section 5.1 for more detail). Massive siltstone (MF1) and pinstriped laminated siltstone (MF2) make up a cluster with the lowest hardness values (Figure 10(b)). Planar to cross-stratified siltstone (MF3) exhibits the largest spread of data across the plots with the highest hardness values, whereas bioturbated siltstone (MF4) consists of a tight clustering and medium to high hardness values (Figure 10(b)).

The accuracy of the discriminant analysis can be assessed through a confusion matrix whereby the original dataset is compared to the transformed dataset in the discriminant analysis (Table 2) to test how well the model is at grouping the datapoints into the predefined groups and the liner combination of all quantitative variables. Based on the confusion matrix, 93% of the original dataset were correctly transformed and remained in the same microfacies post transformation to the lower dimensional space. Massive siltstone (MF1), pinstriped laminated siltstone (MF2), and bioturbated siltstone (MF4) had over 95% accurate transformations, whereas planar- to cross-stratified siltstone (MF3) had 80% of accurate transformation into the lower dimensional space.

A total of 75 natural calcite-filled fractures were documented within the studied cores (Figure 2); subdivided into three types: brecciated (⁠$n=7$⁠), bed-parallel (⁠$n=43$⁠), and vertical to subvertical fractures (⁠$n=25$⁠) (Figure 11); and described below and summarized in Table 3. Well A has fractures within the Upper Montney Member, whereas well B only has natural calcite-filled fractures in the Middle Montney Member.

Brecciated fractures (Figure 11(a)) range in length from 3–9 cm. Apertures range from 1.2–4.8 mm with an average of 2.8 mm. Brecciated fractures occur in pinstriped laminated siltstone (MF2) and a bioturbated siltstone (MF4).

Bed-parallel fractures range in length from 0.5 cm to the width of the core (9 cm) and most likely beyond the length of the core width; therefore, these lengths are noted to be minimum. Apertures range from 0.5–6 mm with an average of 1.8 mm. Bedding plane surfaces commonly exhibit drusy calcite or to patchy calcite mineralization (Figure 11(b)). Bed-parallel fractures notably occur in massive siltstone (MF1; Figure 11(c)) and along the interface of fabric change from MF1 to MF2 vice versa in $F2$ (Figure 11(d)).

Vertical to subvertical fractures range in length from 3.3–40 cm with an average of 12 cm. Apertures range from hairline to 0.5–5.3 mm with an average of 1.6 mm. Notably, vertical fractures exclusively occur in bioturbated siltstone (MF3), except for one occasion in which it occurs in a pinstriped laminated siltstone (MF2) and is either bounded by a change in sedimentary fabric, where they intersect bed-parallel fractures or taper within the same sedimentary fabric (Figure 11(e)). Additionally, vertical fractures occur along with postdepositional, small-scale, thrust, and normal faults that were observed only in well A (Figure 11(f)).

Based on the sedimentary fabric, compositional, and geomechanical properties described in the gridding analysis above, both inter- and intrasample heterogeneities are observed. Results from the discriminant analysis further quantify the relationship between sedimentary fabric and composition to mechanical hardness by clustering based on microfacies, rather than by individual samples. Massive siltstone (MF1) primarily occurs in $F1$⁠; exhibits on average the lowest HU values (avg. 2736), normalized Ca (avg. 12%), Si/Al (avg. 5), highest normalized Al (avg. 11%), Fe (avg. 6%), and S (avg. 4%); and collectively exhibits the lowest hardness (avg. 742) (Table 1). Figure 9 further illustrates this trend, where micas (Al is corresponding proxy) and pyrite (Fe and S are corresponding proxies) are much higher in MF1 as compared to MF3 and MF4. Moreover, pinstriped laminated siltstone (MF2) comprises similar HU (avg. 2721) to MF1, normalized Ca (avg 15%), Si/Al ratio (avg. 5.9), normalized Al (avg. 10%), Fe (avg. 6%), and S (avg, 2%), and similar hardness values (avg. 768). However, differences arise in the distribution of HLD, in which MF1 consists of a narrow, unimodal distribution, whereas MF2 broad, bimodal distribution (Figure 7) reflects data acquisition in both dark and light laminae in MF2 (Figures 5 and 7). Planar- to cross-stratified siltstone (MF3) exhibits the highest HU values (avg. 2851), Ca (avg. 18%), Si/Al ratio (avg. 6.5), and lowest Al (avg. 9.6%) and has the highest hardness (avg. 798). Bioturbated siltstone (MF4) has lower HU values (avg. 2808) and slightly lower normalized Ca (17%), Si/Al ratio (avg. 6) and Al (avg. 9%), Fe (avg. 5%), and S (1%). Only one sample that mainly comprises MF4 exhibited higher HU (avg. 2984) and normalized Ca (avg. 45%) with the highest hardness in contrast to other MF4 samples (avg. 804).

In this study, HU values from CT scans corroborates with both the handheld XRF and the petrographic SEM characteristics and provides a novel approach to quantify sedimentary fabric. Higher HU occurs in microfacies that exhibit higher normalized Ca% (calcite and dolomite cement), lower Al% (micas), and are overall comprises coarse silt (e.g., MF3 and MF4). Lower HU occurs in microfacies that exhibits lower normalized Ca% (negligible calcite cement), higher Al% (micas), and is overall fine to medium silt (e.g., MF1 and MF2). Notably, the abundance of the matrix component (micas and clay-sized quartz, dolomite, and feldspars) negatively correlates to authigenic, pore-filling calcite cement in addition to how well-connected overgrowth cement, which likely play a role on the hardness values for each microfacies.

As described above, elemental composition is similar between massive siltstone (MF1) and pinstriped laminated siltstone (MF2) and between planar- to cross-stratified siltstone (MF3) and bioturbated siltstone (MF4) (Figures 5(b)–5(d)). However, the hardness distributions are dependent on the microfacies (Figure 5(e)). Two endmembers can be discerned when comparing the dispersion of hardness values (Figure 5(e)). Hardness distribution in MF1 and MF4 is narrow (values are close to the mean), and therefore, suggesting that centimeter-scale facies $F1$ and $F3$ are more homogeneous (Figure 4) and less likely to have significant mechanical interfaces. This contrasts the distribution of hardness values in MF2 and MF3 (stratified siltstone; $F2$⁠), which have a similar broad distribution, and therefore, mechanical interfaces in this facies occur due to differences in laminae properties (i.e., laminae thickness, cementation proportion, and matrix variation). These mechanical interfaces are critical to understand as they can impose mechanical barriers that could inhibit fracture propagation in reservoir intervals while undergoing hydraulic fracturing. Additionally, the interbedding between more mechanically homogeneous facies (⁠$F1$ and $F3$⁠) with more mechanically heterogeneous facies (⁠$F2$⁠) could further influence fracture characteristics (height, arrest, aperture, etc.) at the meter- to kilometer-scale.

Hardness distributions as a rock mechanical proxy are also supported by the characteristics of natural calcite-filled fractures in the different sedimentary facies in the studied cores (Figure 11), which highlight the relationship between mechanical and fabric heterogeneity and collectively illustrate how the different fabrics influence natural fracture progradation and potentially under hydraulic stimulation. For example, fractures in massive siltstone (⁠$F1$⁠) are preferentially bed-parallel, but only occur in isolated light, coarser-grained laminas. Whereas fractures in stratified siltstone (⁠$F2$⁠) are bed-parallel and bounded by the abrupt transition between the dominant microfacies, namely, pinstriped laminated siltstone (MF2) and planar- to cross-stratified siltstone (MF3). Fractures in bioturbated siltstone (⁠$F3$⁠) are dominantly vertical to subvertical. Both centimeter- and microscale analyses suggest that the facies and microfacies matrix composition, intergranular cement (carbonate), and overgrowth cement (dolomite and quartz) in addition to the degree of fabric interbedding play a primary role on the variability in mechanical hardness and the geometry and termination of natural calcite-filled fractures along the cored intervals.

Similar findings have been reported on the importance of facies interbedding, mechanical boundaries, and fracture growth patterns in the bioclastic-rich Altares Member of the Middle Montney Member [36], in which contrasting bitumous siltstone and bioclastic packstone/grainstone form stair-step fracture growth patterns, whereas geomechanical homogeneous bituminous siltstones consist of simple dendritic fracture networks. Similarly, to this example from the Altares Member, centimeter-scale facies and natural calcite-filled fractures found in this study further suggest that facies interbedding plays a critical role in fracture geometry and termination.

The geomechanical properties have been well documented in other fine-grained formations in the US and Canada (Figure 12). Hardness values from this study align with reported values from the Montney Formation in other studies (e.g., [19, 54]). As described above, hardness distributions vary within the microfacies. Overlap between sedimentary fabric and their respective geomechanical properties is a consequence of diagenetic homogenization and similarities between detrital composition (silt- and clay-size quartz, dolomite, and feldspar) in addition to pore-filling calcite cement, quartz, and dolomite overgrowth cement. Hardness values in this study are high, ranging from 572–855 HLD. This is likely affected by the presence of the calcareous cement even in the finer, more matrix-rich microfacies (e.g., MF1 and MF2) and presence of pyrite, which collectively enhance the hardness of the rock. In contrast to other fine-grained formations, data from this study plot in the upper half on the hardness axis, like calcareous-, pyritic-, and dolomitic-rich facies in other formations. The two dominant eigenvector directions of the discriminant analysis (Figure 10) in which higher normalized Ca% and Si/Al have positive trends with HLD values, whereas normalized Al% exhibits a negative trend with HLD have also been reported in other fine-grained formations [1, 10, 55]. Notably, the Woodford and Perdrix formations comprise hardness that is facies-dependent (Figure 12), and this is likely due to contrasting compositional make up (e.g., carbonate-rich, organic-rich, or siliciclastic-rich) in these formations [10, 13, 14], which is unlike the Montney or other formations as detrital and diagenetic components (i.e., differences in authigenic minerals) between facies may be subtle.

One key aspect to note also is the sampling methodology across the different studies, and current methodology consists of sampling 10’s–100’s meter along the core for each respective formation (except for the Perdrix Formation study, which sampled representative centimeter blocks every few meters in outcrop). While this methodology offers insight into meter-scale mechanical and compositional heterogeneities, a recent study by Venieri et al. [3] on the Perdrix Formation (Duvernay outcrop-equivalent) in addition to this study have demonstrated that even at the micro- to centimeter-scale, there is a significant under sampling in unconventional reservoirs due to the fine-scale fabric. Well-established relationships between mechanical and compositional properties (e.g., [1, 56]) would suggest that intuitively, more sampling should lead to narrowing and clustering of the data, especially when documenting facies. However, broad distribution of hardness values and overlapping irrespective of facies occurs in almost all formations (Figure 12). As pointed out by Venieri et al. [3], these likely occurs due to the higher resolution data acquisition (increase in sampling points), which, instead, leads to a higher dispersal of the data and, therefore, reduces the strength of the relationship between elemental composition, hardness, and facies as more heterogeneity is documented. This is likely due to microfacies heterogeneities that can comprises centimeter-scale facies and the introduced intrasample heterogeneity that occurs from the spatial resolution of data acquisition. For example, handheld XRF obtained data at 8 mm radius and thereby averaged the bulk composition and included multiple microfacies, whereas hardness values were collected at 3 mm radius and are likely representative hardness of an individual microfacies. Another potential explanation for the reduction of strength of relationship between facies, elemental composition, and hardness could arise from the occurrence of the different mineral composition (i.e., matrix, cement, or as framework grains) and the interaction between the occurrence of the minerals and the handheld piccolo tool. Nevertheless, trends and relationships from the discriminant analysis in conjunction with microscale characterization suggest that mechanical and compositional heterogeneities are facies- and microfacies-dependent and highlight the importance of understanding the small-scale heterogeneities of unconventional fine-grained reservoirs below well log resolution. However, the dispersion of both mechanical and elemental data (Figure 5 and Table 1) obtained from the centimeter gridding analysis raises the question of how accurate upscaled models are when only averaged values are used to model unconventional fine-grained heterogeneous reservoirs.

Additionally, studies focusing on the impact that composition (XRF and/or XRD) has on hardness commonly use scatter plots for trends between hardness and composition (e.g., [3, 55, 57–59]). Figure 13 illustrates the most common elements (Ca, Al, Si, and Si/Al) plotted against hardness and general trends can be discerned. For example, hardness has a positive relationship with Ca% and Si/Al, negative trend with Al%, and a neutral trend with Si%. However, with the added complexity of the microfacies in these plots, there is significant overlap, and relationships between microfacies, elemental composition, and hardness are not straightforward. This is problematic as this study demonstrated that differences in composition and sedimentary fabric play a role in mechanical boundaries and natural fracture termination styles, based on distribution plots of microfacies (Figure 5) gridding analysis (Figures 6–8), elemental mapping under the SEM (Figure 9), and natural fracture characterization (Figure 11). Therefore, there is a need to quantify and understand these elemental and mechanical relationships while also considering fabric heterogeneities in unconventional reservoirs. In this study, applying a discriminant analysis (Figure 10) to illustrate sedimentary fabric, elemental composition, and hardness data provides a novel method to quantify the relationship between geological and geomechanical properties.

Ten samples from the unconventional siltstone reservoirs of the Triassic Montney Formation in Northeastern British Columbia, Canada, were analyzed to understand the influence that sedimentary fabric and composition have on geomechanical properties. Three facies were analyzed and comprise a range of centimeter-scale heterogeneities observed in the study area and included massive siltstone (⁠$F1$⁠), stratified siltstone (⁠$F2$⁠), and bioturbated siltstone (⁠$F3$⁠).

• (1)

  Based on a centimeter-scale gridding, sedimentary fabric, composition, and geomechanical properties were examined, and significant inter- and intrasample heterogeneity exists in the Montney Formation. This sample heterogeneity is a result of microfacies interbedding, which included massive siltstone (MF1), pinstriped laminated siltstone (MF2), planar- to cross-stratified siltstone (MF3), and bioturbated siltstone (MF4)

• (2)

  Microfacies that are matrix-poor have well-interconnected, pore-filling calcite cement in addition to quartz, dolomite, and feldspar cement overgrowths (MF3 and MF4), whereas microfacies that are more matrix-rich have comparatively less cement (MF1 and MF2). This likely plays a primary role in the fabric, elemental composition, and hardness distribution obtained in centimeter-scale facies and collectively highlights the influence of fabric and composition in mechanical properties in the Montney Formation

• (3)

  Hardness data distribution is dependent on the microfacies and is a practical proxy to understand mechanical boundaries and fracture properties. A narrow hardness dispersal suggests negligible mechanical boundaries in a homogenous fabric and is thereby likely to be associated with tall, vertical fractures. Broader hardness dispersal, on the other hand, occurs in more heterogenous, interbedded fabric and suggests the presence of mechanical interfaces

• (4)

  Based on a discriminant analysis trend between sedimentary fabric, composition, and geomechanical properties, sedimentary fabric with the highest centimeter-scale aluminum (Al) content from XRF comprises microfacies that are comparatively matrix-rich that consist of mica and clay minerals (observed in SEM) and exhibit the lowest hardness values. In contrast, sedimentary fabric with a higher elemental calcium (Ca) component comprises microfacies that are matrix-poor, highly cemented by carbonate (calcite and dolomite) and quartz, and overall exhibit a positive trend with hardness measurement

• (5)

  A discriminant analysis proved to be a novel approach to illustrate the effect that sedimentary fabric and compositional have on mechanical properties in more compositionally homogenous formations. Future work could aim on exploring these relationships in more complex facies in the Montney Formation, in less compositionally homogeneous formations, and including other reservoir properties such as organic matter content, porosity, and permeability

• (6)

  Sampling representative facies and their respective microscale properties can be a useful first step that can be done in other fine-grained deposits, namely, in key reservoir facies or prior to upscaling

Collectively, this dataset improves our understanding of the relationship between sedimentary fabric, composition, and geomechanical properties bellow well log resolution and a novel approach of integrating multiple datasets to better statistically model and asses the heterogeneity of fine-grained deposit characteristics, which are key for development strategies and reservoir modelling.

The authors confirm that the data supporting the findings of this study are available within the article and the supplementary materials.

The authors declare no conflicts of interest.

This work was supported by sponsors of the Tight Oil Consortium (TOC) and a Natural Science and Engineering Research Council of Canada (NSERC) Alliance Grant (ALLRP 548652-19) to Dr. Chris Clarkson and Dr. Per K. Pedersen. Additional funding from a Grants-In-Aids 2021 AAPG J. Ben Carsey Sr. Memorial Grant and a NSERC CGS-D to the first author, Patricia Fraino, also allowed for this research to be conducted. Thanks are due to Dr. Amin Ghanizadeh and Dr. Adnan Younis from the TOC for their help with XRF and piccolo tool data collection, Dr. Chris Debuhr for data acquisition using the scanning electron microscope, and Dr. Marco Venieri for providing feedback on an earlier version of the manuscript.