Effect of Diagenesis on Geomechanical Properties of Organic‐Rich Calcareous Shale: A Multiscale Investigation

This study investigates the nano to core‐scale geomechanical properties of a maturity series of organic‐rich, calcareous shales buried to 100°C–180°C, with a focus on: (a) the mechanical properties of organic matter; (b) the elastic response and anisotropy of the shale composite at micro and core scale; and (c) the creep response. Atomic force microscopy was used to target kerogen at nanoscale resolution, and it was found that the elastic stiffness increased with thermal maturity from 5.8 GPa in an immature sample to 11.3 GPa in a mature sample. Nanoindentation testing of the shale matrix showed that diagenesis is a key factor in determining the bulk elasticity, with increasingly intense carbonate cementation at higher thermal maturities contributing to a stiffer response. A multiscale model was formulated to upscale the elastic properties from nanoscale solid clay minerals to a microcracked composite at core scale, with good predictions of the micro and core‐scale stiffness in comparison to indentation and triaxial results. A negative correlation was found between the creep modulus and clay/kerogen content, with greater creep displacement observed in nanoindentation tests in the immature clay‐ and kerogen‐rich sample compared to samples of higher thermal maturity.


Introduction
The robust geomechanical characterization of shale underpins the understanding and prediction of many processes of importance to both geology and engineering. These include the development and nature of folds, faults and fracture systems as a result of tectonic deformation (Albertz & Lingrey, 2012;Hesse et al., 2010;Mitra, 2003;Obradors-Prats et al., 2017) and the flow and leakage of fluids in the context of CO 2 sequestration, the storage of both hydrogen and nuclear waste and unconventional hydrocarbon production (Allen et al., 2020;Dewhurst & Hennig, 2003;Gale et al., 2014Gale et al., , 2007Holt et al., 2020;Nygård et al., 2006;Imber et al., 2014). In addition, the time-dependent mechanical behavior of shale is also very poorly constrained but has important implications for both fluid flow and wellbore stability on engineering timescales (Herrmann et al., 2020;Rybacki et al., 2017;Sone & Zoback, 2014).
"Shale" is the term used here to describe fine-grained sedimentary rocks which are effectively micro-composite materials comprising particles with diameters ranging from a few nanometers to 100 microns, and with highly variable mineralogical compositions including phyllosilicates, quartz, feldspar, carbonates, organic matter, and pyrite. The mineralogical variability and multiscale heterogeneity of shale means that predicting its mechanical properties is inherently difficult. Although correlations between mineralogy and the macroscopic mechanical response have been observed in some experimental studies, the relationship is complex and dependent on additional factors such as microstructure and the grain-scale diagenetic processes which cement, lithify and strengthen the granular mix without necessarily changing the bulk mineralogy (Herrmann et al., 2018(Herrmann et al., , 2020Sone & Zoback, 2013a, 2013b. Good quality core samples of shale suitable for conventional rock mechanics experiments are difficult and expensive to retrieve (Kumar et al., 2012). To fill gaps in core-scale data, it would be very useful to be able to predict the macroscopic mechanical response of shale from the intrinsic mechanical properties of its constituent minerals, and the interactions between mineral grains. This approach requires robust and high-resolution data of the mechanical response of shale at scales which represent the granular building blocks. A key challenge is thus to characterize the mechanical properties of sub-micron to 100-micron size particles, and the way in which the particles interact to generate the macroscale behavior of the composite material.
The mechanical behavior of organic matter is particularly difficult to constrain and is the focus of much current research. Due to its relative softness, the mechanical behavior of shale is greatly influenced by even a few percent of this constituent (Kumar et al., 2015;Sayers, 2013;Vernik & Milovac, 2011). Various methods have been used in attempts to constrain its mechanical properties, including nanoindentation (Ahmadov et al., 2009;Zargari et al., 2013;Zeszotarski et al., 2004) and ultrasonic measurements on samples of isolated kerogen (Yan & Han, 2013). However, the small size of patches of organic matter in the shale matrix, and its internal heterogeneity, requires resolutions below that of nanoindentation to measure in-situ properties (Eliyahu et al., 2015). Recent studies have focused on the use of atomic force microscopy (AFM), and in particular the PeakForce Quantitative Nanomechanical Mapping (QNM) mode, to map the mechanical properties of organic matter at nanoscale (Eliyahu et al., 2015;Emmanuel et al., 2016;Goodarzi et al., 2017;Li et al., 2018). Thermal maturation also affects the geomechanics of organic matter, and combined geochemical/geomechanical studies (Khatibi et al., 2018;Jakob et al., 2019) and molecular simulations (Bousige et al., 2016;Kashinath et al., 2020;Zhang & Jamili, 2015) are just beginning to reveal the detailed link between chemical composition, thermal maturity, and mechanical response.
The challenges of characterizing the mechanical response of shale in the laboratory, at a range of scales, carry forward into upscaling schemes. In principle, knowledge of the properties of shale constituents allow poromechanical models to be constructed to predict macroscopic behavior, for example based on multi-level homogenization techniques that use mineralogical characterizations to determine the contribution of various phases at several scales (Abou-Chakra Gué ry et al., 2010;Hornby et al., 1994;Ortega et al., 2007;Shen et al., 2013;Vasin et al., 2013). However, the complexities of diagenesis, microstructural interactions, and variability in the shale composite lead to difficulties in identifying the governing physics when attempting to link the mechanical response from the material building blocks to core-scale results relevant to geomechanical modeling of geological units. Even in complex multiscale frameworks, incorporating mechanical aspects ranging from weakened interfaces around individual silt particles to millimeter-scale microcracks, a high level of uncertainty exists at each level which may be best assessed by a probabilistic approach (Dubey et al., 2019).
A further difficulty in characterizing the geomechanical behavior of shale is its time-dependent nature (Sone & Zoback, 2014). Knowledge of the creep response of shale is important to determine how hydraulically-induced fractures behave over engineering timescales, with fracture closure having implications for production rates during hydrocarbon extraction and for fracture healing in underground storage of CO 2 or H 2 (Cerasi et al., 2017;Rybacki et al., 2017;Sone & Zoback, 2013b). The time-dependent response of shales is also relevant over a range of scales and researchers have begun to investigate the multiscale nature of creep in shales, for example by nanoindentation tests with hold periods (Mighani et al., 2019;Slim, 2016;Slim et al., 2019).
The objective of this study is to use multiscale testing to constrain the mechanical response of shale at a microstructural level before comparing with core-scale results. It focuses on: (a) the mechanical properties of organic matter at nanoscale; (b) the elastic response and anisotropy of the shale composite, and how microscale measurements relate to triaxial test results; and (c) the creep response of shale at different scales. A maturity series of organic-rich shales buried to 100°C-180°C (immature to overmature) is used here to investigate the link between diagenesis, microstructure, and geomechanical behavior. The experimental methods include PeakForce QNM and high load nanoindentation testing and several practical issues in the application of these techniques to shale are discussed.

Sample Description
For this study, a set of samples from the Posidonia shale was selected with different maturity levels. The Posidonia shale formation from the Lower Saxony Basin is a fine-grained calcareous shale that contains high amounts of organic matter (Type II kerogen). In this study, we use the term "kerogen" as a catch-all for different forms of insoluble solid organic matter, encompassing both primary (derived from initial biogenic precursors) and secondary material (resulting from thermal alteration of primary organics) (Vandenbroucke & Largeau, 2007); in many practical contexts the distinction is unimportant as downhole logs cannot distinguish between them (Craddock et al., 2020). This formation is among the earliest hydrocarbon-producing basins in Europe and has been considered one of the most important source rocks in Western Europe (Rexer et al., 2014). Here, three samples were selected from the Wickensen (WIC), Harderode (HAR), and Haddessen (HAD) boreholes located in the Hils half-graben, northern Germany. For the geomechanical experiments, thin sections of 200 μm thickness were prepared both parallel and perpendicular to the bedding plane and ion-beam polished to ensure a flat surface. The moisture content of the samples was not controlled.

Composition
The mineralogical composition was determined by X-ray powder diffraction (XRD) and is summarized in Table 1. Total organic carbon (TOC) was measured with a LECO C-S analyzer. The bulk density,  b, was obtained from immersion in mercury at 0.172 MPa and grain density, rho g , was measured by the Small Pycnometer Method (Rexer et al., 2013). Sample porosity is then ϕ = 1 − (ρ b /ρ g ). Typical mineral densities ) were used to convert from weight to volume fractions. The kerogen volume fraction was estimated through the empirical relationship (Vernik & Milovac, 2011): where f k is the volumetric kerogen content, TOC is the weight % of organic carbon, ρ k is the kerogen density and C k is an empirical coefficient relating kerogen to measured TOC. The factor C k accounts for the fact that kerogen consists not only of carbon but also contains oxygen, nitrogen, and sulfate, etc. Realistic values range between 0.6 and 1.0, with C k increasing with thermal maturity from around 0.7 for immature shales toward 1.0 in the gas window (Rudnicki, 2016). The three samples have vitrinite reflectance (R o ) values of 0.53% (WIC), 0.89% (HAR), and 1.45% (HAD) and are considered immature, mature, and overmature respectively (Mathia et al., 2016). Previous measurements of kerogen density in Posidonia shale by Rexer et al. (2014) give average values of ρ k = 1.23 g/cm 3 for WIC and 1.35 g/cm 3 for HAD, which agrees with data from Kimmeridge Clay (Type II kerogen) at similar temperatures reported by Okiongbo et al. (2005). A value of 1.25 g/cm 3 was assumed for HAR because, according to chemical structure modeling (Walters et al., 2007) and back analysis from weight % TOC and grain density measurements (Rudnicki, 2016), only a slight increase in kerogen density is expected as % o R approaches $\sim$ 0.9 before a step change in density occurs through the remainder of the oil window.

Microstructure and Diagenetic Alterations
The samples show substantial maturity-related changes in texture at the microscale. To reveal these alterations and allow a link to be made with the geomechanical response, polished thin sections were carbon coated and imaged using a Hitachi SU-70 High Resolution Analytical Scanning Electron Microscope, equipped with an Oxford Instrument Energy Dispersive X-ray microanalysis system (INCA Energy 700   Figure 1 shows representative BSE images. Further discussion on the diagenetic evolution of the Posidonia maturity series can be found in detailed studies by Mathia et al. (2016 and. The immature WIC-7129 sample is a microlaminated calcareous mudstone with bulk mineralogy (Table 1) indicating similar proportions of calcite and clay. Figure 1a shows the carbonate fraction is mostly of biogenic origin (fecal pellets). Nannofossils, generally broken coccoliths of size ≤5 μm, are dispersed in the matrix. The sample is organic-rich (11.8 weight % TOC) and the organic matter is most prominent in the BSE image as thick bodies of algal liptinite. The alginite macerals are generally associated with clay-rich areas, forming a porous clay/kerogen matrix which also contains amorphous bituminite (Mathia et al., 2016). Siltsized quartz grains and pyrite framboids are scattered throughout the matrix and there is limited evidence of diagenetic alterations to the microstructure.
The mature HAR-7046 sample is of similar composition to the immature shale. Less organic matter is present (6.0 weight % TOC) and Figure 1b indicates that alginite lenses have collapsed due to mechanical compaction, appearing much finer in comparison to the voluminous bodies observed in WIC-7129. Thermal maturation has led to a network of non-solvent extractable solid bitumen in the matrix, filling pore space in clay and fossil-rich areas (Bernard et al., 2012;Mathia et al., 2016). The presence of small patches of calcite cement and partial recrystallization of the nannofossil-rich matrix are evidence of diagenesis. The diagenetic processes mainly involve precipitation of calcite cement between grains, and in places filling algal bodies, alongside local replacement of the original fabric through dissolution and reprecipitation of carbonates.
Carbonate diagenesis becomes more intense into the gas window. This is demonstrated in Figure 1c, which shows that the microstructure of the overmature HAD-7127 sample bears little resemblance to the less CHARLTON ET AL. mature samples. Nannofossil structures are lost and the shale matrix consists almost entirely of authigenic dolomite and calcite cement. The diagenetic transformation has obliterated any lamination. Dolomite predominates and the bulk mineralogy shows a substantially higher dolomite content than WIC-7129 or HAR-7046. The dolomite content is also higher than typical HAD samples (Mathia et al., 2016), suggesting the local availability of a source of magnesium. There is much less organic matter (2.3 weight % TOC) in comparison to the samples of lower maturity and its appearance is not easy to identify from the BSE images, likely trapped in small pockets in the cemented matrix. Previous studies have found the organic matter in HAD samples dominated by a tight network of solid bitumen in the clay-carbonate matrix with no alginites remaining (Bernard et al., 2012;Mathia et al., 2016).

Atomic Force Microscopy
PeakForce QNM is an advanced mode of the conventional AFM test, which enables the production of both topographical and mechanical maps of a sample surface (Bruker, 2012). In this non-destructive test, a probe is tapped over the surface of the sample and the induced peak force in the probe, which is kept constant during scanning, is continuously monitored to track the surface. For every approach and withdrawal of the tip to the surface, a force-separation curve is generated from which the reduced modulus, E*, of the surface can be obtained from the unloading (withdrawal) curve (Trtik et al., 2012). The standard Derjaguin-Muller-Toporov (DMT) contact model (Derjaguin et al., 1975) was used as the basis for the calculation of reduced modulus: where F is the applied force, adh F is the adhesive force during contact, r tip is the tip radius, and h is the indentation depth.
The reduced modulus is related to the sample Young's modulus, E, through: where ν is the Poisson's ratio of the sample, and E tip and ν tip are respectively the Young's modulus and Poisson's ratio of the tip.
The resolution of data collection depends on r tip , which is between 10 and 40 nm. In addition, the reliable range of measurement for each probe is correlated with the tip stiffness: the stiffer the tip, the higher the range of measurement. Here, a diamond tip (DNISP) was selected which can be reasonably assumed as rigid due to its stiffness (1,100 GPa) in comparison with that of the shale constituents. The standard HOPG-12M was used for calibration (Trtik et al., 2012), with values of 1 kHz and 50-150 nN chosen for oscillation frequency and peak force respectively.
The focus of the AFM measurement was the organic matter, which mainly appears as elongated alginite in the clay matrix of the immature WIC-7129 sample and collapsed algal lenses and pore-infilling solid bitumen in the mature HAR-7046 ( Figure 1). Sections parallel to the bedding plane were investigated and statistical analysis of the nanomechanical maps was used to estimate the reduced modulus of kerogen, denoted * k E . In addition, fictitious low stiffness values produced as a result of existing holes and cracks on the surface were carefully investigated in order to avoid any influence on the estimated properties. Figure 2 shows the map of reduced modulus for two target areas of WIC-7129 alongside the histogram of modulus measurements. A kernel density estimate was used to identify the peaks in the distribution, which correspond to different phases in the shale. It can be observed that both target areas show a clear peak at low values of elastic modulus, which can be attributed to kerogen. In both target areas, the peak corresponding to the organic material occurs at around * k E = 5.8 GPa, which confirms the consistency of the measurements.

Immature Sample (WIC-7129)
In order to increase the level of confidence that the measured values belong to the kerogen, energy-dispersive X-ray spectroscopy (EDS) analysis was performed. Figure 3 shows a BSE image and EDS chemical analysis of the region around the nanomechanical target area shown in Figure 2b. The target area is shown as a blue rectangle on the BSE image. Large areas rich in carbon were detected and agree well with both the dark lens on the BSE image, which corresponds to organic carbon, and the low value of elastic modulus on the nanomechanical map shown in Figure 2b.
Several surface fractures or cracks are visible on Figure 3, which can also appear on the nanomechanical map, showing low stiffness values. However, these are easily distinguished from organic matter by their high aspect ratios and morphology, as well as by EDS mapping, and would have little influence on the location of histogram peaks. The potential issue of sharp concavities, which may resemble in shape true areas of organic matter, was also considered. Figure 4 shows a topographical map of the AFM target area presented in Figure 2a. A cavity can be observed, which may have been caused by silt inclusions being pulled out during the polishing procedure. The reduced modulus along two cross sections over the cavity are shown in Figure 5. The measured values over the cavity are less than 2 GPa, indicating that the peak found at 5.8 GPa belongs to kerogen in the shale microstructure and is not an artifact of topography.

Mature Sample (HAR-7046)
In order to observe the effect of thermal maturity on the mechanical response of organic matter, AFM was also conducted on the mature HAR-7046 sample. Figure 6 shows a nanomechanical map of a target area CHARLTON ET AL.
10.1029/2020JB021365 6 of 23 alongside the histogram of the reduced modulus. Again, a peak for a soft phase (around 9 GPa) can be seen, which could relate to kerogen. However, the peak is not as prominent as observed in the immature sample, which is due to the organic matter being present in much smaller patches.
To focus on a smaller area of the mature sample, the diamond tip was replaced with a silicon nitride (TAP525) tip. The radius of the latter is reported at 8 nm (we measured 13.5 nm by probing the tip evaluation sample, made of polycrystalline titanium), which is substantially finer than the 40 nm for the diamond tip, allowing an improved resolution. The silicon nitride tip is less stiff (the spring constant for the cantilever was measured at 80 Nm −1 compared to 272 Nm −1 for the diamond tip) but the accurate range is 1-30 GPa, which is sufficient for the organic phase.
The silicon nitride tip was used to scan an area with dimensions less than 5 μm. Figure 7 shows the results for this test and a very clear peak can be identified for the organic matter. The reduced elastic modulus corresponding to the peak value is at * k E = 11.3 GPa, indicating that increased thermal maturity has led to a stiffening of the kerogen. The nanomechanical map demonstrates that it is difficult to distinguish between the clay matrix and silt inclusions due to the limited range of accurate measurement of the silicon nitride tip.

Assessment of Elastic Properties of Organic Matter
Due to its nanoscale resolution, AFM has become the most popular method used to characterize the elastic properties of organic matter in shale. Table 2 summarizes the reduced modulus obtained in this work alongside those reported in the recent literature. The reported Young's modulus in AFM studies is generally based on a modal value describing a peak in the histogram of measurements, as described previously. This tends to fall in the region of 2-20 GPa, although modal values as high as 30 GPa have been reported (Emmanuel et al., 2016).
The wide range of stiffness values measured on the organic phase and the geomechanical heterogeneity of organic matter stems from chemical compositional differences combined with the physical and chemical effects induced by thermal maturation. Eliyahu et al. (2015) observed a bimodality in the distribution of Young's modulus and linked this to the presence of a soft bitumen and stiffer kerogen phase. However, kerogen also consists of a range of macerals, the molecular structures of which are transformed during thermal maturation. With increasing tempera-CHARLTON ET AL.   ture, kerogen experiences chemical and physical alterations as carbon chains are broken down, releasing petroleum and other high molecular compounds such as asphaltenes. The remaining kerogen changes progressively from an aliphatic-rich mixture to a more structured aromatic complex with increased density (Patience et al., 1992;Okiongbo et al., 2005). The geomechanical implications have recently been investigated using molecular modeling of kerogen at different maturities, which indicates that the increased density, accommodated by greater aromaticity, leads directly to the stiffer mechanical response (Bousige et al., 2016;Kashinath et al., 2020). This reflects the conclusions of combined geochemical and AFM studies where aromaticity has been linked to an increase in the in-situ Young's modulus (Fender et al., 2020;Khatibi et al., 2018). Our findings on Posidonia shale show an increase in stiffness with maturity, from * k E = 5.8 GPa at R o = 0.53% to * k E = 11.3 GPa at R o = 0.89%. Increasing aromaticity has been observed in the Posidonia shale by Bernard et al. (2012), who reported the aliphatic-aromatic ratio of kerogen changing from 0.96 (WIC) to 0.76 (HAR) to 0.52 (HAD). Our results therefore support the link between increased aromaticity and higher stiffness. Table 2 indicates that the increase in the stiffness of organic matter from immature to mature has been broadly observed in several shale formations. It is more challenging to relate Young's modulus to maturity across different formations, that is, to identify a shale-invariant relationship. This may be due to the complex biological origins and amorphous nature of organic matter and the range of physical and chemical processes that can occur with burial, such as the development of secondary organic porosity (related CHARLTON ET AL.   Figure 2a; (b) measured reduced modulus along CS-1 and CS-2 from left to right. The low stiffness values, close to zero, along each cross section show the presence of a cavity and indicate that the histogram peak at 5.8 GPa is not associated with surface topology but represents a material phase (organic matter). to gasification) that could reduce the effective stiffness measured on the organic phase by AFM toward the gas window (Mathia et al., 2016;Suwannasri et al., 2018). Further exploration of AFM-based infrared microscopy on shales (Jakob et al., 2019;Yang et al., 2017), which allows both mechanical and chemical material characterization at a nanoscale, would enable additional insights into the source of mechanical heterogeneities in organic matter and the influence of thermal maturity. Given its variable composition, using a single value to describe the stiffness of organic matter in shale is a simplification but is often necessary for modeling purposes.
CHARLTON ET AL.     (Kashinath et al., 2020), which results in a variation of 9% in E k , as seen in Figure 8. The relationship between Poisson's ratio and maturity is still uncertain; for example, modeling by Bousige et al. (2016) found no clear link whereas Kashinath et al. (2020) observed a higher Poisson's ratio with increasing kerogen density. Both studies showed an exponential increase in Young's modulus as kerogen density increased from 0.9 to 1.5 g/cm 3 .

Nanoindentation
To investigate the mechanical response of the shale composite, nanoindentation testing was conducted. In this test, an indenter with known mechanical properties is pushed into a material with unknown properties. The continuous loading and unloading curves in terms of force, F, and indentation depth, h, are extracted to determine mechanical parameters (Oliver & Pharr, 1992). The indentation is conducted on a grid such that, as for AFM, a statistical characterization of the mechanical response is available (Ulm & Abousleiman, 2006).
Like AFM, the load-displacement curve and surface deformation can be used to obtain a reduced modulus E*, often referred to as the indentation modulus in nanoindentation testing. The reduced modulus is calculated from the initial slope of the unloading curve,   ( / ) h h max S dF dh , which is a purely elastic response: where A c is the contact area, which depends on indentation depth, and β is a correction factor related to the tip shape, equal to 1 for a spherical tip and 1.034 for a Berkovich tip (Fischer-Cripps, 2011). Again, the reduced modulus is related to the sample Young's modulus by Equation 3; a diamond tip is used, meaning the effect of tip deformation can be ignored.
Hardness, which contains both the elastic and plastic response of the sample, can also be obtained and is defined as (Oliver & Pharr, 1992): where F max is the maximum force applied during the test.
For each sample, thin sections parallel (defined as direction-1) and perpendicular (direction-3) to the bedding plane were tested. The indentation was conducted using a Berkovich tip on 7 × 7 grids with a spacing of 100 μm. In this study, a high load of 500 mN (the maximum machine load on a NanoTest Vantage) was applied with the intention of capturing a homogenized response representative of the composite microstructure. Such an approach is sometimes referred to as microindentation to better describe the length scale of the experiment (e.g., Slim et al., 2019). If the goal is to obtain the mechanical properties of the shale constituents, low load nanoindentation with F max ∼ 5 mN is more appropriate.
Each indent affects a characteristic voxel of size 3-5h max (Abedi, Slim, Hofmann, et al., 2016) and given h max did not exceed 10 μm in our experiments, the spacing of the grid was sufficient to avoid interactions CHARLTON ET AL.
10.1029/2020JB021365 10 of 23 It is normal procedure to hold the load at the peak for a few seconds before the unloading stage to allow the system to reach thermal and mechanical equilibrium. Hold periods of 10 s (WIC-7129 and HAR-7046) and 120 s (HAD-7127) were applied. By extending the holding stage, the creep response of the samples may be assessed, as reported in Section 6. ( / ) E E of 1.21. Stiffness anisotropy, with a stiffer response parallel to bedding, is commonly observed in clay-rich shales (Ulm & Abousleiman, 2006) while, in contrast, hardness shows a less obvious relationship to bedding orientation and is often relatively isotropic Bobko & Ulm, 2008). In this sample, H 1 = 0.35 GPa and H 3 = 0.43 GPa (H 1 / H 3 = 0.81).

Reduced Modulus and Hardness
With increasing thermal maturity, the geomechanical response generally becomes stiffer and stronger. The mature HAR-7046 has a similar composition to WIC-7129, although less organic matter is present. The load-displacement plot for HAR-7046, Figure 9d, shows more scatter than recorded on WIC-7129 with indentation depths ranging from 4.5 to 10 μm. The variability is also evident in the coefficient of variation (COV) of E* and H, both increasing significantly when compared to the immature sample. This could be the result of localized diagenetic strengthening and stiffening of the shale matrix, which may also be a cause of the positive skew in the histograms in Figures 9e and 9f, and indicates that the nanoindentation was less successful in capturing a homogenized response. For HAR-7046, * 1 E = 24.5 GPa and * 3 E = 19.8 GPa, both representing an increase of ∼70% from WIC-7129; the anisotropy is only marginally higher ( * * 1 3 / E E = 1.24). In contrast, H 1 increases to 0.54 (a 54% rise) but H 3 drops to 0.38 and H 1 /H 3 is 1.42.
HAD-7127 has significantly lower clay and kerogen content than the samples of lower thermal maturity and is dominated by authigenic carbonate cements. The composition has a clear impact on the nanoindentation results. Figure 9g shows that the maximum indentation depth, between 3 and 6 μm, is substantially lower than in the less mature samples. As reported in Tables 3 and 4, both the reduced moduli ( * 1 E = 51.5 GPa and * 3 E = 39.5 GPa) more than double from HAR-7046 and an even greater increase occurs in hardness. A wide range of values was recorded, as can be seen in the histograms of Figures 9h and 9i, although the COV of both measurements is similar to that of the mature sample. In terms of anisotropy, * * 1 3 / E E shows a small increase to 1.30, which indicates that the stiffness response becomes slightly more anisotropic across the maturity series. In contrast, hardness anisotropy is marginally lower in HAD-7127 (H 1 /H 3 = 1.39) than HAR-7046 and no clear relationship can be discerned with thermal maturity, perhaps due to the variability in the hardness measurements. Figure 10 shows the positive correlation between E* and H, where the relationship is expected to be approximately  * E H (Ulm & Abousleiman, 2006). A specific power law fit to all indents is very similar to that proposed by Mashhadian et al. (2018) on calcite-rich phases of Eagle Ford shale, and the high values of E* and H measured on the carbonate-rich HAD-7127 sample compare particularly well. The soft WIC-7129 sample sits below the power law fit and there is a much clearer difference between indents in parallel and perpendicular directions compared to the more mature shales, being in two distinct groups for the WIC sample.

Homogenization Methods
To obtain the macroscopic stiffness tensor of the anisotropic shale composite, linear micromechanics can be used to homogenize the elastic response through various levels (Zaoui, 2002   where f r is the volume fraction, r  is the stiffness tensor and r  is the strain localization tensor of the r-th phase.
The strain localization tensor is given by: where  is the fourth-order identity tensor and s  is the Hill tensor, related to the well-known Eshelby tensor  : Esh    (Eshelby, 1957), which is a function of particle shape and orientation.
Two commonly used homogenization techniques are the Mori-Tanaka (MT) method (Mori & Tanaka, 1973) and the Self-Consistent Scheme (SCS) (Hill, 1965). For the MT approach, the inclusions are embedded in a matrix phase such that  0 mat   , where mat  is the stiffness tensor of the matrix. For the SCS, the inclusions are embedded into the homogenized medium itself, leading to  0 hom   ; this morphology is described as a granular texture , with no phase dominating the composite, in contrast to the matrix-inclusion morphology of the MT scheme.

Multiscale Model
The formulation of the homogenization scheme follows that of Goodarzi et al. (2017), with the addition of a further level of upscaling to include microcracks. Shale is assumed to consist of silt inclusions randomly embedded into a clay matrix comprising solid clay, kerogen and pores. The multiscale model has four levels: Level 0, the solid unit of clay (at a length scale of ∼1 nm); Level I, the porous clay matrix (100 nm); Level II, silt inclusions embedded in the clay matrix (10-100 μm); Level III, microcracks (1 mm). The model is illustrated schematically in Figure 11.

Level 0 (Solid Clay)
The overall anisotropy of the shale response is included by assuming a transversely isotropic (TI) stiffness tensor for the solid unit of clay. The five elastic components are taken from Ortega et al. (2007): C 11 = 44.9, C 33 = 24.2, C 13 = 18.1, C 44 = 3.7, and C 66 = 11.6 GPa. These values are intended to represent the overall stiffness of the solid clay and are ostensibly insensitive to specific clay mineralogy. This is due to the nanogranular morphology described by Bobko and Ulm (2008), in which the mechanical response of the clay phase is mainly determined by chemo-mechanical interactions between individual clay particles or agglomerates rather than the fundamental behavior of a single particle. The clay fraction of the Posidonia samples is dominated by illite and illite/smectite with secondary kaolinite; this compares well with several of the validation samples used by Ortega et al. (2007). However, there is still uncertainty related to the elastic properties of the solid clay that should be adopted for homogenization schemes. Recently,  and Dubey et al. (2019) have back analyzed the TI elastic constants on various shales and found much higher stiffness values, closer to those obtained for individual clay particles (Wang et al., 2001).

Level I (Porous Clay Matrix)
The porosity of the shale is assumed to exist solely in the clay and kerogen matrix at Level I, meaning ϕ mat = ϕ shale /f mat where f mat is the volume fraction of the matrix components. The elastic properties of the organ-CHARLTON ET AL.

10.1029/2020JB021365
13 of 23  ic matter, assumed isotropic, are taken from Table 2 and converted to Young's modulus using ν = 0.3 (WIC-7129: E k = 5.3 GPa; HAR-7046: E k = 10.3 GPa). Kerogen in the overmature sample HAD-7127 is assumed to have the same properties as measured in the mature HAR-7046 (HAD-7127: E k = 10.3 GPa). This follows the findings of Emmanuel et al. (2016) who reported an increase in Young's modulus from immature shale into the oil window but little change in the stiffness of the organic phase with a further increase in maturity, perhaps due to generation of secondary porosity and hydrocarbons.
Numerical investigation by Goodarzi et al. (2016) showed that SCS performs well at homogenization of the porous clay matrix, and so is adopted here.

Level II (Silt Inclusions)
At the second level, silt particles within the shale composite are included in the homogenization. The particles are assumed to be spherical and isotropic. The insignificant effect of particle shape on the anisotropic macroscopic response is supported by experimental and modeling studies (Ortega et al., 2010). Recent AFM investigations have found the elastic properties of quartz and calcite in shale rocks are generally lower than those normally reported in the literature. Assuming ν = 0.06 for quartz and 0.32 for calcite , mean Young's modulus values of 69 and 55 GPa for quartz and 53 and 40 GPa for calcite were measured by Eliyahu et al. (2015) and Graham et al. (2021) respectively. This is ∼40% less than the 95 GPa (quartz) and 78 GPa (calcite) found in Mavko et al. (2009). The difference may be related to the presence of microporosity within minerals in sedimentary rocks, in contrast to the highly crystalline phases that are typically measured. Given the carbonate diagenesis evident in Posidonia shale, the Young's modulus of dolomite, reported as 117-124 GPa by Mavko et al. (2009), was reduced by the same percentage. The elastic properties are summarized in Table 5.

Level III (Microcracks)
Microcracks, or microfractures, are open fractures that occur in shales because of local failures in the microstructure. Microcracks are characterized by a small aspect ratio, ϵ = c/a, with the crack aperture (half width c) ranging from nanometers to tens of micrometers and length a being one to four orders of magnitude greater, reaching mm-scale. The orientation of microcracks is commonly parallel to the bedding plane, cracks tending to form along weak planes in the rock (e.g., Vernik, 1993), but preferred orientation can also present at various angles, including perpendicular to the bedding plane, due to a range of factors such as stress history, mineralogy and kerogen content (Ougier-Simonin et al., 2016).
The presence of microcracks has been shown to affect the mechanical properties of the rock mass, with a particularly significant effect being the increase of anisotropy (Vernik, 1993(Vernik, , 1994. The effect of cracks can be incorporated into the homogenization schemes described in Section 4.1 through the Hill tensor. As discussed by Qi et al. (2016), a cracked medium can be considered a matrix-inclusion morphology and it is therefore appropriate to adopt the MT scheme. Assuming a single set of open, penny-shaped microcracks, the crack stiffness tensor  c 0  and the homogenized stiffness is: where d is the crack density parameter,  3 d a  and  is the number of cracks per unit volume (Dormieux et al., 2006). Solutions for the Hill or Eshelby tensor have been derived for the case of penny-shaped cracks in a transversely isotropic background medium, the cracks aligned in the direction of anisotropy, by Withers (1989). This solution can be extended to multiple sets of randomly oriented cracks (Giraud et al., 2007;Qi et al., 2016) but attention here is limited to horizontal microcracks, as microCT scans on Posidonia shale indicate this is a realistic idealization (Kanitpanyacharoen et al., 2012).  Table 5 Elastic Properties of Minerals Used in the Multiscale Model (Brown et al., 2016;Mavko et al., 2009;Whitaker et al., 2010) The cracks are assumed to be penny-shaped (ϵ = 0.01) and a crack density of d = 0.07 is used, based on typical values determined by Sarout and Guéguen (2008) during triaxial tests on clay-rich shale. Back analysis from UPV data, including carbonate-rich shales, by Dubey et al. (2019) found a similar crack density (d = 0.043) for the horizontal crack case, indicating the chosen parameters are likely to be representative. The effect of crack fluid is neglected in the multiscale model.

Level II
The high load nanoindentation undertaken in this study is considered to be representative of the homogenized response at Level II. The components of the TI stiffness tensor at this level can be related to * 1 E and The predicted indentation moduli are compared with the experimental values in Figure 12. The plot shows the mean experimental modulus ± one standard deviation. SCS is used at the first level of homogenization (SCS I ), with the results of both SCS and MT shown at the second level (SCS II , MT III ). Figure 12 indicates that the homogenization is relatively successful in predicting the experimental indentation moduli, particularly for WIC-7129. The volume fractions of silt inclusions are 38% (WIC), 56% (HAR), and 79% (HAD); SCS II predicts a higher stiffness than MT II because, as found in numerical simulations by Goodarzi et al. (2016), the predicted stiffness of the two methods diverges as the volume fraction of silt inclusions increases above 40%. The immature WIC sample is a collection of clay, organic matter and carbonate grains and so is well-suited to the formulation of the multiscale model. The effects of carbonate diagenesis in samples of higher maturity pose challenges for upscaling. In HAD-7127, the porous clay matrix has been replaced to a large extent by dolomite and calcite cements. SCS II gives the best estimate for this sample, perhaps because the shale texture can no longer be described as a matrix-inclusion morphology. Although the effects of cementation are not explicitly captured by the homogenization, appropriate choices for the elastic moduli of the various constituents can still enable reasonable estimates of the bulk elastic response to be made. Table 6 shows that the homogenization can predict the anisotropy of the shale response * * 1 3 ( / ) E E quite well for WIC-7129 and HAR-7046, with SCS I -MT II being the more successful. The multiscale model assumes that the anisotropy originates in the solid clay unit. This appears to be true for the clay-rich WIC sample but as the clay fraction reduces from HAR to HAD the estimate of anisotropy gets progressively worse. This suggests that the source of anisotropy must lie elsewhere and may be linked to carbonate diagenesis forming a preferred orientation, invalidating the assumption of isotropic silt grains.   by Rybacki et al. (2015). Young's modulus in directions parallel and perpendicular to the bedding plane can be obtained by Mavko et al. (2009):
As shown in Figure 13, the multiscale model tends to over-predict the elastic modulus at triaxial scale. This may partly be due to the use of a secant modulus, given that in the uniaxial tests the stress-strain response shows an initial stiffening due to closure of existing discontinuities in the sample. In addition, the overestimation may also be attributed to damage inflicted as the sample is loaded. The anisotropy at Level III is shown in Table 7. It is evident that the anisotropy is much higher at triaxial scale than observed during nanoindentation testing, which might be due to the influence of microcracks and fractures. The major effect of including horizontal microcracks in the homogenization scheme is to increase the anisotropy of the response by reducing E 3 . The anisotropy of HAR and HAD is well-captured but a crack density of 0.24 is required to predict the anisotropy of WIC observed in the experiment (E 1 /E 3 = 2.71; SCS I -MT II -MT III ). This agrees with the description of the WIC sample developing abundant microcracks during deformation, whereas HAR and HAD samples failed in a brittle manner along a single localized shear plane (Rybacki et al., 2015).

Nanoindentation Creep
The creep behavior of Posidonia shale was analyzed through nanoindentation testing. The creep tests were conducted as described in Section 4.1 but with a hold period of 120 s at peak load for all samples. A spherical tip was used; the presented data also includes the Berkovich grids on HAD-7127 with a 120 s hold period.
Typical displacement-time responses during the holding stage are shown in Figure 14, from individual indentations. The creep displacement is defined as Δh(t) = h(t) − h 0 , where h 0 is the displacement at the start of the holding stage (t = 0 s). In nanoindentation tests, a logarithmic function is commonly used to fit the contact creep compliance function introduced by Vandamme and Ulm (2013), such that a creep modulus can be obtained. The compliance function describes the creep response due to a unit input of stress and corrects for the geometry of the indenter. The creep compliance function is (Zhang et al., 2014):

Table 7
Experimental (Rybacki et al., 2015) and Predicted Anisotropy of Young's Modulus (Level III) where a u is the radius of the projected contact area between indenter and sample,

and
* 0 E is the indentation modulus at t = 0 s, which could also be measured from the unloading phase as previously. A logarithmic function can be used to fit the creep compliance function as follows: where C is the creep modulus and τ is a characteristic time (Vandamme & Ulm, 2013). The creep modulus describes the rate of creep strain that accumulates, with a higher C indicating a lower creep rate; τ quantifies the time at which the long-term logarithmic creep behavior begins.
The creep compliance function for 7 × 7 indentation grids on each Posidonia shale sample is shown in Figure 15. The variability in the creep response is indicated by a shaded area corresponding to ± one standard deviation of the fitted creep modulus. It is noted that the mean characteristic time is in the range 0.1-0.5 s across all samples, proving the suitability of the creep modulus to describe the longer-term behavior, given that the rate of change of the creep compliance is: The compliance function shows that there is less creep response as thermal maturity increases. The immature sample WIC-7129 has a much greater tendency to creep than the samples of higher maturity. This can be seen in the statistics of C given in Table 8, where a general increase in the creep modulus occurs across the maturity series.
A correlation is evident between the creep modulus and clay/organic content, as shown in Figure 16. Low load nanoindentation creep testing of oil window Vaca Muerta shale (Mighani et al., 2019) has shown that the clay matrix tends to have a low creep modulus, indicating creep deformations are concentrated in this phase. At a larger scale, BSE images on samples of Posidonia (HAR) shale after triaxial creep experiments have shown that (primary) creep is mainly accommodated by pore volume compaction and bending of clay minerals (Herrmann et al., 2020). Sone and Zoback (2013b) note that the correlation is generally evident in samples from a single shale unit but cannot be extrapolated to all shales, as other factors also influence the propensity for a shale to creep. For example, Rassouli and Zoback (2018) suggest that in carbonate-rich samples microporosity within carbonate grains and intergrain microfractures could be sources of creep. In Posidonia shale, carbonate diagenesis may be a key control on creep behavior and the rigid matrix formed by cementation in samples of higher maturity, particularly HAD-7127, may explain the non-linearity of the correlation in Figure 16, and the high creep modulus measured on the HAD sample.
Sample orientation also affects the creep behavior. An anisotropy in creep response of the HAD-7127 sample is evident, with a greater creep compliance (lower creep modulus) when loaded perpendicular to the bedding direction, which again agrees with the conclusions of Sone and Zoback (2013b) in triaxial tests of various North American shales. The  anisotropy ratio C 1 /C 3 is slightly higher than that observed for indentation modulus and hardness (Tables 3 and 4).
It is also clear that the immature WIC-7129 shows a much less variable creep response than the samples of higher maturity, showing a lower COV. This mirrors the results for E* and H and indicates that the creep modulus is representative of a homogenized response for WIC-7129, not greatly influenced by individual constituents or localized diagenetic features such as cementation.
The relationship between creep modulus and hardness is plotted in Figure 17. It is generally observed that C is linearly related to H (e.g., Slim et al., 2019). Despite the high variability in the measured values of H and C on HAD-7127, evident in the large error bars, fitting C = αH through the mean of each indentation grid tested as part of this work gives α = 458, with R 2 = 0.86. This linear scaling can also be seen to provide a good fit with creep data on various US shales of different maturity obtained from nanoindentation tests with applied loads ranging from 4.8 to 50 mN. This confirms the observation of Slim et al. (2019) that α is scale-independent, given the maximum load of 500 mN in this study. The relationship between C and H can be explained by analogy with the logarithmic creep behavior of soils during an oedometer test (Vandamme & Ulm, 2013). The change in void ratio, e, in an oedometer test is linearly related to log(t), with a coefficient of secondary compression defined as C αe = −de/d(log(t)) (Mitchell & Soga, 2005). In one-dimensional compression, the strain rate       (1 )e and, recognizing that the applied vertical stress σ is equivalent to H, the rate of creep compliance is: By comparison with Equation 15: which reveals the expected linear scaling. Given a porosity of ∼0.1, this leads to C αe = 0.0024 which is close to estimates of this coefficient on shale by Mesri and Castro (1987).
The relationship between C and E* has been previously fitted by a power law on US shales, although the correlation is close to linear (Slim, 2016). Figure 18 indicates that there is some difference between Posidonia and US shales, with the Posidonia shale showing a steeper increase in C as the indentation modulus increases beyond 30 GPa. This may be related to the different diagenetic alterations that can occur with increasing thermal maturity, for example the overmature HAD-7127 shows a significant transformation, combined with differences in mineralogy.

Comparison With Triaxial Creep
Long-term triaxial creep tests are normally carried out by holding a constant deviator stress over a time period of up to several days (Herrmann et al., 2020;Rybacki et al., 2017;Sone & Zoback, 2013b. These experiments have shown that shale first exhibits a phase of decreasing creep strain rate (primary creep), with some samples then reaching a constant minimum level (secondary creep) before showing an accelerating strain rate with time (tertiary creep) leading to failure (Rybacki et al., 2017). Primary creep is predominately due to pore space reduction and deformation of weak phases (clay and organic matter) while secondary and CHARLTON ET AL.  tertiary creep is associated with microcracking and ultimately a localized macroscopic failure plane (Herrmann et al., 2020). The occurrence of secondary and tertiary creep depends upon sample composition and test conditions (confining stress, deviator stress and temperature) (Herrmann et al., 2020) but is generally associated with the application of a deviator stress close to the ultimate strength of the shale (Rybacki et al., 2017). Saturation and fluid chemistry also strongly affect the creep response, particularly in clay-rich materials (Cerasi et al., 2017).
To compare the creep response of shale in triaxial and nanoindentation tests, it must be acknowledged that the samples are subjected to different stress states. Zhang et al. (2014) state that in ideal linear viscoelastic materials that creep deviatorically with no volume change, the relationship between the triaxial and indentation creep modulus (triaxial, C txc ; indentation C ind ) is C txc = (1 − 0.5 2 )C ind = 0.75C ind , which results from Equation 3. Shale exhibits a more complex response involving (visco)plastic behavior meaning it is difficult to establish a direct relationship between indentation and triaxial creep. Scale effects will also be important but the creep response during high load nanoindentation, which captures a relatively homogenized response, may be expected to result from similar physical mechanisms as observed during the primary creep phase of triaxial tests, since the applied load is in place for a relatively short period of time and does not cause extensive damage (fracturing) to the sample. As discussed in Section 6.1, the tendency for creep deformations to be associated with the porous clay matrix has been observed in both (low load) nanoindentation (Mighani et al., 2019) and triaxial testing (Herrmann et al., 2020).
The development of creep strain during the primary phase is commonly described by a power law or logarithmic function (Herrmann et al., 2020;Rybacki et al., 2017;Sone & Zoback, 2014), much as for the nanoindentation creep described previously. In a triaxial test, the creep compliance can be written: where ɛ j is a measured strain (axial, radial or volumetric) and σ j is the conjugate loading stress (confining or deviatoric stress). Fitting a logarithmic function to the creep compliance, cf. Equation 14, leads to: where C j and τ j are the creep modulus and characteristic time respectively. The subscript j is dropped in the following, where axial moduli are reported.   Figure 19 shows that the creep response during high load nanoindentation gives similar values of the creep modulus as recorded at triaxial scale. The general trend suggests a positive correlation between the creep and elastic modulus indicating that stiffer shales tend to creep less, a finding also observed by Sone and Zoback (2013b) and attributed in part to stress partitioning effects between soft and stiff material phases.
For the HAR shale, loaded perpendicular to the bedding, the mean triaxial creep modulus is 92 GPa. This is close to the value obtained from nanoindentation testing on the same shale (C = 115 GPa, parallel to bedding), suggesting that the high load nanoindentation gives a similar creep response to the triaxial test. The anisotropy C 1 /C 3 would be 1.25, which is comparable to the anisotropy observed from nanoindentation testing on HAD-7127.

Conclusions
This study reports a multiscale investigation of carbonate-and organic-rich Posidonia shales with thermal maturities ranging from immature to overmature. Nano and microscale mechanical tests show that the macroscale geomechanical properties of shale are the result of the interplay of mineralogy, microstructure and diagenesis. The dissolution of biogenic carbonate and its reprecipitation as cement at higher maturities is a particularly important control on the diagenetic evolution of mechanical properties in these samples.
Peakforce QNM successfully determined the elastic modulus of organic matter at a resolution of a few nanometers. The reduced modulus of immature kerogen was * k E = 5.8 GPa, increasing to 11.3 GPa in the oil window.
High load nanoindentation testing of the shale generated the indentation modulus and hardness of the shale composite at a microstructural level. The immature WIC sample gave the softest response while diagenetic cementation with increasing maturity had the effect of increasing the stiffness and hardness of HAR and HAD samples. A bulk response could be obtained on the softer, immature shale with the maximum applied load of 500 mN. However, the indentation results showed greater variability for the more thermally mature samples as the diagenetically-strengthened microstructures prevented a truly homogenized response.
Upscaling of the elastic properties through a multiscale model showed a good prediction of the indentation modulus. A particularly successful estimate was achieved for the immature shale, for which the bulk mechanical properties can be related straightforwardly to the individual constituents. Diagenetically modified and cemented mature samples are less simple to describe but good estimates of the bulk elastic response can still be made with appropriate choices of the elastic properties of the constituents, incorporating to some extent the diagenetic effects. To improve confidence in the mechanical properties, a more complex (but less general) formulation would be needed, specifically accounting for diagenetic features such as carbonate cementation. Horizontal microcracks were included in a third level of homogenization and improved the prediction of Young's modulus in uniaxial tests, especially the anisotropy of the response.
Creep behavior was investigated by holding a constant load for a specific time period in nanoindentation tests. The concept of the creep modulus was used to interpret the time-dependent response. Less creep was observed with increasing thermal maturity and there was a negative correlation between clay/organic content and the creep modulus, indicating that the deformation of these soft phases has a strong influence on the creep response. This agrees with observations made from microstructural imaging of triaxial creep tests on Posidonia shale. The creep response during high load nanoindentation and the primary creep phase of triaxial tests generated similar values for the creep modulus, suggesting the same deformation mechanisms are active. Further investigation of the effect of in-situ conditions (pressure, saturation, temperature etc.) at different scales may help to better constrain the creep response. Figure 19. Creep modulus versus Young's modulus. Solid markers: perpendicular to bedding. Empty: parallel to bedding. Data from Mighani et al. (2019) and references therein, except Bowland (Herrmann et al., 2020) and WIC/HAR/HAD (this study; Herrmann et al., 2020;Rybacki et al., 2015). Black markers indicate C from nanoindentation; gray markers indicate C from triaxial test.