Nanoindentation of Horn River Basin Shales: The Micromechanical Contrast Between Overburden and Reservoir Formations

We present a micromechanical characterization of shales from the Horn River Basin, NW Canada. The shales have contrasting mineralogy and microstructures and play different geomechanical roles in the field: the sample set covers an unconventional gas reservoir and the overburden unit that serves as the upper fracture barrier. Composition and texture were characterized using X‐ray diffraction, mercury injection porosimetry, and scanning electron microscopy (SEM). Grid nanoindentation testing was used to obtain the mechanical response of the dominant phases in the shale microstructure. Samples were indented parallel and perpendicular to the bedding plane to assess mechanical anisotropy. Chemical analysis of the grids with SEM‐EDS (energy dispersive X‐ray spectroscopy) was undertaken and the coupled chemo‐mechanical data was used in a statistical clustering procedure (Gaussian mixture model) to reveal the mechanical properties of each phase. The results show that the overburden consists of a soft clay matrix with highly anisotropic elastic stiffness, and stiffer but effectively isotropic inclusions of quartz and feldspar; the significant anisotropy of the overburden has been previously observed on a much larger scale using microseismic data. Creep displacement is concentrated in the clay matrix, which is the key phase for fracture barrier and seal applications. The reservoir units are harder and have more isotropic mechanical responses, primarily due to their lower clay content. Despite varied compositions and microstructures, the major phases of these shales (clay/organic matrix, quartz/feldspar, dolomite, and calcite) have unique mechanical signatures, which will aid identification in future micromechanical characterizations and facilitate their use in upscaling schemes.

Since the mid-2000s, researchers have begun to investigate the mechanical properties of the shale constituents, the grains of phyllosilicates, quartz, feldspar, carbonates, organic matter (OM) and pyrite that variously make up the micro-composite shale material. The microstructure is significantly heterogeneous and the scale of the shale "building blocks" ranges from sub-micron sized clay plates to lenses and particles of organic material at a few microns across to silt-size minerals. The experimental methods have included nanoindentation testing Charlton et al., 2021;Du et al., 2022;Goodarzi et al., 2017;Ortega et al., 2007;Ulm & Abousleiman, 2006;Ulm et al., 2007;Veytskin et al., 2017) and atomic force microscopy (AFM) (Charlton et al., 2021;Eliyahu et al., 2015;Emmanuel et al., 2016;Goodarzi et al., 2017;Khatibi et al., 2018). These micromechanical tests can be conducted on small pieces of shale, such as drill cuttings or core fragments, instead of relying on the recovery of intact core; this also facilitates the characterization of shale anisotropy.
Nanoindentation describes the technique of instrumented grid indentation, whereby an indenter is pushed into the shale in a grid pattern and mechanical properties extracted from the load-displacement curves (Constantinides et al., 2006;Fischer-Cripps, 2011;Ulm et al., 2007). The length-scale of the experiment depends mainly on the magnitude of the applied load, and in shale can target either the properties of individual minerals (with a maximum force of around 2-50 mN), or a homogenized response representative of the bulk material (100-500 mN) (Graham et al., 2021). In contrast, AFM involves rapidly tapping a probe over the sample surface with a peak force of ∼100 nN (Bruker, 2012), with the approach and withdrawal of the tip continuously tracked to obtain a high-resolution nanomechanical map. This technique has often been used to target specific, difficult-to-reach phases in shale, such as the OM (Charlton et al., 2021;Eliyahu et al., 2015;Fender et al., 2020).
One advantage of nanoindentation over AFM is that by holding the indenter in the material at maximum load for a period of time, creep behavior can be assessed. The geomechanical response of shale is well known to be time-dependent (Sone & Zoback, 2013a), meaning that under a constant stress the rock will undergo continuous deformation, closing fractures over time. In unconventional reservoirs, the implication is reduced production rates (Sone & Zoback, 2014), while in fluid storage applications creep deformation may result in self-healing of seal lithologies and prevention of leakage pathways, particularly around injection wells (Cerasi et al., 2017). The creep behavior of shale is complex, and in long term triaxial creep tests, the response can encompass several phases (primary-secondary-tertiary creep; Rybacki et al., 2017) and is also strongly affected by saturation, temperature, and fluid chemistry (Cerasi et al., 2017). Triaxial creep experiments take hours, days (Herrmann et al., 2020;Rybacki et al., 2017;Sone & Zoback, 2013b), or even years (Gasc-Barbier et al., 2004) to complete and, as a result, available data on the core-scale creep behavior of shale is limited. Nanoindentation creep is measured on much shorter time periods, generally 60-180s (Mighani et al., 2019;Slim et al., 2019). Many studies have shown that the creep deformation in both nanoindentation and in the primary creep phase of a triaxial test under constant deviator stress can be described by a power or log relationship. Indeed, Charlton et al. (2021) showed that the creep modulus of Posidonia shale measured from high-load nanoindentation testing is comparable to results from core-scale triaxial creep experiments (lasting up to around 1 day), suggesting similar mechanisms are governing the response, although the creep relationship across scales is still an area of research (K. Liu et al., 2021). Low-load nanoindentation tests allow for greater understanding of where the creep deformation is concentrated, for example, in the softer phases of clay and OM (Mighani et al., 2019), and hence how changes in mineralogy in the field may influence the potential for helpful self-sealing mechanisms to occur.
Correlations between mineralogy and macroscopic mechanical properties have been reported (Herrmann et al., 2018;Sone & Zoback, 2013a) and the bulk mineralogy often determines the phase contributions in homogenization models, where the aim is to fill the gaps in core-scale data by upscaling from the mechanical prop erties of the shale constituents (Abou-Chakra Guéry et al., 2010;Goodarzi et al., 2016;Hornby et al., 1994;Ortega et al., 2007). However, factors such as texture and diagenesis can also influence mechanical behavior without significantly altering the mineralogy Charlton et al., 2021). Constructing an appropriate upscaling model, which identifies and captures the governing mechanical processes, is therefore a significant challenge. Recent efforts have begun to account for the inherent variability of the shale response by propagating uncertainty through a multiscale homogenization model using a probabilistic approach (Dubey et al., 2019).
Despite the challenges of upscaling, micromechanical characterization enables fundamental insights into the response of shales, beyond that possible from standard core tests or wireline logs. This paper presents a comprehensive micromechanical characterization of several shale lithologies from the Horn River Basin. The shales have contrasting mineralogy and contrasting geomechanical roles in the field: the sample set covers both an unconventional gas reservoir and the overburden unit that forms an upper fracture barrier. The shale samples are first described in terms of mineralogy, porosity, and texture. Nanoindentation testing is then used to obtain the mechanical response of the shales at a microstructural level. The mechanical data is coupled with chemical analyses and the chemo-mechanical data is used to constrain the behavior of the individual material phases in terms of elastic stiffness and the time-dependent creep response.

