The Assessment of Organic Matter Young's Modulus Distribution With Depositional Environment and Maturity

Quantification of risk to seal integrity in CCS, or gas extraction from hydraulic fracturing, is directly affected by the accessibility of organic pores within organic rich mudrocks. Knowledge of the host organic matter's mechanical properties, which are influenced by depositional environment and thermal maturity, are required to reduce operational risk. In this study we address the effect of both depositional environment and maturity on organic matter Young's modulus by means of Atomic Force Microscopy Quantitative ImagingTM, which is a nondestructive technique capable of nanomechanical measurements. Shales from varying marine depositional environments covering kerogen Types II (Barnett), IIS (Tarfaya), and II/III (Eagle Ford/ Bowland) are analyzed to capture variance in organic matter. The findings show organic matter has a Young's modulus ranging between 0.1 and 24 GPa. These marine shales have a bimodal distribution of Young's modulus to some degree, with peaks at between 3–10 and 19–24 GPa. These shales exhibit a trend with maturity, whereby Young's modulus values of <10 GPa are dominant in immature Tarfaya shale, becoming similar to the proportion of values above 15 GPa within the oil window, before the stiffer values dominate into the gas window. These peaks most likely represent soft heterogeneous aliphatic rich kerogen and stiff ordered aromatic rich kerogen, evidenced by the increase in the stiffer component with maturity and correlated with 13C NMR spectrocopy. These findings enable increased realism in microscale geomechanical fracture tip propagation models and may allow direct comparison between Young's modulus and Rock‐Eval parameters.