Horn River Basin
The Horn River Basin is located in northeast British Columbia, Canada (Figure 1). It covers an area of approximately 12,000 km 2 in the northwestern part of the Western Canadian Sedimentary Basin where predominantly carbonates and marine shales were deposited in the Middle and Late Devonian (Oldale et al., 1994). As illustrated in Figure 1a, the carbonate Presqu'ile reef platform bounds the Horn River Basin to the south and east while to the west the Bovie fault (maximum displacement 1,200 m) separates the Horn River Basin from the Liard Basin (Ross & Bustin, 2008).
The Horn River shale sequence comprises the Evie and Otter Park members of the Horn River formation, and the Muskwa formation ( Figure 1b). Thermal maturity is in the dry gas window across the basin at around 2.5% Ro (vitrinite reflectance; Reynolds & Munn, 2010;Ross & Bustin, 2009), equivalent to a maximum burial temperature of approximately 200°C, and these formations have been exploited for unconventional gas production (BC Oil and Gas Commission, 2014). The Muskwa formation is a gray to black, siliceous shale that ranges in thickness from 30 to 80 m, thickest in the west toward the Bovie fault (Ross & Bustin, 2008). It is organic-rich, with Dong et al. (2015) reporting an average 3.4 weight percent (wt%) total organic carbon (TOC). The Otter Park member has less OM, averaging 2.4 wt% TOC, and is described as a dark-gray, non-calcareous to calcareous or siliceous shale (Dong et al., 2015(Dong et al., , 2017. It has a maximum thickness of 270 m in the southeast part of the Basin, thinning to the north and west (McPhail et al., 2008). The more clay-rich lithofacies of the Otter Park member have been reported to act as fracture barriers in hydraulic fracturing operations (Yu & Shapiro, 2014). The Evie member is 40-75 m thick, with very high TOC content: 3.7 wt% average TOC in Dong et al. (2015). It is a dark gray to black, calcareous, siliceous shale (McPhail et al., 2008).
The Horn River shale sequence is underlain by the shallow water carbonates of the Keg River Formation and overlain by the argillaceous, organic-lean shale of the Fort Simpson formation, which reaches thicknesses of over 1,000 m (Ross & Bustin, 2008). The Fort Simpson shale has a gradational contact with the Muskwa formation (McPhail et al., 2008) and has been shown to perform as an effective fracture barrier (BC Oil and Gas Commission, 2014). Microseismic data from the Horn River Basin also shows that the Fort Simpson shale has a very strong seismic anisotropy, with evidence of shear wave triplication (Baird et al., 2017).

Composition
Core was obtained from two wells (A-A100-B/094-O-09 and D-094-A/094-O-08) drilled in the eastern part of the basin (see Figure 1a). Overburden (Fort Simpson) samples were taken from well A-A100-B/094-O-09 from depths 2,537.4 m (sample A3) and 2,542.2 m (A6), approximately 50 m above the Muskwa formation. A sample of the Muskwa formation (D1) was taken from well D-094-A/094-O-08, at depth 2,489.25 m, while Otter Park and Evie samples were obtained from well A-A100-B/094-O-09 at depths of 2,665.25 (A16) and 2,705.5 m (A20) respectively. Bulk composition was assessed by X-Ray Diffraction (XRD) with combustion for TOC. Mercury injection porosimetry (MIP) was used to estimate porosity. The MIP tests were conducted on 1-1.5 cm 3 shale cubes cut from the core. The samples were dried in a humidity controlled oven to 90°C and were periodically weighed until their weight became stable. They were then analyzed in a Micromeritics Series V mercury injection porosimeter. The MIP tests were conducted in 40 pressure increments between 2 and 55,000 psi. Initial pressure increments were spaced at 3 psi and these gradually increased to 5,000 psi increments at the higher-pressure steps. The MIP data was conformance corrected and then the porosity (ϕ) was calculated using: where V Hg is the total volume of mercury injected into the sample and V b is the bulk volume of the sample determined during the MIP test. It should be noted that MIP could underestimate porosity due to the presence of disconnected pores that are not accessible to the mercury. Table 1 shows the bulk composition of the sample in terms of volume percent (vol%). Volume is preferred to weight because volumetric composition drives mechanical behavior (Herrmann et al., 2018;Rybacki et al., 2015;Sone & Zoback, 2013a), for example, by determining, together with grain arrangement, the existence of any load-bearing framework (Rybacki et al., 2016). In addition, in heterogeneous materials such as shales, volume fraction is also surface fraction  and so is directly related to the mineral phases expected to be encountered in grid nanoindentation. The values in Table 1 are converted from the measured weight percentages  (Okiongbo et al., 2005).
The data shows that the overburden samples are clay-rich (A3 is 63.1 vol% clays and A6 54.2%), with the clay fraction dominated by illite with lesser amounts of illite-smectite (having less than 20% smectite layers), chlorite and kaolinite. The Fort Simpson formation is poor in OM and the samples considered here have OM around 2 vol%. In wt% terms this is 1.3% (A3) and 0.9% (A6), values which are slightly higher than the general TOC of <0.5% reported by Ross and Bustin (2008). This may be due to the samples' proximity to the Muskwa contact. Silt particles are mostly quartz (∼30%), with small amounts of albite (plagioclase feldspar), pyrite, and carbonates (dolomite and siderite). Porosity measured by MIP is 1.7% in A3 and 1.4% in A6; this is slightly lower than previous MIP data reported by Ross and Bustin (2008), who gave a porosity range of 1.9%-4.65% across eight Fort Simpson samples.
The reservoir members have notably different compositions to the overburden, with much greater amounts of OM. The Muskwa sample (D1) is primarily quartz (85.6%) with a minor fraction (7.4%) of feldspar (albite, microcline) and an absence of clay minerals. OM is 4.3% in volume terms (2.2 wt%) and MIP-porosity is low (1.1%); both are quite similar to results reported by Ross and Bustin (2008). Sample A16 is from a carbonate-rich lithology of the Otter Park member, dominated by dolomite with minor quartz/feldspar and a small presence of muscovite. TOC is 2.3 wt% while the measured porosity is minimal (0.1%). As expected, the Evie sample (A20) has the highest organic content (5.4 wt%). It is a mixture of quartz (43.7%), clay (22.7%), and carbonate (13.5%) with a significant 10.1% OM in volume terms.

Texture
Sample texture was analyzed by scanning electron microscopy (SEM). Polished thin sections were carbon coated and imaged using an FEI Quanta 650 field emission SEM, which is fitted with an Oxford Instruments INCA 350 EDX system/80 mm X-Max SDD detector. Back-scattered electron mode was used with a working distance of 11 mm and 20 keV accelerating voltage.  Figure 2e) is quartz-rich, and the SEM imaging indicates that much of the quartz is likely to be diagenetic. Figure 2f shows a framework of authigenic quartz cement, resulting from recrystallized biogenic silica. Silt-size grains of dolomite are dispersed alongside smaller, predominantly euhedral, pyrite grains with OM mixed deeply into the quartz framework. Sample A16 is from a carbonate-rich lithofacies of the Otter Park member. The shale mainly consists of rhombs of diagenetic dolomite (Figures 2g and 2h) and approximately half of the dolomite has a high iron content. Pyrite is euhedral and quite sparse, with the remainder being scattered quartz grains and a clay/organic phase. Figures 2i and 2j show that sample A20, from the Evie member, has a more mixed mineralogy, and the microstructure is a collection of quartz, carbonate (calcite and dolomite), clays, OM and pyrite (euhedral and framboidal) grains. It is noticeable in Figure 2j that there is very little orientation of the clay minerals, especially when compared with the overburden samples. Alongside the relatively low clay content, this indicates that the sample is rigid-grain supported, which disrupts clay alignment during compaction and diagenesis (Day-Stirrat et al., 2010). There is some evidence of authigenic quartz cements, but this is much less apparent than in D1 and cannot be conclusively identified without cathodoluminescence data (Milliken et al., 2012). Figure 2j suggests dispersed fragments of micron-scale OM of indeterminate origin.

Nanoindentation
Nanoindentation testing was carried out using a NanoTest Vantage system to investigate the mechanical response of the shale at a microstructural level. Samples were cut both parallel and perpendicular to the bedding plane and the sections were prepared by broad ion beam polishing to minimize surface roughness. A diamond Berkovich indenter was pushed into the samples to a maximum load (F max ) of 5 mN, selected so that the indentation measures the distinct responses of the individual constituents. The load/unload rate was 0.5 mN/s and the load was held at F max for 60s to allow for system equilibrium and to investigate the creep response. The indentation was conducted on 15 × 15 grids, with indents spaced 5 μm apart given expected indentation depths of a few hundred nanometers and an affected volume 3-5 times this depth (Abedi, Slim, Hofmann, et al., 2016).
The testing protocol and typical responses are illustrated in Figure 3, where the indentation direction relative to the bedding plane is also defined. The indentation (reduced) modulus, E*, and hardness, H, were extracted from the force-displacement (F − h) curves ( Figure 3a) (Oliver & Pharr, 2004) where A c is the contact area ( ≃ 24.5ℎ 2 for a Berkovich indenter) and S = (dF/dh) is the initial unloading gradient.
The time-dependent response was characterized in terms of a creep modulus, C. This is defined by the contact creep compliance function (Vandamme & Ulm, 2013), which describes the creep response due to a unit input of stress and corrects for the geometry of the indenter. A logarithmic fit to the creep displacement is used: where t 0 is the time and h 0 is the displacement at the start of the holding stage and c f and τ are fitting coefficients which can be used to characterize the creep response. The suitability of the logarithmic fit is shown in Figure 3b, which shows a typical response obtained during testing. The coefficient τ represents the characteristic time at which long-term logarithmic creep behavior starts. Converting to creep compliance, the creep modulus is: where is the radius of the projected contact area between the indenter and sample, = √ ∕ . As for hardness (Equation 3), the contact area at the end of the hold period is used to determine .

Chemo-Mechanical Coupling
Interpretation of nanoindentation testing generally relies on statistical approaches due to the "big data" nature of massive grid indentation campaigns. Several recent studies have shown the benefit of coupling chemical testing with mechanical nanoindentation data to improve the identification of material phases in the shale microstructure (Abedi, Slim, Hofmann, et al., 2016;Deirieh et al., 2012;Y. Liu et al., 2022;Prakash et al., 2021;Veytskin et al., 2017). This is particularly important when the geomechanical properties are linked to upscaling schemes: a more robust characterization of the fundamental micromechanics should allow a better prediction of properties at larger scales.
To this end, SEM-EDS was conducted on a selection of nanoindentation grids (post-indentation), and a Gaussian mixture model (GMM) used to extract the mechanical properties of the material phases from the coupled chemo-mechanical data. Gaussian distributions have regularly been applied in the deconvolution of nanoindentation results (usually E* and H) on shale (e.g., Du et al., 2022;Li et al., 2019;Ulm et al., 2007). Here, a lognormal distribution is used in the mixture model for the creep modulus, which stems from the logarithmic fit in Equation 4. In this case, the natural logarithm of C is normally distributed: ln ∼  ( ) , where the parameters α and β are respectively the mean and variance of ln C. Therefore, ln C can easily be used in a standard GMM clustering algorithm. The mean ( ) and variance ( 2 ) of C are respectively: The GMM procedure implemented in Matlab (MATLAB, 2020) was applied, in which likelihood is optimized using the iterative expectation-maximization algorithm.
The methodology is illustrated in the flowchart in Figure 4 and can be described as follows: 1. Low load nanoindentation grids are conducted to obtain mechanical properties (E*, H, C).
2. The grids are located and imaged using SEM and chemical analysis is undertaken with EDS, guided by the mineral composition determined by XRD analysis. 3. The EDS element intensities at each indent are combined with the mechanical data and GMM clustering is used to identify the material phases.

Experimental Program
The experimental program was carried out in two sets ( Table 2). All nanoindentation testing was carried out according to the protocol described in Section 3.1. Chemo-mechanical coupling was carried out as described in Section 3.2. Set 1 was mechanical-only, and grids were conducted on polished sections of all samples in parallel and perpendicular indentation directions. In Set 2, additional sections of samples A3, A6, and A20 were prepared and extra nanoindentation grids conducted in both parallel and perpendicular directions, with SEM-EDS undertaken after indentation. The samples from the Fort Simpson and Evie shales were targeted for chemo-mechanical coupling due to their more varied mineralogy.

Nanoindentation Results (Set 1)
The mechanical results obtained from Set 1 of nanoindentation tests are shown in Figure 5. Indents for which the load-displacement behavior was strongly affected by surface damage, resulting from existing fractures or particle and projected on to the other plots as reference lines. In (a-ii) and (b-ii), the shaded area indicates a soft zone containing the likely properties of the clay/organic matrix. The mechanical properties indicated by the red circle on (c-ii) and blue circle on (d-ii) were measured by Yang et al. (2020) through nanoindentation testing of quartz and dolomite respectively. breakage and flaking, have been removed from the analysis, as have indents where the creep behavior did not conform to a log fit. The retention rate was at least 88% for each grid.
The load displacement curves (Figures 5a-i to 5e-i) demonstrate typical load-hold-unload stages and reveal that the overburden and reservoir samples show very different mechanical behavior. The maximum displacement is around 750 nm for overburden samples A3 and A6 in comparison with <500 nm for the reservoir samples, indicating that the reservoir samples give a generally harder response. In addition, the overburden shows a clear anisotropy, with a distinct response when indented parallel (maximum indentation depth 508 nm in A3 and 509 nm in A6) compared with perpendicular to the bedding plane (maximum indentation depth 675 nm in A3 and 730 nm in A6). The reservoir lithologies are more isotropic in terms of maximum indentation depth in the parallel and perpendicular directions, with only A20 showing a minor difference. The measured depths in parallel (resp. perpendicular) indentation directions were: 386 (389) nm in Muskwa D1; 395 (389) nm in Otter Park A16; and 435 (499) nm in Evie A20.
The calculated indentation modulus and hardness (which is a function of the maximum indentation depth) are plotted for each indent in Figures 5a-ii to 5e-ii. In general, a positive correlation is expected, with E* approximately proportional to √ (Ulm & Abousleiman, 2006). The nanoindentation data follows a positive correlation but the data tends to plot in groups, with the shale constituents having different mechanical signatures. Power laws were fitted to the quartz-rich D1 and dolomite-rich A16 samples to help identify the different clusters (Figures 5c-ii and 5d-ii); note that the power exponent is very close to 0.5 in both cases. The fitted power laws pass through the mechanical properties identified as characteristic of quartz (red circle on Figure 5c-ii) and dolomite (blue circle on Figure 5d-ii) by Yang et al. (2020), who used nanoindentation testing in a dynamic mode to obtain continuous depth-dependent measurements of E* and H on shale samples from the Lower Silurian Longmaxi formation in the Sichuan Basin, China.
In terms of bulk mineralogy, each sample possesses at least 20% quartz, and this can be recognized in Figures 5a-ii to 5e-ii as those indents which closely follow the D1-fitted power law. In a similar way, the nanoindentation grid on sample A6 locates a cluster of indents on dolomite. The range of hardness values is likely due to a substrate effect in which mineral grains interact with softer phases or the pressure-bearing framework at larger indentation depths, hence the general relationship of stiffness and hardness reducing toward a soft phase (Yang et al., 2020). Note that the D1-fitted power law tends to be slightly stiffer than the quartz phases in the other samples, which is likely to be a result of the diagenetically stiffened framework of authigenic quartz cement in sample D1. Using wireline logs, Dong et al. (2017) found that authigenic quartz was associated with "brittleness" and high stiffness in the bulk geomechanical response of core samples from reservoir lithologies in the Horn River Basin; in contrast, detrital quartz had little effect on the bulk mechanical properties due to the grains being situated in a load-bearing clay matrix.
In the sample set considered here, OM typically appears as small particles, often intimately mixed into the shale microstructure (see Figure 2) and it is not possible to resolve such fine particles with nanoindentation testing; instead, we consider a fine-grained phase of porous clays and OM. Indents in the softer region (E* < 60 GPa and H < 2.5 GPa) are likely to belong to this matrix of clay and OM, and this phase is particularly evident in the clay-rich overburden samples A3 and A6 (see shaded areas in Figures 5a-ii and 5b-ii). These clusters show an obvious anisotropy with samples indented perpendicular to the bedding plane showing a softer response. This phase is crucial as the clay matrix has been found to be the most significant factor in determining brittleness in the Horn River reservoir units, with high clay content associated with ductility (low brittleness and low stiffness) (Dong et al., 2017). Clay-rich facies of the Otter Park shale have been observed to act as fracture barriers (Yu & Shapiro, 2014).
The creep modulus is also expected to be positively correlated with H (Vandamme & Ulm, 2013), and such a relationship can be observed in the data in Figures 5a-iii to 5e-iii. The creep modulus is plotted on a logarithmic scale to aid comparison. For sample A3, the characteristic time τ is typically ≤ 1s and C = 10 2 , 10 3 , and 10 4 GPa correspond to respective creep displacements of approximately 60, 15, and 5 nm over the 60s hold period. It is notable that in samples A3 and A6, there is a cluster of indents with low creep modulus, approaching 100-200 GPa, which most likely represents the soft (clay/organic) phase, while the harder phases possess higher C. Reservoir samples show generally higher values of creep modulus. While there is little evidence of anisotropy in the reservoir samples in terms of E* and H, both D1 and A16 show an anisotropic creep response (Figures 5c-iii and 5d-iii). For D1, the creep modulus is clearly higher when indented in the parallel direction compared to the perpendicular direction, for which the creep modulus is surprisingly low given the quartz-dominated composition. This is discussed in Section 4.3. In contrast, the anisotropy is reversed in sample A16, with C tending to be higher in the perpendicular direction. Figure 6 shows SEM-EDS analysis of the nanoindentation grids on samples of the Fort Simpson overburden (A3 and A6) in both parallel and perpendicular directions. The chemical analysis reveals the distribution of material phases targeted by the nanoindentation grids: the matrix phyllosilicates are rich in Al and Si while quartz grains (high Si) are frequent inclusions, approximately 5 μm in size. The grids also fall on isolated grains of dolomite (high Ca and Mg), which are larger (∼10-20 μm diameter) and more frequent in sample A6. In both samples, framboidal pyrite (high Fe) and albite (high Na) are scattered in the matrix as small particles, a few microns across. The clay minerals are strongly aligned, and the bedding direction can be clearly observed in Figures 6a  and 6c. The SEM images show that in the perpendicular direction, the indentation led to flaking of clay particles, particularly around the harder mineral inclusions; this is indicated on Figure 6d where clay minerals have broken from the sample surface around large dolomite grains. This particle breakage tends to show as "pop-ins" in the indentation load-displacement curves, or as large initial displacements if the surface is damaged, and those indents have been removed from further analysis to avoid influencing the results (the retention rate was at least 84% for Set 2).

Overburden (Fort Simpson)
The nanoindentation data was linked to EDS chemical analysis to identify the mechanical response of the different material phases in the shale, following the procedure described in Section 3.2. The nanoindentation grids on A6 were positioned on areas with varied mineralogy, and so this sample is used to illustrate the analysis procedure. In the parallel indentation direction, three material phases were assumed: a clay matrix, a quartz/feldspar phase, and dolomite inclusions. Figure 7 shows the results of the clustering approach using only mechanical data and using coupled chemo-mechanical data with Al, Si, and Ca intensity values. Using only mechanical data, shown in Figures 7a and 7b, a soft phase is identified, which is likely to be the clay matrix, and two stiffer clusters are identified, separated in terms of hardness and creep modulus. We note that the plot of ln C -H conforms well to the Gaussian error ellipses, indicating that a lognormal distribution is suitable for C. Chemo-mechanical clustering is shown in Figures 7c and 7d with the normalized intensity values for each phase shown in Figure 8. This approach identifies a very similar clay matrix (rich in Al and Si) to the mechanical-only clustering but finds quartz (Si > Al) and dolomite (high Ca) phases with a strong positive correlation between E* and H. This is evidence of a "substrate effect" which agrees with the interpretation of the Set 1 data (Section 4.1). The values of creep modulus for the dolomite and quartz phases are mixed, with a less clear distinction between these phases in terms of C than E*.
The clusters are plotted on a composite SEM-EDS map in Figure 9 and are successfully located almost directly on the corresponding particles. The GMM volume fractions (67% clay matrix, 19% quartz + feldspar, and 15% dolomite) are comparable to the bulk XRD but with a greater proportion of dolomite due to the grid being positioned on large dolomite grains. The group of small pyrite framboids at the bottom of the grid causes some difficulty as the indents tend to fall at the grain edges; these indents tend to be grouped into the dolomite phase.
The chemo-mechanical clustering of the indentation grid on sample A6 perpendicular is shown in Figure 10. In this direction, the clay matrix is shown to have mean mechanical properties (± one standard deviation) of * 3, = 22.6 ± 4.1, H 3,cly = 1.1 ± 0.5, and C 3,cly = 424 ± 350 GPa. In contrast, in the parallel direction, the mean values are * 1, = 39.7 ± 8.7, H 1,cly = 1.5 ± 0.5, and C 1,cly = 458 ± 231 GPa. This reveals the intrinsic anisotropy of the clay matrix, which shows a far softer mechanical response when indented perpendicular to the bedding plane. This exacerbates the "substrate effect" and to account for this, an additional "mixed" phase was included in the clustering procedure. As shown in Figure 10, the mixed phase sits between the mechanical properties of the clay matrix and those of the quartz and dolomite phases. Figure 11 plots the GMM clusters on the composite SEM-EDS map and the mixed phase often lies at grain boundaries, for example, along the edge of the dolomite rhomboid at the top of the grid. In contrast, the GMM-clustered dolomite and quartz indents are mostly located on large particles of the corresponding minerals; by isolating these phases, the mechanical properties are shown to be comparable to those in the parallel direction which implies a relatively isotropic behavior in comparison with the strongly anisotropic matrix.
Sample A3 is primarily a clay matrix with quartz/feldspar inclusions. Dolomite was undetected using XRD but the EDS analysis indicates that the grid in the parallel indentation direction was located on a small dolomite grain, which can be revealed in the chemo-mechanical clustering (Figures 12a and 12b) and the SEM-EDS map (Figure 13a). The identified chemo-mechanical clusters are very similar to sample A6. In the perpendicular indentation direction (Figures 12c and 12d), the clay matrix again gives a much softer response and a "mixed" phase is identified to separate the matrix and inclusions. Figure 13b indicates that the clustered quartz/feldspar indents are located on intact grains; this phase shows a lower hardness compared to the parallel direction, which is most likely due to the softer substrate. Figure 14 presents SEM-EDS analysis of the Set 2 indentation grids on sample A20. The mineralogy is more varied than the Fort Simpson samples, particularly in terms of the carbonate content. The EDS maps show small grains of calcite (generally <5 μm across) mixed within the clay matrix, with less

10.1029/2022JB025957
14 of 24 frequent dolomite particles. Quartz and albite grains dominate the inclusions alongside euhedral pyrite. In the parallel-indented sample, some orientation of the clay particles is evident although this is not as clear as in the Fort Simpson shales. Flaking of the clay phase can also be observed in the perpendicular-indented sample, but the damage is less extensive than in the overburden samples.
The chemo-mechanical clustering of the A20 indentation data is shown in Figure 15. In the parallel indentation direction, one indent (indicated on Figure 15a) can be easily identified as pyrite, having stiffness and hardness much greater than other phases and corresponding well with previous estimates, for example, E py = 265 GPa (Whitaker et al., 2010). This data point was excluded from the GMM algorithm. The clustering reveals a similar pattern to the Fort Simpson samples, with the addition of a calcite phase that is mixed closely with the clay matrix but generally distinguished by a higher stiffness and lower hardness. In the parallel indentation direction, the mean mechanical properties (± one standard deviation) for calcite are * 1, = 48.9 ± 11.5, H 1,cal = 2.0 ± 0.6, and C 1,cal = 631 ± 292 GPa, and for the clay matrix are * 1, = 42.5 ± 8.9, H 1,cly = 2.2 ± 0.7,  and C 1,cly = 702 ± 298 GPa. Comparing Figures 15a and 15c, the clay matrix becomes much softer in the perpendicular indentation direction whereas the calcite phase shows less anisotropy (mean * 3, = 42.0 ± 11.6, H 3,cal = 2.1 ± 0.8, and C 3,cal = 882 ± 546 GPa; * 3, = 31.0 ± 5.6, H 3,cly = 1.8 ± 0.6, and C 3,cly = 875 ± 532 GPa). This result agrees with a recent micromechanical investigation of several calcareous shale samples (Graham et al., 2022), where fine (clay-sized) calcite grains were shown to possess anisotropic elastic properties, being stiffer when indented parallel to the bedding plane, but with a less pronounced anisotropy than observed for clay minerals. However, the anisotropy of the Evie clay phase is also lower than observed in the matrix of the Fort Simpson samples. Figure 16 shows the GMM-clustered indents on composite SEM-EDS maps. Again, the clustering algorithm is quite successful in identifying the different phases in the shale composite. The Evie sample (A20) is generally finer-grained than the Fort Simpson samples analyzed previously. Underneath the grids, the calcite is mixed into the clay matrix and appears in small grains, generally less than 5 μm across. Similarly, pyrite (both euhedral and framboidal) tends to also appear as small particles. The fine particle size is particularly evident in the perpendicular indentation direction shown in   Figure 16b, where an Fe/Mg-rich phase is identified that includes both dolomite and pyrite. The small grain size relative to the indent depth (∼500 nm) and grid spacing makes it challenging to distinguish between the different phases as some indents will inevitably be located at grain boundaries or impact on several different phases, leading to a degree of homogenization in the measured responses.

Assessment of Elasticity and Creep
The mixture models generated from the coupled chemo-mechanical data (Set 2) can also be used to extract the mechanical properties of the different phases from the Set 1 data. This is achieved by first carefully selecting the components of the GMM based on the indentation direction, sample composition (Table 1), and a visual interpretation of the results (such as described in Section 4.1). Then, each Set 1 data point is assigned to the GMM component with the highest posterior probability for that data point, in a process known as hard clustering (Murphy, 2012). As an example, Figure 17 shows the results of sample A6 (Set 1) in the parallel indentation direction, clustered using the chemo-mechanical GMM fitted in Set 2; the figure demonstrates that the results of Set 1 and Set 2 are very similar (compare Figure 7) and the GMM is able to satisfactorily distinguish between the material phases. Phase statistics of Set 1 can then be calculated from the clustered indents (rather than relying on the individual Gaussian distributions, which were fitted to the Set 2 data).
The stiffness and creep properties of the dominant material phases (clay matrix, quartz/feldspar, dolomite, and calcite) across the reservoir and overburden units are shown in Figure 18, incorporating both Set 1 and Set 2 data. The Fort Simpson shale samples consist of a soft, highly anisotropic clay matrix (mean * 1 ∕ * 3 = 1.8) with comparatively homogeneous inclusions of quartz/feldspar and dolomite. The significant seismic anisotropy measured in the field (Baird et al., 2017;Yu & Shapiro, 2014) has been linked to the orientation of the clay fabric, and the micromechanical results confirm that the stiffness anisotropy of the clay phase is a dominant factor in the overall anisotropy of the Fort Simpson shale. A highly anisotropic clay fabric is also consistent with previous studies showing how clay minerals in mudstones become increasingly aligned during burial due to diagenetic recrystallization, with the degree of alignment moderated by the bulk clay:silt ratio (e.g., Day-Stirrat et al., 2017Ho et al., 1999). Creep is also concentrated in the matrix, where the lowest creep modulus is found, indicating that this is the key phase for potential self-healing of fractures in shales; this is important in the context of both shales acting as seals and as hydraulically fractured reservoirs. This is demonstrated in Figure 19, where the average creep displacement is shown for both the clay matrix and quartz/feldspar phase, and it can be seen that the creep displacement is much greater in the soft matrix. The time-dependent response of the matrix also shows some anisotropy, particularly in sample A3 where C 1 /C 3 = 1.85, compared to C 1 /C 3 = 1.10 for sample A6. The creep modulus of the inclusions (quartz/feldspar and dolomite) is much higher than the matrix values with little evidence of anisotropy, although C qtz tends to be more variable due to the greater range in hardness of the quartz/feldspar phase compared to dolomite (see e.g., Figure 17b).  8 mN) nanoindentation grids on both immature (Antrim, Barnett, and Woodford) and mature (Marcellus and Haynesville) samples of U.S. shale gas formations and extracted the mechanical properties of the matrix by clustering. The figure shows that the Fort Simpson matrix is broadly comparable with that of the mature U.S. shales but shows a greater anisotropy both in terms of elasticity and creep. Indeed, the U.S. shales show no clear creep anisotropy. This may be due to the influence of OM, which is generally considered to have an isotropic geomechanical response; the Fort Simpson formation is organic-lean so that creep behavior is likely to be controlled by porous clay, whereas the samples tested by Slim et al. (2019) mostly had TOC higher than 2.5 wt%. The Evie sample is also organic-rich, but the clay/organic phase shows a lack of creep (high creep modulus) compared to the other shales. This may be due to microstructural effects, such as quartz grains acting as a rigid framework to limit creep displacement, or the influence of fine calcite grains mixed in the clay matrix.
The reservoir units have contrasting mineralogies. Sample D1 (Muskwa) is predominantly quartz and this shows an essentially isotropic stiffness that is the highest amongst all samples for the quartz/feldspar phase; this can be linked to the quartz-cemented nature of the microstructure in D1, in contrast to the small grains of detrital quartz that appear in other samples. The average creep displacement-time curves from the nanoindentation grids are plotted in Figure 19b, which also shows the creep displacement from the quartz/feldspar phase of the Fort Simpson samples. The figure indicates that the apparent anisotropy in the creep modulus of sample D1, where the creep modulus in the perpendicular indentation direction is close to the C values of the clay matrix, is slightly misleading. The magnitude of the quartz/feldspar creep displacement in the 60s hold period across these samples is actually quite similar, between 10 and 20 nm. The creep displacement in the perpendicular indentation direction of D1 tends to increase steeply over the hold period, without showing a decreasing strain rate, such that fitting the logarithmic function (Equation 4) gives an average characteristic time τ of 11.3s (compared to 3.1s for parallel indentation) and a low creep modulus. A longer hold period is needed to assess whether the creep displacement would continue to increase.
A16 is a carbonate-rich sample from the Otter Park unit and the micromechanical characterization reveals a stiff dolomite phase, which shows relatively minor creep, and little evidence of significant anisotropy. As discussed in Section 4.2.2, sample A20 (Evie) is a fine-grained mix of clays, quartz, feldspar, and carbonates, and, due to the resolution of the nanoindentation testing, it is a challenge to distinguish the various phases. The small grain sizes contribute to relatively low stiffness values for quartz and dolomite, particularly when indented perpendicular to bedding, and the lack of data points, which leads to a high uncertainty in the mean value, likely contributes to any apparent anisotropy. It is notable that in the parallel indentation direction, the clay matrix stiffness is comparable with the values measured in the Fort Simpson samples (around 40 GPa), but the anisotropy is much lower, with * 1 ∕ * 3 = 1.3. Figure 20 indicates that this level of stiffness anisotropy is quite typical for a porous organic-clay phase, but the creep modulus is substantially larger than may be expected, perhaps as a result of the relatively low clay content of the Evie sample.