Introduction
Quantification of mudrock geomechanics has become increasingly important in order to inform the design of hydraulic fracturing treatments. This is driven by the need to reduce economic risk and balance environmental impact. Additionally, in order to assess the suitability of conventional oil and gas reservoirs along with saline aquifers for storage of carbon dioxide, geomechanical models are required to constrain the conditions under which the top seal may begin to fail through fracturing (Ingram & Urai, 1999). Laboratory measurements and geomechanical modeling become paramount, as production from reservoirs causes changes in the effective stress regime, triggering seismicity and inducing new fractures and fault reactivation, potentially compromising seal integrity. In an effort to minimize risk and improve understanding, A large proportion of shale geomechanical modeling has been dedicated to predicting bulk rock properties. The majority of hydrocarbons are however stored in interparticle pores within the matrix, or within intraparticle pores associated with the organic matter (Wang & Reed, 2009). The most likely migratory route is through these pores for nonfracture-related seal release. The degree of anisotropy within whole rock shale models is in a direct correlation with the organic matter content, which is expected to have a key control on shale geomechanical properties (Sone & Zoback, 2013). The bulk matrix properties of a shale can be determined from triaxial testing or high load nanoindentation tests. Due to the limited resolution of these experiments (greater than 5 cm 3 and 5 μm 2 , respectively), and the size of interspersed organic matter (which can be less than 1 μm in diameter), obtaining representative measurements of organic matter Young's modulus is challenging and relies on measurements of larger organic particles. Therefore, novel techniques must be applied to determine the mechanical properties of these individual components.
The Young's modulus of organic matter has previously been measured by means of nanoindentation and Atomic Force Microscopy (AFM), with results indicating a range of 1-32.2 GPa (Table 3). Nanoindentation studies have found a range of values of approximately 1-20 GPa for kerogen Young's modulus (Table 3 and references therein). The aforementioned limited resolution of nanoindentation precludes it from being used to resolve the mechanical properties of most dispersed bitumen along with the majority of internal heterogeneity in kerogen. To further complicate the interpretation of nanoindentation measurements, organic matter porosity has been inferred to cause variability in the Young's modulus (Kumar et al., 2012), since porosity will cause a reduction in Young's modulus relative to a solid unit (Bobko & Ulm, 2008). Submicron scale organic porosity is abundant in organic matter across the maturity range (Curtis et al., 2012;Loucks et al., 2009) but generally increases with maturity.
Kerogen undergoes both chemical and physical changes with maturity, becoming increasingly aromatic, losing a large proportion of its thiophenes, amines, and other nitrogen, sulphur, oxygen (NSO) compounds. This is coupled with kerogen becoming increasingly porous while also reducing in overall volume (Curtis et al., 2012;Okiongbo et al., 2005;Suwannasri et al., 2018). Pores within the kerogen are primarily nanometer intraparticle pores, with secondary porosity created where the kerogen volume has decreased (Bernard et al., 2012). Several studies have identified that organic matter develops secondary porosity with maturity, which can range to nanoscale intraparticle circular pores (Bernard et al., 2012) or interparticle elongate fissures (Allan et al., 2014(Allan et al., , 2016(Allan et al., , 2018. This secondary porosity accounts for the majority of porosity increase within shale samples with maturity and can comprise up to 6-7% volume (Suwannasri et al., 2018). The onset of this porosity increase occurs during the oil window, with a large proportion already generated by the wet gas window. Modeling work by Goodarzi et al. (2017) shows that porosities up to 20% are possible at the scale of a nanoindent and could reduce the measured Young's modulus by up to 60% of its solid value, which was supported by nanoindentation measurements across a series of maturities by Zargari et al. (2016).
AFM studies mainly using Bruker's AFM PeakForce™QNM™ have measured a maximum organic matter Young's modulus of 32.2 GPa (Table 3 and references therein). Organic matter has been shown on occasion to comprise a stiff outer coating with a softer porous center (Eliyahu et al., 2015;Emmanuel et al., 2016a).
Effective elasticity is affected directly by porosity in studies using both acoustic wave velocity (Allan et al., 2014(Allan et al., , 2016Suwannasri et al., 2018) and nanoindentation, whereby derived organic matter Young's modulus decreases from 20 to 7-12 GPa (Zargari et al., 2016). However, in higher-resolution studies using AFM, an increase in organic matter Young's modulus is observed with maturity (Emmanuel et al., 2016a;Goodarzi et al., 2017;. Further variability in organic matter Young's modulus has been identified using AFM, with a bimodal distribution forming at higher maturities, believed to reflect the difference between the stiff kerogen component and softer bituminous component (Emmanuel et al., 2016a). This phenomenon was also observed using nanoindentation of artificially matured Bakken shale by Zargari et al. (2011), who attributed a significant amount of Young's moduli values below 10 GPa to bitumen formed during hydrous pyrolysis.
In this study we use Atomic Force Microscopy Quantitative Imaging (AFM QI™) for nanomechanical measurements. The benefit of this technique is both the increased resolution (10-50 nm vs. 1 μm) compared with nanoindentation measurements and the nondestructive nature of AFM QI™. This allows imaging with scanning electron microscopy after measurements have been acquired. The range of Young's modulus values in AFM nanomechanical studies of organic matter shows greater heterogeneity, both between organic particles and also within an individual organic particle, resulting in a greater variance than that obtained through nanoindentation (1-32.2 GPa vs. 1-11.1 GPa; see Table 3), and thus more likely represents the true variability of organic matter Young's modulus.
As mentioned previously, the degradation of organic matter with progressive thermal maturation (burial) is well understood and documented; it progresses from an aliphatic rich organic mélange into a better structured aromatic complex. The rate of this degradation is dependent on the primary composition of the organic matter, which in turn is a function of depositional environment and preservation efficiency (Pepper & Corvi, 1995a, 1995bPepper & Dodd, 1995). This raises the question of whether the rate and extent of thermal modification during maturation are determined by organic precursor input and if this could affect the mechanical properties of organic matter. Although it is noted that geochemical maturity changes also impact Young's modulus, little research has been focused on the impact of different geochemical input characteristics on Young's modulus.
We hypothesize that geochemical variability from subtle differences in depositional environment and thermal maturity cause intrinsic patterns within the Young's modulus distribution of organic matter. Using this technique, we are able to investigate trends in Young's modulus that represent the transition of organic matter through different thermal maturities. These distributions may allow for increased realism in microscale geomechanical models, which will assist in the interpretation of fracture tip propagation in the subsurface. This information is a first-order control on the fracture mechanics of organic pores, which are known to house the hydrocarbons or CO 2 enabling increasingly cost effective hydraulic fracturing and CCS seal integrity studies.

Depositional Environment and Geological Setting of Studied Samples
In order to conduct this study, shale samples were chosen from a range of marine depositional environments, which include some of the most prolific black shale depo-types known in the petroleum industry.

Barnett Shale
The Mississippian Barnett Shale Formation was deposited by hemipelagic processes in a shallow epeiric sea  covering the majority of north and central Texas. The Barnett Shale is, at its maximum ∼300 m of relatively homogenous black shale, except for a debris flow deposit in the center (Ellenburger), itself with a maximum thickness of 400 ft in the northeast of the basin. The Barnett Shale was deposited at approximately 14 μm/year (Loucks & Ruppel, 2007) and has average TOC content of 2.5-3.5% (Jarvie et al., 2005). Studies into the sedimentary structures in the Barnett Shale suggest a shallow paleo-water depth, below 120 m (Loucks & Ruppel, 2007). The Barnett shale was suited for this study as a classical shale composed almost entirely of hemipelagic deposits giving rise to a characteristic Type II kerogen . The Barnett shale also has a maturity signature not completely dominated by present day overburden. The Barnett is anomalously mature in the northeast and southeast of the Fort Worth Basin, due to the impact of hydrothermal fluids during the time of the Ouachita Thrust in the Late Carboniferous-Early Permian (Bowker, 2007). This allows an insight into the shale at similar depositional ages, from an almost identical environments but at different levels of thermal maturity. Samples in this study were selected accordingly (  .

Bowland Shale
The Bowland Shale in England is also of Mississippian age and deposited in the Bowland Basin (covering large parts of northwest Lancashire), which is one of a series of small basins occupying tilted fault blocks formed during rifting in the late Devonian/mid Chesterian. These fault blocks were tectonically inactive during the deposition of the Bowland Shale in the middle Mississippian (Fraser et al., 1990). The Bowland Shale is mainly composed of hemipelagic deposits, although a southward prograding prodelta, responsible for deposition of the Millstone Grit (Pennsylvanian), meant an increase in terrestrial organic matter present. Facies changes within the Bowland are related to basin topography during the Late Mississippian/Early Pennyslvanian, which is observed in either orbitally driven gradual changes or abrupt gravity slump deposits (Gawthorpe, 1987). The Bowland Shale interval studied here is from the Lower Bowland Shale (mid-Mississippian) and represents a maximum flooding surface, where anoxia is believed to be widespread within the basin (Gawthorpe, 1987). This choice of interval maximizes the potential of measuring organic matter during nanomechanical testing, with the range of TOC of the Bowland in the Preese Hall-1 well sampled being from 0.5-6.1wt% (Fauchille et al., 2017), compared to 6.47wt% in the sample tested (Table 1). Records suggest that the Lower Bowland shale within the Bowland Basin is oil to gas window mature (Andrews, 2013).

Eagle Ford Shale
The Eagle Ford Shale was deposited during the Cenomanian-Turonian Ocean Anoxic Event (OAE II), in the Western Interior Seaway ranging from Canada across the majority of the western USA to Mexico. The Eagle Ford was deposited in three major depocenters: one to the east and one to the west of the San Marcos Arc paleohigh, along with a further depocenter in the rapidly subsiding Maverick Basin (Condon & Dyman, 2006). This study focuses on samples of marginally mature Eagle Ford from the Cinco Saus well in the north of the Maverick Basin and mature Eagle Ford shale from the Carol Well in Lavaca County from high on the east limb of the San Marcos Arch, south of the Karnes trough (Condon & Dyman, 2006;Fan et al., 2011). These samples represent different thermal maturities but also subtly different deposition environments, with the Eagle Ford deposited closer to the shore in the Carol well, which is consistent with the increased Type III kerogen content reported in the Eagle Ford around the San Marcos Arch (Condon & Dyman, 2006), which is a similar depositional environment to the Bowland Shale.

Tarfaya Shale
Like the Eagle Ford, the Tarfaya shale was also deposited during the Cenomanian-Turonian OAE II and is composed primarily of fine-grained calcareous material. The Tarfaya shales were, however, deposited in warm upwelling nutrient rich water off northwest Africa, enhancing primary productivity, and are generally richer in organic matter (Sachse et al., 2011), but this is highly variable at orbital time scales (Kolonic et al., 2002;Wagner et al., 2013). The studied interval shows strong biomarker and trace metal evidence support of temporal fluctuations in bottom water and photic zone euxinia (Kolonic et al., 2002). The Tarfaya shale is immature in the samples studied but makes an excellent potential source rock further offshore where it is mature, with organic matter concentrations of up to 16% (Sachse et al., 2011), and with the majority of TOC contents falling between 3% and 9% for the section sampled from Kolonic et al. (2005). The organic matter within the Tarfaya is highly labile and known to be sulfurized, fitting into the characteristic Type IIS region of a pseudo-Van Krevelen diagram with HI values of between 400 and 800 mgHC/gTOC reported (Kolonic et al., 2005;Sachse et al., 2011). Organic matter is reported to be predominantly amorphous, with minor concentrations of alginite in the form of lamalginite (Sachse et al., 2011).

Preparation and Analysis
Shale samples, prepared as thin sections perpendicular to bedding, were analyzed using scanning electron microscopy (SEM), with a FEI Quanta FEG 650 SEM. The microscope was operated under low vacuum conditions (0.82 Torr), at between 15 and 20 kV, allowing a resolution of <100 nm on the uncoated samples.

Organic Matter Chemical Characterization
Total organic carbon (TOC) and Rock Eval Pyrolysis measurements were conducted to assist in maturity analysis and organic matter typing. Rock Eval analysis was undertaken using a Rock-Eval 6 configured in standard mode, with a heating rate of 25 • C/min between 300 • C and 600 • C in an inert N 2 atmosphere, and the residual organic matter was then oxidized from 300 • C to 850 • C at 20 • C/min. The hydrocarbons released during pyrolysis were measured using a flame ionization detector and CO and CO 2 measured using an infrared dector. Oil saturation index (OSI), Hydrogen Index (HI), and Oxygen Index (OI) were calculated using the following equations (1): The values of S1, S2, and S3 were measured by the integration of the respective peaks, with S1 representing in situ hydrocarbons, S2 representing the cracking of organic matter to hydrocarbons during heating, and S3 representing the organic matter analyzed as CO and CO 2 . Further details of this method are found in Słowakiewicz et al. (2014).
Solid-state 13 C NMR spectroscopy was undertaken on powdered samples to help characterize the organic composition of each shale. The spectra were recorded at a frequency of 100.62 MHz using a Bruker Advance III HD spectrometer and 4 mm magic angle spinning probe. The spectra were recorded at ambient temperature, using cross-polarization and TOSS spinning sideband suppression with a 1 ms contact time, 2 s recycle delay, and at a sample spinning rate of 10 kHz. For each NMR analysis, 100-165 mg of sample was used. The spectral referencing was obtained by using neat tetramethylsilane carried out by setting the high-frequency line from adamantane to 38.5 ppm. The absolute value of signal varied between each sample, which was a reflection of the variable organic carbon content across the range of shales sampled.
NMR spectra were analyzed using chemical ranges from Chadwick et al. (2004), whereby aliphatic hydrocarbons exist primarily of the spectral range of 0-55 ppm and aromatic hydrocarbons over 95-170 ppm.
Although solid-state NMR is not absolutely quantitative, it can be considered a semiquantitative method of characterizing organic composition.

Palynofacies
Palynofacies preparation was undertaken on approximately 5 g of sample from each shale using hydrochloric acid and hydrofluoric acid to remove carbonate and silicates. The kerogen fraction was then filtered FENDER ET AL. 5 of 24 . After this, a set amount of force is applied (c). It is on this gradient that the Young's modulus is measured (red line). The amount of displacement (d − d 0 ) is then measured at the largest force (F). After which the tip is removed from the surface by the oscillation of the cantilever (d).
through a 10 μm sieve and strewn-mounted on microslides. Palynofacies analysis was undertaken on a Leica DM2700 P polarizing microscope at 400X magnification, equipped with a PetrogLite motorized stage. Three-hundred randomly selected particles were analyzed following the palynofacies classification of Tyson (2012) and McArthur et al. (2017). Slides were also studied in blue light excitation following the standard recommendations for epifluorescence on palynological slides (Tyson, 2006(Tyson, , 2012.

Nanomechanical Measurements
The bedding-normal shale samples were prepared for nanomechanical measurements using a GATAN Precision Etching and Coating System (PECS™) II-Argon Broad Ion Beam Mill on to reduce shale topography. The central 28 mm 2 of the sample sections was polished for 90 min using a 4 kv ion beam to achieve maximum results without excess spluttering.
A JPK NanoWizard 3 (now Bruker JPK BioAFM™) was used to generate nanomechanical measurements using QI™. This technique along with its sister technique PeakForce TM QNM TM has successfully obtained mechanical measurements of biological materials (Chopinet et al., 2013), cements (Trtik et al., 2012), and shales (Goodarzi et al., 2017). Recent research suggests that these two techniques are similar when measuring moduli in the range expected from organic matter (Graham et al., 2020). The AFM QI TM was used with the following settings: Force applied: 500 nN; Z length: 750 nm; and speed: 253 μm/s. Silicon Nitride RTESPA-525 tips were selected, which are the latest variant of the TAP-525 used by Emmanuel et al. (2016a). These tips are indicated to have a nominal working range of between 0.1 and ∼20 GPa (Pittenger et al., 2013). However, Emmanuel et al. (2016a) measured Young's modulus values of up to 32 GPa using similar TAP-525 tips, suggesting that there is flexibility to measure Young's moduli in this region with the RTESPA-525 tips.
The general method for acquiring AFM mechanical measurements involves oscillating a tip/cantilever assemblage onto a less stiff sample surface. The oscillation occurs at a frequency below the resonance frequency with a set force. Movement of the tip-surface interaction is measured by deflection of a lazer on a photo-diode using beam theory. More information on the QI™ is described in QI™ mode-quantitative imaging with the NanoWizard 3 AFM (JPK-Instruments, 2011) and in Figure 1.
New tips were used per four to eight scans due to tip wear, or sooner if the tip was believed to be damaged during a scan. Spring constants and sensitivities were calculated using contact-free calibration, whereby the tip is oscillated away from the sample and the flextural movement of the cantilever is recorded before use. This is the method suggested by Bruker and in a comparison with contact-mode calibration with a sapphire standard produced more realistic spring constants and sensitivities. Young's modulus calibration was obtained before scanning with a new tip by scanning a sample of Highly Ordered Pyrolitic Graphite (HPOG-12M), adjusting sensitivities and the spring constants to obtain a value of 18 GPa ± 0.5 GPa.
Random areas of each shale were scanned within the 28 mm 2 ion beam milled area to ascertain the Young's modulus of organic matter. Areas scanned were generally 250 μm 2 at a resolution of 512 × 512 pixels. Each sample was successfully scanned a minimum of 10 times. This approach has been chosen as a means of upscaling measurements of individual organic matter particles from literature (Emmanuel et al., 2016a;Goodarzi et al., 2017;Khatibi et al., 2018) and attempt to align the-high resolution data from the AFM with other bulk sample measurements (Rock Eval). To locally upscale, some AFM scans were meshed together to form a larger 1,000 μm 2 area. While this only reflects a small proportion of the 28 mm 2 milled area, the total area scanned of approximately 2,500 μm 2 may offer the first insights into upscaling AFM results and measuring organic matter as a matrix, rather than a grain.
An advantage of the AFM QI™ is the user has a constant feedback of force curves, which allows understanding of tip functionality. If during scanning a tip was thought to be damaged, the tip was replaced immediately and the area rescanned if possible.
Once the force curve has been produced, the reduced Young's modulus is derived from the fitting of the extend curve, in accordance with the operation manual for the Bruker BioFilm AFM QI™. This discounts the minor amounts of adhesion on the samples (adhesion is corrected for in processing of results). These samples were then analyzed using the Hertz-Sneddon model (Sneddon, 1965): where F is the peak force, E ′ is the reduced Young's modulus, R is the tip radius (approximately 10 nm for RTESPA-525 tips), and (d − d 0 ) is the sample displacement ( Figure 1). Using Equation 2, Young's modulus can be calculated from the reduced Young's modulus, using the Poisson's ratio of the sample ( s ) and tip ( tip ) along with the Young's modulus of the tip-cantilever array (E tip ).

Modulus Data Analysis
A four-step approach was undertaken in the processing of the raw AFM QI™ data (Figure 2. • The initial reduced Young's modulus data are processed through a 30 GPa cutoff, which represents the reported absolute limit of the tip from Emmanuel et al. (2016a) (Stage 1 in Figure 2). At this stage any data below reported lower limit of the tip are removed (0.1 GPa; Pittenger et al., 2013). • The residual data are then checked for a maximum value of organic matter Young's modulus, achieved by comparing the AFM QI™ scans with their equivalent SEM scan and EDX data obtained. This was undertaken first by enhancing the SEM through color thresholding using the method proposed by Camp (2013). The pixel values of the enhanced SEM are then manually compared with visible values for organic matter and porosity. The mask created is combined with a mask of the colors associated with the peak carbon values from EDX (Stage 2 Figure 2). This allows removal of dark colors in the SEM, which are associated porosity, as these values will show no peak in carbon on the EDX. • The data are then put through an adhesion filter, which removes values where adhesion is greater than 5% of the total force exerted. This level of adhesion makes the value suitable for the Derjaguin-Muller-Toporov model (Derjaguin et al., 1975). • Lastly, any areas of topographic variability above 5 • are filtered out in all shales apart from the Tarfarya, in which the requisite polish from ion beam mill could not be achieved due to the high proportion of organic matter and contrast with calcite. Further, ion beam milling may have caused artificial maturity effects of the organic matter.  Upon filtering the Young's modulus data of each scan and each shale, the results were visualized as a Kernel Density Estimate (KDE) plot in order to identify trends using the GetDist package for Python (Lewis, 2019), which contains inbuilt functionality to undertake boundary correction on Kernel Density Estimates.

Deconvolution
The deconvolution technique first described by Bobko and Ulm (2008) has been used frequently to understand the complex distributions in nanoindentation data. This technique uses the central limit theorem, whereby the distributions of a series of means will always be normally distributed (Bauer, 2011). In the case of nanoindentation, each indent measures a mean Young's modulus across the indented area (∼1-5 μm 2 ). Therefore, once a series of nanoindentation results are plotted, the resulting probability density function represents a composite of the normal distributions for each mechanically active phase.
These phases can be described by their proportion on the surface, along with their mean and standard deviation. This technique is useful when dealing with effective moduli, such as at the scale of a single nanoindent. However, it requires some adjusting when downscaling the measurements to local moduli values, which are not necessarily normally distributed.
Here we use the deconvolution technique to understand the potential different distributions within the organic matter Young's modulus. Initially, the experimental cumulative distribution function (ECDF) of the organic matter Young's modulus data is calculated using Equation 3.
where N is the number of data points in the modulus data set, sorted into ascending order, and i is an incremental point on the data set. First, the ECDF is compared to a two-phase normally distributed system, the CDF of which is generated using Equation 4.
where is the mean and is the standard deviation. The superscript j represents the jth phase in the system.
The CDFs generated of the two phases are minimized for the smallest root-mean-squared error (RMSE) against the ECDF using Equation 5. 10.1029/2020JB020435 In this each phase (F) occupies a fraction( f ) of the sample surface. This minimization is constrained by the sum of all phases equating to one (Constraint A), which represents a physical constraint on the model. A second analytical constraint (Constraint B) is selected to avoid overlapping phases. This is important in systems whereby the mechanical phases are unrelated. However, in this deconvolution we present results with both this constraint in place but also relaxed, as the phases of organic matter may be fundamentally related (kerogen hosting bitumen or aliphatic and aromatic end members of carbon).
Previous research suggests ) that the normal distribution may not be the most suitable fit for coal macerals, which are similar in structure to the organic matter in shales. As such, in accordance with those results, minimization has been undertaken using the CDF for the gamma distribution (Equation 7) to indicate if this prevails within organic matter within shales.
4. Results The TOC of the Tarfaya shale sample is the highest observed at 17.5%, which is indicative of a TOC rich layer in this section (Poulton et al., 2015). The Bowland shale TOC value, which at 6.5%, is also among the richest intervals in the entire section logged (0.51-6.1%wt reported in Fauchille et al., 2017).

Organic Content Typing and Analysis
The inherent variability between shales is demonstrated in mineralogy, with the Tarfaya and Eagle Ford displaying the characteristics of a marl rather than the typical black shale. Tmax values indicate that the Tarfaya sample is thermally immature (<435 • C), while the Bowland has a Tmax greater than 460 • C indicative of the gas window. The Barnett TP Sims has been reported to be in the gas window (Zhang et al., 2014), and although the Tmax is below 460 • C, the low S2 and OSI values suggest that this sample is also within the gas window. The Eagle Ford Carol shale is at approximately late oil maturity, represented by high OSI, due to little expulsion. The Barnett and Cinco Saus Eagle Ford shales have reached marginal oil window maturity (Tmax 445 • C and 443.5 • C, respectively).
Palynofacies characterization (Table 2) supports previous depositional environment assignments. The Barnett and Tarfaya samples composed solely of amorphous organic matter, typical of algal-dominated hemipelagic deposition (Kolonic et al., 2002;Loucks & Ruppel, 2007). The Bowland sample contains the most terrestrially derived organic matter with 2% vitrinite and 3.3% inertinite. Both the Eagle Ford samples contained similar maceral compositions with the Carol slightly enriched in vitrinite (three counts per section compared to two). The proportion of terrestrial organic matter in the Carol and Bowland samples may be lower than anticipated due in part to sampling from organic rich intervals of the core, representing periods of less continental runoff and increased preservation of marine organic matter.
SEM imaging perpendicular to bedding indicates a subtle difference in the appearance of the organic matter between the shales studied ( Figure 3). The grain fabric in the Barnett Lake Davis sample (Figure 3a) appears to be random, and the organic matter distribution is almost entirely bimodal and occurs either among clay-rich interlayers or as scarce single agglomerations within the well-cemented quartz rich layers. The organic matter within the TP Sims (Figure 3b) is subtly different to that of the Lake Davis sample and occurs as either larger (∼50μm) or smaller (∼10μm) irregularly shaped lenses. The fabric of TP Sims sample shows little lamination within the organic matter or other minerals, similar to the Lake Davis sample, indicative of possible early cementation. The organic matter in the Bowland sample ( Figure 3c) appears to be either in the form of large (>10μm) spores or smaller lenses of AOM, distributed parallel to bedding. Across the Bowland, the organic matter appears to be strongly laminated with the fabric.
The habit of organic matter within the Eagle Ford is similar in both the Carol sample ( Figure 3d) and the Cinco Saus sample (Figure 3e). Organics generally appear as either as 10-40 μm clumps within a fecal pellet or as thinner wispy lamellae within the matrix of the micrite interlayers.
The organic matter appearance in the Tarfaya is similar to the Eagle Ford; both are primarily formed from fecal remains (Figure 3f). While similar in morphology of organic matter, the Tarfaya has a significant proportion existing as ∼10μm lamellae. The Tarfaya is also finer grained and has an organic-dominated matrix with interspersed micrite.
Results from Solid State NMR spectroscopy, displayed in Figure 4, show the proportion of aliphatic to aromatic hydrocarbons within the kerogen of each shale. The proportion of aliphatic to aromatic kerogen within each shale has been determined by integrating the values at 0-55 and 90-145 ppm, respectively, from the   NMR results. The thermally immature Tarfaya shale is significantly enriched in aliphatic carbon (70%) compared to aromatic carbon within the kerogen results shown in Figure 4. The Cinco Saus, Carol, and Barnett Lake Davis all comprise approximately 68-69% aromatic kerogen, whereas the Bowland and Barnett TP Sims contain >80% aromatic kerogen.

Organic Matter Young's Modulus
The Young's modulus values obtained within this study along with other studies are included Table 3, which shows that the range of values found within this study are similar to other AFM studies. Random areas of each shale sample were scanned using the AFM QI™ and the overall reduced Young's modulus results are displayed in Figure 5. The conversion of E ′ to E from Equation 2has been used on a scan from the Barnett Shale Lake Davis and Eagle Ford Carol samples (Figure 6), with a change in resulting in a shading represents the 95% confidence interval for each shale, and the number in blue or black refers to the amount of AFM scans used to create the distribution once filtered. There appears to be two distributions for the Barnett Lake Davis and Eagle Ford Cinco Saus shales, which are at peak and early oil window maturities, respectively.  (Eliyahu et al., 2015;, but it is likely that Poisson's ratio will vary in a similar manner to Young's modulus between shales, between organic particles, and within organic particles. Therefore, further work must be undertaken to better constrain the Poisson's ratio of organic matter, and because of this, we have chosen to report our distributions in terms of reduced Young's modulus (E ′ ).
The modulus distributions in Figure 5 show the average across each shale sample and 95% confidence interval. Even though a minimum of 10 AFM scans have been undertaken on each shale, not all scans were suitable for use after filtering, owing to topographic issues or large amounts of adhesion. Thus, the number of scans is included in Figure 5 at the top of each distribution. As can be seen, the Bowland and Tarfaya were most affected by adhesion and topography, with the Bowland showing significant adhesion believed to be from clay mineral interactions. The Tarfaya samples also showed moderate adhesion within areas dominated by organic matter, along with a polishing issue mentioned previously.
The Young's modulus range of the Barnett shale was between 0.1 and 23 GPa. The Lake Davis scans ( Figure 5a) exhibit two different probability distributions: Distribution 1, bimodal with peaks at 6.0 and 20.0 GPa, but with the stiffer distribution dominant, and Distribution 2, bimodal with the initial distribution showing a subpeak at 6.6 GPa but a major peak at 9.8 GPa and a second distribution with a modal value of 20.0 GPa. This distribution has a high proportion of values within the less stiff (9.8 GPa mode) peak. The Barnett TP Sims sample also exhibits a bimodal distribution, but with a significantly lower proportion of values less than 10 GPa, centered on a modal value of 8.9 GPa. The modulus peak in the stiffest phase is at 22.8 GPa. The TP Sims organic matter also has a greater maximum stiffness and appears to be more uniformly distributed than the organics in the Lake Davis, with a small 95% confidence interval (Figure 5b).
The Bowland shale Young's modulus distribution ranges from the cutoff minimum (0.1 GPa) to 24 GPa and comprises two main peaks, with a small proportion of organic matter data situated around a modal value of 8.8 GPa, and another more significant peak at 24.0 GPa (Figure 5c). The Bowland appears less bimodal than the Lake Davis Barnett, but with a better defined initial peak than exhibited in the TP Sims distribution.
Scans of the Carol well of the Eagle Ford shale indicate organic matter Young's modulus values from 0.1 to 20 GPa (Figure 5d). The Carol, similar to the Bowland and Barnett TP Sims, displays a bimodal distribution in the main dominated by a stiff peak with a modal value of 20.0 GPa and a lesser peak at 3.3 GPa. This value is below the lower modal value in the Bowland and TP Sims. The most significant peak in the Carol well organic matter is centered on a modal value of 20.0 GPa contributing approximately 42% of the overall data. The Cinco Saus well of the Eagle Ford has a different distribution to the Carol and the other aforementioned distributions with an equal bimodal distribution (Figure 5e), with a similar number of values within soft and stiff subdistributions: 43.7% of data >10 GPa and 56.3% of data <10 GPa. This is represented in Figure  5e as the two distributions and accompanying 95% confidence intervals. The modal values for the softer Figure 7. The resulting phase distributions of a two phase deconvolution on each shale distribution both with Constraint B in place (constrained) and Constraint B relaxed (unconstrained). The distribution of the phases in constrained and unconstrained deconvolution is identical using normal distributions. There is a large disparity of fit between the constrained and unconstrained gamma distributions. The gamma distribution appears to poorly fit the Tarfaya shale in all cases, while the normal distribution only fits poorly for the stiffer phase indicating that this may not be present in the Tarfaya.
distribution are between 4.7 and 6.1 GPa (Figure 5e), and modal values of the stiff distribution are between 19.0 and 19.6 GPa.
The Tarfaya shale has a distribution dominated by a peak centered around a modal value of 5.0 GPa (Figure 5f), which was different from the other distributions observed. This probably corresponds to the stiffness of the organics within both the shale matrix and the fecal pellets, suggesting that these two phases of organic matter are geomechanically similar. There is a second small peak centered around 12.5 GPa, which may be the start of a bimodal distribution, although only represents 17.74% of the data.

Deconvolution
The deconvolution results indicate that the two-phase system cannot be described perfectly as a combination of two Normal or two gamma distributions (Figure 7). However, the Normal distribution deconvolution Note. The mean value for the normal distributions is shown along with the modal values for the Gamma distributions on each shale. a This is an unrealistic second phase in the Tarfaya due to being stiffer than the maximum stiffness in moduls; the curve itself as observed in Figure 7 has a frequency of <0.001. b There is no modal value in the Tarfaya due to the nonfit of the data. The Tarfaya single phase indicates the error associated with fitting a single normal and gamma distribution to the curve.
leads to lower RMSE than the gamma ( Figure 10; Table 4) and appears to fit the data better. Both the normal and gamma distributions identify two definitive peaks: a first peak with a modal value of between 5 and 13 GPa and a second with a modal value of between 16 and 22 GPa in the deconvolution constrained by nonoverlapping means (Figure 7). Table 4 and Figure 8 demonstrates that the gamma distribution has a lower modal value than the equivalent normal distribution for each of the fractions, which is most likely due to the different shape of these distributions. The differences in modal value for the gamma and normal distributions are greater in the softer phase (mean difference = 2.47 GPa) than the stiffer phase (mean difference = 0.79 GPa).

Analysis of the phases in
The mean RMSE associated with the use of the gamma distribution on the constrained data is greater than the RMSE associated with the normal distribution (Table 4). This average trends toward a similar value when the Tarfaya shale is removed (normal RMSE = 0.014; gamma RMSE = 0.018). Neither of the deconvolution methods successfully fit the Tarfaya shale to the same accuracy as the other shales; however, the fit is significantly less accurate using a gamma distribution leaving an RMSE of ∼0.29 (Figure 10; Table 4); this appears as a flat line with a frequency of zero in Figure 4. Minimization using one-phase deconvolution leads to errors in Gamma fit of the same order of magnitude of the other errors. The error associated with each other shale is relatively similar with no obvious other outliers.
Deconvolution has also been undertaken relaxing Constraint B (Equation 6), which indicates an identical result when deconvolving into two normal distributions (Figure 8; Table 4). The error associated with unconstrained means on a gamma-based deconvolution yields a poorer average RMSE but does allow fitting of the Tarfaya shale with a smaller error. The unconstrained gamma fitting also results in many values not deviating from the initial estimate of 20.0 GPa, which may be due to the model finding a local minima.

Discussion
The organic matter Young's moduli values found within this study have two defining characteristics: a range of 0.1 to 24 GPa, fitting within the limits of the RTESP-525 tip, and for most shales a degree of bimodality. This range of Young's moduli has been reported by various studies using AFM (Eliyahu et al., 2015;Emmanuel et al., 2016a;Goodarzi et al., 2017; using single observations of individual samples, wells, or formations. We find that this trend of bimodality and Young's moduli <30 GPa is ubiquitous among shales from different formations and maturities and reflects the distribution of organic matter in the shale matrix.

Depositional Environment
Analysis of organic matter Young's modulus according to depositional environment is relatively difficult in these samples and would be in many other shales, as the main organic component is AOM. The most terrestrially influenced shale (Bowland) within this study comprises ∼95% AOM, with all other shales composed of at least 98% AOM. The scanned regions of each shale in this study contained only AOM, and thus, it is difficult to ascertain the depositional environments from AFM scans alone. Vitrinite and inertinite macerals within these shales had a long axis range of approximately 10-55 μm, which would comprise an entire AFM scan if found. Thus, it may not be prudent to select visible macerals to undertake AFM measurements due to the capacity to misrepresent the organic population as a whole.
However, AOM can be subcharacterized based on its organic precursor (Tyson, 2012). The AOM within Barnett, Bowland, Eagle Ford, and Tarfaya has been mostly derived from a planktonic algal marine precursor, which appears to be relatively uniform in its modulus pattern, with a bimodal distribution centered around peaks of 4-9 and 14-20 GPa. This implies that although the depositional environment may be important, among AOM it is adequate to use Rock Eval pyrolysis measurements to describe organic matter Young's modulus, rather than more time intensive analytical techniques such as palynofacies analysis.
Although geochemically the AOM within each of the shales may be similar, there are subtle differences in Young's modulus distribution. When deconvolving the modulus data, it appears that the softer phase has a lower Young's modulus (between 5.5 and 8.5 GPa) in the marly, calcite-rich shales of the Tarfaya and Eagle Ford, than the equivalent peak (between 9.4 and 12.8 GPa) in the more siliceous Barnett and Bowland shales. This may be reflection of the different appearances of organic matter within these two different facies, as fecally derived organic matter appears less stiff than whispy organic lamellae. This could be due to the carbonate coating and early cementation of the fecal remains, reducing the amount of compaction related stiffening in the organic matter.

Maturity
The evolution of organic matter Young's modulus with maturity has been previously identified by AFM on individual shales (Emmanuel et al., 2016a;Khatibi et al., 2018). This evolution is exhibited in an increase in the maximum and modal values of Young's modulus (Emmanuel et al., 2016a). The data presented here suggest that this is ubiquitous across Type II marine kerogen rich shales. The lowest modal value and maximum value of Young's modulus within this study is found within the Tarfaya shale, which is thermally immature.
The modal values and maxima in Young's modulus increase through the oil window and are at their highest within the gas window, exhibited by the distributions of the Barnett TP Sims and Bowland samples, which have a maxima of approximately 23 and 24 GPa, respectively. These maximum values appear to reflect similar values to those found by Khatibi et al. (2018) but are lower than those in Emmanuel et al. (2016a). The transition from a soft dominant to a stiff dominant distribution occurs within the early/peak oil window, demonstrated within the scans of the Eagle Ford Cinco Saus and Barnett Lake Davis samples, which have an equal distribution and a Tmax in the early oil window (Figure 9).
Although the maximum value of organic matter Young's modulus increases with maturation, the increase is rather modest, reflecting an increase of only 4 GPa from early oil to peak gas window. The distribution of modulus values appears to represent the change through maturity more clearly, becoming increasingly dominated by a peak with a modal value of between 19 and 24 GPa with maturity after the early oil window (Figure 9).
Whereas previous interpretations suggest that this bimodal distribution and evolution are simply due to kerogen and bitumen, this may be too simplistic a description of the observations presented. Two other factors must be taken into consideration: organic porosity created by the migration of bitumen from the system as petroleum and the overall chemical change to the organic phase (bitumen and kerogen).
Nanoindentation studies on shales with maturity suggest that there is a reduction in effective Young's modulus with increased maturity (Shukla et al., 2013;Zargari et al., 2016) caused by the development of organic intraparticle pores from bitumen generation and migration. If this were to be the case, it would be expected that the proportion of higher modulus values and overall modulus values would decrease with maturity. This is not exhibited either within this study or in AFM literature (Emmanuel et al., 2016a(Emmanuel et al., , 2016bKhatibi et al., 2018). This important observation indicates that one should be cautious when comparing AFM and nanoindentation-derived mechanical measurements. It is most likely that the created pores within the organic matter are resolvable by the AFM, (>10 nm in diameter), but not by the nanoindentor tip (<1μm in diameter). The AFM therefore allows insight into how the load-bearing organic matter changes in modulus once pores are created, demonstrating a Young's modulus increase that counteracts any decrease in effective Young's modulus due to porespace generation with maturity. Although this concept has been proposed by Suwannasri et al. (2019) using data-driven modeling, this is the first instance of observations of the general increase in solid organic matter Young's modulus with maturity.
With organic porosity unlikely to be the cause of the increasing bimodality in AFM results, chemical transformations of the organic matter are likely to be the cause in the observed trends. The results of the solid-state NMR undertaken on the kerogens indicate that there is a sharp increase in aromatic structures with maturity ( Figure 4); this trend toward an ordered graphite-like composition away from an organic melange could be the driver for the increase in modal modulus. Statistical tests however between the two distributions in the form of student t tests and Kolmogorov-Smirnov tests indicate little correlation. This is perhaps not unsurprising as the Young's modulus plots do not readily deconvolve into two separate phases using the technique applied in section 3.6, which inhibits the use of statistical tests dependent on normal distributions. Second, the NMR spectra do not readily distinguish between bitumen and kerogen, which further complicates this interpretation, given the changes in bitumen composition during maturation (Liu et al., 2018). Though the relatively low abundance of bitumen in immature shales (1-2% of overall organic matter in Tarfaya Rock Eval results) and large proportion of Young's modulus values within the lower modulus phase indicate that attributing this peak to bitumen alone is incorrect.
It appears that the Young's modulus distributions and evolutions observed within Type II kerogen are heavily dependent on the kerogen composition at the sample interval. The increases in aromatic kerogen observed by NMR (Figure 4) coincide with the increase in the proportion of higher modulus phase ( Figure 5). Raman spectroscopy by Khatibi et al. (2018) on the Bakken shale indicates there may be a correlation between the increase in aromatic hydrocarbons and Young's modulus values, which is further substantiated by the results of the NMR to Young's modulus comparison here.
Hydrocarbon chain length must also play a key role in the Young's modulus distribution, with lower chain length compounds likely to have a lower Young's modulus, whereby it may be that short chain aromatic compounds may overlap with long chain aliphatic compounds. During the early stages of the oil window, most of the C-N and C-S bonds and some of the C-O or bonds are broken, and the resulting hydrocarbon side chains are lost (Tissot & Welte, 1978). This causes the kerogen to become significantly less enriched in aliphatic hydrocarbons by the middle oil window (Figure 4), with a large variability in chain length. Retained bitumen and hydrocarbons generated in the oil window have been shown to exist as a complex mix of aromatic and aliphatic compounds with variable chain lengths (Cornford et al., 1983). This complex mixture of bitumen and kerogen may manifest as the equal bimodal distributions in the Barnett Lake Davis and Eagle Ford Cinco Saus (Figure 4), which appears as a kerogen composition similar to Figure 10b. This most likely has an adsorbed bituminous phase, which may significantly alter the strength of the aromatic kerogen and cause the complex Young's modulus distributions in Figure 5. It is difficult to comment on the and gas window maturity (∼1.3%R o ) after the models of (Behar & Vandenbroucke, 1987;Vandenbroucke & Largeau, 2007). The length of aliphatic chains, number of NSO compounds and general heterogenity decrease with increasing thermal maturty, until all that remains are aromatic rings with short aliphatic side chains. These compositions are similar to those within the Tarfaya (immature), Barnett Lake Davis & Eagle Ford Cinco Saus (oil window) and Barnett TP Sims & Bowland (gas window). exact nature of the bitumen within each shale, due to the changes in bitumen composition with maturity; the likelihood is that the bitumen will decrease in aliphatic carbon with maturity (Craddock et al., 2015), becoming more aromatic, though probably to a lesser extent than the load bearing kerogen at the equivalent maturity.
Deconvolution of the peaks two phases supports this kerogen distribution, in which the model identified two peaks but was unable to fit a normal distribution to them. This is caused by overlap between the two phases, and in the softer phase a large standard deviation. This could be due in part to the complex melange of organic compounds and polymeric material within the kerogen and bitumen. While the kerogen becomes enriched in aromatic components with maturity, these are most likely connected to aliphatic structures with side chains and branches. This would form a complex distribution with little homogeneity but two peaks representing broadly aliphatic and aromatic rich compounds, which we see in Figure 4. This deconvolution technique has been used successfully on many other minerals and composites (Constantinides & Ulm, 2007;Vandamme & Ulm, 2009), but mainly with defined stoichiometry, which is not present in the organic matter.
Although we believe that the two phases here are mechanically interlinked, relaxing Constraint B (Equation 6) does not improve the fit in either the normal or gamma distributions. This suggests that the two phases are intrinsically mechanically different. The inability to deconvolve the Tarfaya may suggest that in an immature shale, the organic matter is too mixed to deconvolve and requires the breaking of hydrocarbon chains to clearly separate out aliphatic and aromatic carbon. The NMR results suggest that the Tarfaya 10.1029/2020JB020435 has very little aromatic carbon, which may be either be the small subpeak in the 12-14 GPa range or be lost within the large soft aliphatic peak.
During the gas window, the majority of long chain aliphatic components have either been lost or are cracked to short chain gas hydrocarbon compounds (<C 5 ) during catagenesis (Tissot & Welte, 1978). This leads to a more ordered, aromatic rich structure (Figure 10c), exhibited by the aromatic dominance in the NMR results (87% in Barnett TP Sims and 88% in Bowland) in Figure 4. Any bitumen present at this maturity would most likely exist as highly volatile condensate, unlikely to contribute significantly to organic matter Young's modulus, or as pyrobitumen. While we cannot rule out the contribution of pyrobitumen (generated as a residue from the in situ cracking of bitumen to less dense components), to the organic matter Young's modulus distribution of the samples in question, no pyrobituminous organic matter was identified or scanned during this study. It is most likely that although pyrobitumen is characterized geochemically with bitumen, due to their common provenance as products of during thermal cracking, geomechanically pyrobitumen is most akin to kerogen, due to the polymeric nature of its composition.
These findings provide further complexity to modeling of hydraulic fracture tip and homogenization (Goodarzi et al., 2016), as models can be adapted for a dependence on both depositional environment and maturity. By correlating organic geochemical properties with geomechanical properties for the first time, we suggest the modeling of previously expensive to measure mechanical properties through an adaptation of relatively inexpensive techniques such as Rock-Eval pyrolysis. This linkage between geochemistry and geomechanics could be readily improved by undertaking geochemical measurements and AFM scans on cuttings obtained from wellbores of both carbon capture and storage seals and hydraulic fracture targets.
This research in combination with the thermodynamics of hydrocarbon generation may offer valuable insight into the initial migration of hydrocarbons from the organic matter host. The amount of pore pressure created during maturation is fundamentally linked to the mechano-chemistry of the organic matter, whereby generation of significant hydrocarbons may cause a pore pressure unsustainable by the host organic matter. The extent to which the organic matter within a shale can accommodate increases in stresses is determined by the mechanical properties and is influenced heavily by Young's modulus. The results here indicate that by the onset of generation of significant quantities of hydrocarbons (early peak oil window), a bimodal distribution of Young's modulus is already established (Eagle Ford Cinco Saus and Barnett Lake Davis). Therefore, further investigation may be necessary to determine which of the two distributions plays the major role in the formations of microfractures, enhancing primary migration. Using the trends in the data within this study, this problem may be approached in a stochastic manner.

Conclusions
Shale geomechanical properties are a key factor in the economic viability of an unconventional shale gas target and play a major role in determining the column height of CO 2 that can be housed in a conventional reservoir. As such, knowledge of these properties is important. Much is known about the mechanical properties of common minerals present in shale. However, the organic matter which acts as the host for the hydrocarbons in place and a diffusive path of the CO 2 , is less well understood.
Analysis of the Young's modulus of organic matter in four prolific black shales, with subtly different marine depositional environments and thermal histories, was undertaken on a nanoscale using AFM QI™, a similar technique to the PeakForce QNM™ used in other studies (Eliyahu et al., 2015;Emmanuel et al., 2016aEmmanuel et al., , 2016b. The results of this analysis indicate a bimodal distribution of organic matter Young's modulus in all of the shales studied, with a maximum value of 24 GPa.
Differences in depositional environments between marine shales have relatively little impact on the organic matter Young's modulus. The Young's modulus of marine organic matter evolves with maturity, initially dominated by values <10 GPa centered on a peak of ∼5-6 GPa while immature, which is believed to represent the aliphatic-rich heterogenous mixture in immature kerogen. At the oil window, the distribution evolves to become equally distributed between a peak at 3-6 GPa and a second, tighter distribution around 19-24 GPa, believed to represent the more aromatic kerogen. With increasing maturity into the gas window, this peak around stiffer values increases in maximum stiffness up to 24 GPa and contributes to 70-80% of the organic matter distribution, reflecting the dominance of aromatic structures within the organic matter in the gas window. This allows modeling of the Young's modulus in a Type II kerogen as a series of bimodal distributions that trend from a soft-dominated to stiff-dominated distribution with increasing maturity.
The Young's modulus distribution of organic matter in most samples cannot be explained using normal Gaussian distributions, and so future research should attempt to determine the Young's modulus distribution of these complex organic phases. Future studies should also investigate different types of organic matter, that is, different assemblages of heterogenous and homogenous amorphous organic matter along with other major macerals and maceral types, for example, spores vitrinite and inertinite.
This modeling of other macerals, combined with future studies on the evolution of Type I kerogen with maturity, will assist in overall modeling of organic matter Young's modulus in both cap rocks and hydraulic fracture prospects by allowing further correlation between the trends well known as part of source rock geochemistry and trends that occur within source/cap rock geomechanics.