Conclusions
This paper has reported the results of a micromechanical characterization of overburden and reservoir shales from the Horn River Basin. The shale samples were first described in terms of mineralogy, porosity, and texture. Low-load nanoindentation testing was then used to obtain the mechanical response of the shale constituents. A selection of nanoindentation grids was combined with SEM-EDS analysis and the coupled chemo-mechanical data used to identify the mechanical properties of the different material phases present in the shale. At the resolution of the nanoindentation grids, which covered an area of 70 × 70 μm with indents spaced at 5 μm, the major material phases were recognized as a clay/organic mixture alongside quartz/feldspar, dolomite, and calcite; these  Where two grids were conducted, the results of Set 1 are plotted to the left (circular markers) and Set 2 to the right (triangular markers). Solid markers show bedding plane-parallel indentation and empty markers show bedding plane-perpendicular indentation. Error bars indicate ± one standard deviation. appeared in different proportions and in different microstructural arrangements across the sample set. The main conclusions are as follows.
• The Fort Simpson shale forms the overburden, acting as the upper fracture barrier during hydraulic fracturing operations. SEM imaging showed that this shale consists of a strongly aligned clay matrix with granular inclusions of detrital quartz and some carbonates. The micromechanical results reveal that the clay/organic matrix is soft (low stiffness, hardness, and creep modulus) and highly anisotropic, particularly in terms of elastic stiffness, while the granular inclusions are much stiffer and show relatively isotropic properties; the clay matrix is therefore the key phase for fracture barrier or sealing applications. The elastic anisotropy of the clay fabric ( * 1 ∕ * 3 ) averaged 1.8 and this response can be linked to the very strong seismic anisotropy observed in microseismic data recorded in the field. • The reservoir units have varied mineralogical compositions and quite different microstructures to the overburden and in general gave a harder and more isotropic mechanical response, with less creep. This also agrees with microseismic data, where reservoir lithologies show less significant anisotropy and higher stiffness than the Fort Simpson shales. The microstructure of the Muskwa sample consisted of a framework of authigenic quartz cement, and this gave a stiffer response than that measured on the detrital quartz grains present in the overburden. The Evie sample was a diverse collection of clay and OM together with relatively fine grains of quartz, feldspar, carbonates, and pyrite. Fine carbonates were mixed into the clay/organic phase and it was notable that the clay showed a lower elastic anisotropy ( * 1 ∕ * 3 = 1.3 ) than the Fort Simpson samples, which may be related to a more random particle orientation.
• The different material phases showed unique mechanical signatures in the nanoindentation results across all samples, despite their different compositions and textures and different geomechanical roles. These signatures can be used to assist in the interpretation of nanoindentation results on shale, in particular to guide the use of automated clustering procedures incorporating only mechanical data, which could identify misleading clusters.
The micromechanical results allow a fundamental insight into the distinct geomechanical behavior of both the ductile overburden and brittle reservoir lithologies. Being able to identify and predict the link between shale composition and mechanical behavior in the field is relevant across many subsurface  energy and decarbonization technologies. Future work should incorporate the micromechanical data into homogenization schemes to upscale the mechanical properties from grain to core-scale.

Data Availability Statement
The nanoindentation and SEM-EDS data used in the study are available at the National Geoscience Data Centre (NGDC) via https://doi.org/10.5285/5c9eb939-ebc3-46f6-8690-ed1d9c46bca6 with an Open Government Licence .