Time-lapse nanometre-scale 3D synchrotron imaging and image-based modelling of the response of shales to heating

The development of pore and fracture networks at the nano-scale as a response to heating can reveal coupled physical relationships relevant to several energy applications. A combination of time-lapse 3D imaging and finiteelement modelling (FEM) was performed on two typical thermally immature shale samples, Kimmeridge Clay and Akrabou shale, to investigate thermal response at the nm-scale for the first time. Samples were imaged using Transmission X-ray Microscopy (TXM) with a voxel resolution of 34 nm at the I13–2 beamline at Diamond Light source, UK. Images were taken after heating to temperatures of 20 C, 300 C, 350 C and 400 C. The initiation of nano-pores within individual minerals and organic matter particles were observed and quantified alongside the evolution from nano-pores to micro-fractures. The major expansion of pore-volume occurred between 300 and 350 C in both samples, with the pores elongating rapidly along the organic-rich bedding. The internal pressures induced by organic matter transformation influenced the development of microfractures. Mechanical properties and strain distributions within these two samples were modelled under a range of axial stresses using FEM. The results show that the overall stiffness of the shale reduced during heating, despite organic matter becoming stiffer. The varied roles of ductile (e.g., clay minerals, organic matter) and brittle materials (e.g., calcite, pyrite) within the rock matrix are also modelled and discussed. The configurations of organic matter, mineral components, porosity and connectivity impact elastic deformation during shale pyrolysis. This work extends our understanding of dynamic coupled processes of microstructure and elastic deformation in shales to the nm-scale, which also has applications to other subsurface energy systems such as carbon sequestration, geothermal and nuclear waste disposal.


Introduction
Shale, a fine-grained sedimentary rock with a complex microstructure, plays an important role in energy and environmental systems. It can be the source and storage for natural gas (Braun and Rothman, 1975;Kobchenko et al., 2011;Chen and Xiao, 2014), can act as the cap rock for subsurface hydrogen storage or carbon sequestration (Lahann et al., 2013;Rohmer et al., 2016;Jia et al., 2019), and could also act as host rock for nuclear waste disposal (Neuzil, 2013;Liu, 2014;Charlet et al., 2017). The responses of shale to elevated temperatures are crucial to understand as they provide direct evidence to track burial history (Tourtelot, 1979), impact gas and geothermal recovery (Gehlin and Hellström, 2003), and have an impact on storage potentials and environmental impacts on nuclear waste disposal (McCarthy et al., 1978).
A range of coupled processes occur during the heating of shale, including hydrocarbon generation, multiphase flow and mineral reactions (Saif et al., 2019). It has been shown that thermal changes in shales can modify the pore structure, which consists of micro-fractures, mineral matrix pores, and organic matter pores (Loucks et al., 2009;Curtis et al., 2012;Tiwari et al., 2013;Ma et al., 2017). These pores range from a few microns to less than 1 nm, and control producible pore volumes and fluid flow rates (Nelson, 2009;Loucks et al., 2012;Ma et al., 2019). Recently, X-ray computed tomography (XCT), a nondestructive imaging technology, has been used to provide 3D structural information on shales from millimetre to sub-micron scale (Landis and Keane, 2010;Ma et al., 2016aMa et al., , 2017. The non-destructive nature of the XCT technique and ability for 3D reconstruction are beneficial for resolving spatial connectivity and quantifying geometric information in shales (Curtis et al., 2010;Ma et al., 2018Ma et al., , 2019. In-situ (Pilz et al., 2017) or ex-situ (Saif et al., 2017) imaging can be performed for characterisation of the dynamic pyrolysis process while heating shales to temperatures over 300 • C. Tiwari et al. (2013) reported the impact of temperature on shale pyrolysis based on XCT imaging at 42 μm voxel resolution. They found that larger pores were produced at higher temperatures, and the distribution of hydrocarbon products was likely to depend on the development of channels during pyrolysis. Saif et al. (2017) characterised the petrophysical properties of the Green River oil shale by ex-situ imaging of pyrolysis (380-420 • C) combining micro-CT (0.8 μm voxel) with Focused Ion Beam Scanning Electron Microscopy (FIB-SEM). They captured fractures nucleating at 380 • C and expanding in volume at 400 • C. Pilz et al. (2017) used 4D synchrotron X-ray tomography to quantify fracture propagation in a Kimmeridge Clay sample during heating, at temperatures between 370 and 380 • C. They documented the complex relationship between clay matrix, principal strains and fracture nucleation, and calculated the coefficients of thermal expansion based on full field Digital Volume Correlation (DVC). Panahi et al. (2018) investigated the effects of heating rate and confinement on fracture development of the Green River Shale by in-situ dynamic XCT imaging with a voxel size of 0.56-2.8 μm. Using 2D digital image correlation (DIC), they concluded that final temperature significantly influenced fracturing, and flake-like kerogen patches facilitated the development of inclined fractures. Saif et al. (2019) performed an insitu dynamic investigation of the Green River shale between 20 • C and 600 • C using the TOMCAT beamline X-ray micro-tomography with a voxel size of 1.63 μm. The results of DVC analysis showed that both the hydrocarbon generation and dynamics of the local strain impacted the expansion and closure of microfractures. However, these previous studies have focused more on the relationship between fracture evolution and shale pyrolysis and have yet to image nano-scale dynamics. The deformation and influence of other minerals are difficult to capture due to resolution limitations. The integrated process of organic matter transformation, mineral alteration and microstructure variation is thus sparsely explored.
With the advancement of X-ray optics and mechanics, Transmission X-ray Microscopy (TXM) has made tomography possible at nano-scale spatial resolutions (Andrade et al., 2016). A combination of condenser and Fresnel zone plate objective lens makes the TXM capable of magnifying radiographs and probing internal structure nondestructively (Gorelick et al., 2011;Anderson et al., 2020). Currently, TXM has been used in shale heterogeneity and anisotropy analysis (Sun et al., 2018), pore network characterisation , multimodal imaging (Anderson et al., 2020), and extraction of shale-gas during hydraulic fracturing (Kiss et al., 2015). TXM is uniquely well-suited for characterising the evolution of rock matrices including small organic matter particles and clay minerals, thereby extending the studies on thermal responses of shale to nm-scale.
Elasticity describes the stiffness of a rock, which significantly impacts fracture initiation and propagation. The commonly used DIC or DVC analyses in 3D imaging are powerful techniques for quantifying strain and displacements on fractures (Pilz et al., 2017;Panahi et al., 2018). The recent development of digital rock physics enables a combination of advanced imaging techniques and numerical simulation, where multiple simulations can run on the same 3D images to explore a series of petrophysical properties (Sakellariou et al., 2007;Dvorkin et al., 2008). Finite-element modelling (FEM) has been proved to be applicable in elastic prediction for 3D digital core images (Arns et al., 2002;Shulakova et al., 2013;Faisal et al., 2019). FEM is especially suitable for simulating the elastic behaviour of materials with multiple components featuring contrasting properties (Garboczi, 1998;Roberts and Garboczi, 2000;Sifakis and Barbic, 2012). Multi-scale X-ray imaging can provide a true 3D model with information about structure and compositions of the specimen for FEM simulation ((Shulakova et al., 2017)). Therefore, a combination of high-resolution non-destructive imaging and FEM can be used to develop a model to accurately reflect the heterogeneous rock structure. This is a significant improvement over using empirical assumptions, which are often forced to attribute deviations between experimental and numerically determined elastic moduli to micropores and intermediate phases that are smaller than the image voxel resolution (Saenger et al., 2011;Soulaine et al., 2016;Zhao et al., 2019). Despite the mature application of FEM in research related to microstructure (Chawla et al., 2006;Bieniaś et al., 2012;Ghosh and Chakraborty, 2013), there are limited studies investigating the relationship between the overall stiffness of shale samples and the dynamic evolution of microstructure and elasticity within individual components.
The research presented here aims to investigate the thermal response of shale at a ultra-high resolution, including the evolution of microstructure, pore-fracture network and elastic properties from room temperature to 400 • C. It utilises a state-of-the-art synchrotron based TXM imaging technique to provide time-lapse images of microstructural changes induced by changing thermal conditions at the highest resolution employed to our knowledge. The evolution of pores and fractures, together with surrounding minerals and organic matter, are studied for the first time at nanoscale in 3D and over time. The contributions of organic matter, pore structure, mineral matrix, connectivity and highmodulus pyrite to elastic properties and strain distributions at different temperatures are then discussed through coupling with FEM. The integrated analysis of high-resolution 3D imaging and simulated elastic deformation is promising and can be applied to wider materials systems to investigate time-dependent dynamic changes.

Sample selection and preparation
Two organic-rich and thermally immature shale samples were collected from the Akrabou formation in the Errachidia Basin, Morocco (Sample E) and the Kimmeridge Clay Formation in the Wessex Basin, the UK (Sample K). The former is a calcite-rich shale with a total organic carbon (TOC) of 16% and the latter is clay-mineral dominated with a TOC of 15%. These samples were chosen to characterize the response to heating in mineralogically-contrasting organic-rich shales.
Mineral components were determined and quantified using X-ray Diffraction (Bruker D8Advance XRD, Billerica, USA) with an uncertainty of ±1%, and are listed in Table 1.
Samples were cut into 50 μm cylinders on 1 mm diameter discs using a laser cutter 3D-Micromac microPREP system equipped with a 532 nm Diode-Pumped Solid State picosecond laser) in the Manchester X-ray imaging Facility (MXIF). and the laser power is only 0.1 W during cutting so the heat effects can be ignored. Samples were then mounted on sample holders using high-temperature ceramic adhesion agent.

Time-lapse imaging with transmission X-ray microscopy (TXM)
The imaging experiment was conducted using phase contrast TXM at I13-2 beamline at Diamond Light Source (DLS) (Storm et al., 2020). Phase contrast imaging using Zernike phase rings can significantly enhance the signal of weakly absorbing samples and improve statistics in the reconstruction (Storm et al., 2018). The energy was 12 KeV and the exposure time was 1.8 s. 1498 images were acquired per tomogram at 34 nm voxel edge-length.

Ex-situ heating procedure
Previous studies have suggested the major thermal response in shales occurs at temperatures above 300 • C (Pilz et al., 2017). Therefore, due to the time-limited nature of synchrotron beamtime, tomographic images were not collected between 20 • C and 300 • C. Furthermore, owing to the experimental constraints of space and stability, in-situ scans during heating were not possible. The working distance of the X-ray optics is on the order of 10s of millimetres and limits the available space for the installation of a furnace. In addition, the alignment of the X-ray optics is very sensitive and drifts must be limited to below 100 nm to preserve the alignment of the X-ray optics for the phase contrast imaging. The stability argument also holds for the sample. In order to reach resolutions on the order of 100 nm, drifts must be kept to below this number. This means that a nearly perfect thermal insulation between sample and the mounting pin would be required and the temperature of the sample and mounting pin would must not change for a scan. This means that a constant temperature ramp would not be possible and after each temperature change, the sample would need to equilibrate again completely. This is very time consuming and considering the limited amount of available beamtime, offline heating and cycling allowed to measure several samples in parallel.
Samples were heated in an oven to specific temperatures (300 • C, 350 • C and 400 • C) at a heating rate of 5 • C/min. The temperature was then held for 1 h before cooling at the same rate of 5 • C/min. A tomogram was collected at room temperature for each sample as a reference before pyrolysis. A tomogram was then collected after heating to each temperature point, and subsequent cooling to room temperature. Sample K broke into pieces during heating to 400 • C owing to the abundant fractures generated inside the sample, therefore the image at 400 • C was not able to be collected.

Image processing and quantification
The tomography images were reconstructed using the Savu processing pipeline (Wadeson and Basham, 2016). Ring artifacts were removed using a method developed by Vo et al. (2018), and the Gridrec algorithm was used for the tomographic reconstruction (Dowd et al., 1999;Marone and Stampanoni, 2012).
The reconstructed slices for different temperature points were registered first using Avizo software (Version 9.3, Thermo Fisher, USA) to align the centre and allow investigation of the evolution between tomograms. A non-local means filter (spatial standard deviation 3, intensity standard deviation 0.2, search window 10 pixels, local neighbourhood 3 pixels) was then applied to reduce noise (Fig. 1). Different phases (e.g. pore, minerals, and organic matter) were segmented using a watershed algorithm. Volume fractions of different phases are displayed in Table 2. Segmentation errors were calculated by left-and rightdeviating 1 minimum unit in grayscale thresholding of each phase and the results are shown in Table 2. Quantification of different phases, such as equivalent diameter, local thickness and orientation (Z direction and XY plane), were performed using Avizo. Specifically, equivalent diameter is the diameter of a sphere with the same volume (Jennings and Parslow, 1988); local thickness is the diameter of the largest sphere  which lies within an object (Hildebrand and Ruegsegger, 1997); orientation indicates the direction of the longest axis. Pores with at least one common vertex are considered connected (axially connected along Z direction). Shape factor is defined in Avizo as A 3 /(36πV 2 ), where A is the surface area of the pore and V is its volume (Youssef et al., 2007).

Finite element simulation
Elastic properties of shale sample at each temperature were simulated by Abaqus FEA (Version 6.14, SIMULIA, USA). TXM images were pre-processed and finite-element meshes were constructed before elastic simulation. The original meshes were generated in Avizo and were then simplified and remeshed to increase the efficiency of the simulation (Shulakova et al., 2013). The best isotropic (BI) algorithm (Zilske et al., 2007) with a 50% reduction coefficient was applied to triangulated surfaces of the meshes to perform the FEM. Edge flipping, edge collapsing and vertex translation were utilized to achieve a good tetrahedral grid generation. Uniaxial deformation was then simulated by applying a normal pressure to a certain plane (XY plane in this work). Elements of samples were assumed to be isotropic and linear elastic, which is a common assumption of the FEM method (Shulakova et al., 2013;Faisal et al., 2017). The ratio of mean stress and mean strain was calculated to obtain the elastic modulus.
Quantitative analyses were performed at different scales to determine the minimum representative size of typical features, using similar procedures as (Ma et al., 2016a) and . The minimum trialled subvolume size was 50 voxels in three stochastic centroids, and was increased by 50 voxels in each dimension subsequently (Fig. S1). The deviation of volume fraction of different phases compared to the entire bulk sample was calculated. The representative subvolume was defined as that where a relatively stable variation of elastic modulus (within 5 GPa), and a small deviation in volume fraction (within 5%) is reached. This same subvolume was selected to research the strain evolution with increasing temperature. According to Hui et al. (2019), the formation pressures of shale reservoirs are generally in the 15-50 MPa range. Some overpressured shale formations display very high formation pressures up to ~74 MPa, e.g. shales in the Dongpu Depression of the Bohai Bay Basin (Luo et al., 2016). These pressure ranges were utilized to model the strain distribution under different pressure conditions. It should be stressed that shale is generally not isotropic, but deviations here are acceptable in these nano-to micro-scale models owing to difficulties in the measurements of thermal properties at nanoscale (Hopkins, 2013).

Results
Four phases were identified for quantification: pores (and fractures), clay matrix with organic matter, granular minerals (e.g. calcite, quartz) and heavy minerals (e.g. pyrite).

Evolution of pore-fracture system
In sample E, fractures developed, with porosity increasing by 3.7% from 20 • C to 300 • C, and these continued to grow at 350 • C and 400 • C, resulting in a final porosity of 12.1%. Small pores were discrete at 20 • C ( Fig. 2 A1), and 60% of the pore volume was occupied by pores with an equivalent diameter smaller than 2 μm (Fig. 3 A1). At 300 • C, the pore equivalent diameter was in the range of 0.2-5.8 μm, and the pores that had an equivalent diameter < 2 μm developed rapidly in volume compared with the original sample ( Fig. 3 A1 develop between pores after heating to 300 • C. The percentage of connected pores are increased to 5.4% at 350 • C (Fig. 2 C2) and the fracture widening was also observed after heating to 350 • C (Fig. 3 A2). A further increase in temperature enhanced the connectivity of pores which propagated ultimately into a large interconnected network occupying a volume fraction of 9.2% at 400 • C (Fig. 2 D2). At this temperature the pore equivalent diameter varied from 0.4 μm to 15.9 μm (Fig. 3 A1). The largest connected pore equivalent diameter increased to 13.4 μm at 350 • C with a local thickness of 18.9 μm, and to 16.0 μm at 400 • C with a larger local thickness of 21.6 μm. Pores and microfractures mainly developed within organic matter and at the interface of organic matter and minerals (Fig. 4). A temporary closure of microfractures was also observed after 300 • C, which is likely due to the effect of localised stress ( Fig. 4 A1-A4). By contrast, there was no well-connected porosity formed within sample K between 20 • C and 350 • C although the porosity increased from 0.9% to 4.1% over this temperature range (Fig. 2 E1-G1). Between heating from 350 • C to 400 • C sample K broke into pieces, making imaging at 400 • C impossible. This is likely due to proliferation of connected fractures through the whole sample during this heating step. At 300 • C, the pore equivalent diameter was in the range of 0.3-1.7 μm, increasing to 0.4 μm to 6.0 μm at 400 • C (Fig. 3 B1). The local thickness of most microfractures showed a trend of reducing first and then slightly increasing during heating. This may reflect the variation of internal pressures during this process. Microfractures in sample K underwent propagation and compression of their smallest dimension while heating to 300 • C, before gradually widening as temperature increased further (Fig. 3 B2). Mineral deformation was obvious during heating and cooling in sample K (Fig. 4 B1-B3). Due to the formation and elongation of new pores, it is difficult to track the development of the geometries and locations of some pores. Compression and collapse occurred when the compaction caused by mineral deformation exceeded the pore expansion during heating, and some of these pores deformed further during the cooling (Fig. 4 B1-B3). The small pore indicated by red arrows in Fig. 4 B1-B3 also appears to display temporary compression driven by expansion in nearby materials. However, the pore volume was actually found to increase slightly, measuring about 2.159, 2.162 and 2.200 μm 3 at 20 • C, 300 • C and 350 • C, respectively ( Fig. 4 B1-B3).
Changes of fracture orientations were analysed relative to the plane parallel to the Z-axis (0-90 o ) and the XY plane (0-360 o ). There is a clustering of fractures oriented at around 75 o to the Z direction (Fig. 5). The orientations of pores are diverse in the XY plane for sample K, whereas pores within sample E display an alignment at 320-337.5 o which is weakened during heating. The orientation distribution further demonstrates that major microfractures preferentially propagated along the bedding plane.

Evolution of organic matter distribution
Due to the high clay mineral content and the limitation of distinguishing the grayscale boundary between clay minerals and organic matter in sample K, both clays and organic matter are defined as the clay phase. In order to facilitate subsequent elastic simulations, we considered the phase with grayscale between pores and minerals to be mainly clay minerals mixed with organic matter in sample K. Therefore, changes related to organic matter alone for sample K are not discussed. There is a large and well-connected organic matter particle existing in sample E (Fig. 6). The connectivity of the large organic particle was damaged during pyrolysis. The volume percentage of connected organic matter of sample E fell from 37.2% at 20 • C to 28.3% at 400 • C. The scatter diagrams in Fig. 6 indicate that in the samples' original state, disconnected organic particles are near spherical and elliptic in shape. The organic matter regions of sample E with an equivalent diameter smaller than 0.8 μm are near elliptic, whereas scattered organic matter regions with a larger equivalent diameter (> 0.8 μm) became elongated in shape as temperature increased.

Evolution of elastic properties
Based on XRD measurements (Table 1), sample K has high calcite and clay mineral (kaolinite and illite) content accounting for 46% and 22% weight respectively. By contrast, sample E is much more heavily dominated by calcite, which has a weight percentage over 85%. Four phases of sample E were used in our simulations, pores, organic matter, calcite and pyrite. We ignored the impact of clay minerals on stiffness because of the very limited clay content (Table 1). The Young's modulus and Poisson's ratio of the different minerals are required in the elastic simulation and were sourced from the scientific literature (Table S1). It should be stressed that several researchers have observed a stiffening of organic matter with increasing thermal maturity (Emmanuel et al., 2016;Li et al., 2018). Thus, Young's modulus values of 10 GPa, 13 GPa, 15 GPa, and 16 GPa were selected for organic matter at 20 • C, 300 • C, 350 • C and 400 • C respectively, to mimic this stiffening process. Pores were assumed to be gas-saturated with a modulus of 0.2 GPa. Calcite and pyrite were assigned Young's modulus values of 56 GPa and 300 GPa, respectively (Eliyahu et al., 2015;Emmanuel et al., 2016). The input Poisson's ratio of sample E was 0.3 for calcite, organic matter and pore phases (Shulakova et al., 2013;Devarapalli et al., 2017), and 0.15 for pyrite (Schön, 2015). By contrast, an "average" elastic modulus of 60 GPa and Poisson's ratio of 0.3 were used for the mineral phase in sample K considering the relatively high content of carbonates than quartz ( Table 1). The clay phase was assigned a Young's modulus of 20 GPa.
Within these subvolumes, volume percentages of each phase remained within 5% of those found for the whole sample, and relatively stable moduli were observed for subvolumes larger than 150 3 voxels (Fig. S1). Since the distribution of phases in sample E is more uniform than that in sample K, a larger subvolume of sample K was used to generate a more accurate result. Considering the necessary model size Fig. 4. 2D cross-sections of TXM images. (A1)-(A4) display microfracture variation of sample E, where red arrows show the closure of a microfracture due to effects of local stress. (B1)-(B3) show the phenomenon that pores develop (yellow arrow) and reduce (red arrow) during heating of sample K due to local difference in the coupling process of expansion and compression. Although the small pore (red arrow region) appears to be closed at the same slice position, the whole volume displays a bit increase. (OM: organic matter; Cal: calcite; Py: pyrite). (For interpretation of the references to colours, the reader is referred to the web version of this article). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) and computation time, 10 cubes with 150 3 voxels of sample E and 200 3 voxels of sample K at each temperature were extracted randomly to investigate the impacts of compositions and pore structure on elastic deformation during shale pyrolysis. Shale stiffness was found to reduce during unconfined heating, despite organic matter becaming stiffer in sample E (Fig. 7). The elastic modulus of sample E ranges from 47.0 GPa to 65.1 GPa with a mean value of 56.8 GPa at 20 • C, and from 7.6 GPa to 59.6 GPa at 300 • C (mean 32.0 GPa). Due to shale pyrolysis and development of a connected pore network, the elastic modulus lies in the 3.6-34.8 range GPa (mean 12.3 GPa) at 350 • C, but falls to 2.9-24.5 GPa (mean 7.3 GPa) at 400 • C (Table S2). The average unheated Young's modulus of sample K is lower than that of sample E due to a relatively high clay content (Fig. 7, Table S3), which falls from 35.5 GPa at 20 • C to 19.1 GPa at 300 • C and finally to 13.8 GPa at 350 • C. Above 300 • C, the trend of modulus reduction in sample K slows down.

Evolution of strain distribution with increasing vertical stress
The strain distributions within the selected sub-volumes were mapped under different simulated vertical stresses to investigate strain at subsurface conditions. Fig. 8 A1-D2 show the 3D strain evolution within a subvolume of sample E at constant axial stress (40 MPa) but varying pyrolysis temperatures. The highest strain is concentrated in the newly developed pore-fractures, with smaller strains accumulating within the organic matter regions. In sample K, the fractured areas between clay and mineral phases display larger strains (Fig. 8 A3-C4). Under increasing axial stresses (0.1 to 70 MPa), large strains accumulate in the fractured areas between clay/organic matter and minerals, and then spread to these dark phases (Fig. 9). These pores and fractures close up under higher confining pressures potentially inducing pore collapse in high strain regions in the subsurface.

Advantages of time-lapse TXM in thermal response studies of shales
Time-lapse TXM extends observations of the thermal response of shales to the nm-scale. To our knowledge, this is the first time the thermal response has been studied for shales at 34 nm voxel size. Micropores with an equivalent diameter smaller than 1 μm occupy 24.1% and 29.5% of porosity in sample E and K respectively. This nm-scale imaging therefore allows investigation into the influence of a much larger proportion of the true porosity of the material than those at μm scale. Compared with micro-scale studies, the nm-scale observation in this study can characterize and quantify the evolution of pores within individual minerals or organic matter particles in 3D (Fig. 4). Local stresses within organic matter particles can cause nano-pore initiation before accumulating toward sudden failure in the surrounding minerals, driving the formation of micro-scale fractures. The application of TXM in this study allows for the local stress and nano-pores to be visualised before the presence of fractures is observed. This nano-scale imaging provides important evidence of the earlier stages of deformation for 20   The increase in porosity that we observe during heating is primarily attributed to the transformation of organic matter into mobile hydrocarbons. During the transformation process, fluid pressure will increase within the organic matter. When the internal vapor pressure exceeds the local mechanical strength of grain boundaries, microfractures will elongate and grow rapidly (Kobchenko et al., 2011;Saif et al., 2019). The expansion and compression caused by accumulation of internal pressure can transform the shape of organic matter, which may change the contact relationship between particles, potentially affect the overall stiffness (Fig. 8). All of these sub-micron perturbations may be overlooked in common micro-scale X-ray imaging. Therefore, in order to reveal local changes it is necessary to use nm-scale technique as a supplement of micro-CT. Similar elastic upscaling work combining XCT and FEM have been performed on sandstone and carbonate rock by Shulakova et al. (2013) and Faisal et al. (2017). Each of these studies considered two phases in their simulation due to resolution limitations, i.e., single mineral phases and pore. We are able to segment more phases based on the high-resolution TXM images, improving the accuracy of simulation. In order to assess the influence of voxel size on elastic simulation results, we re-sampled the images with decreasing voxel resolution. A 5 μm × 5 μm × 5 μm region of interest was selected to perform re-sample at voxel sizes of 34 nm, 68 nm and 102 nm (Fig. S2A1-C1). The dark phase contains pores, clay minerals and organic matter, whereas other minerals are assigned as the bright phase. Decreasing resolution resulted in mosaic edges of dark phase, leading to a relatively low gray value of the edge and a reduced volume fraction of the dark phase (Fig. S2A2-C2). The increase of segmented bright minerals lead to an increasing elastic modulus, which was found to be 52.9 GPa, 54.4 GPa and 55.9 GPa at 34 nm, 68 nm and 102 nm voxel size respectively. This demonstrates that high-resolution TXM imaging is an advantage in improving the accuracy of simulations based upon imagebased modelling.

Contribution of compositions and pore structure to elastic properties
Shale brittleness reflects its stability and ease of fracture (Rickman et al., 2008;Carmichael, 2017;Huo et al., 2018). Both Young's modulus and Poisson's ratio are important parameters reflecting the ability to initiate and propagate a fracture under stress, and influence the rock brittleness (Rickman et al., 2008). The higher the Young's modulus, the less elastic deformation will occur under stress. Since Young's modulus is only material-dependent, the petrophysical properties of shale play a key role in stiffness. Relationships between the composition, porosity and simulated elastic modulus of sample E are displayed in scatter diagrams (Fig. 10). Pores are undeveloped at 20 • C for sample E. Only 4 of the 10 extracted cubes at 20 • C can be quantified for porosity, and fall between 0.1 and 1.4% (Table S2). When pores are undeveloped, the influence of calcite and organic matter on elastic properties is pronounced (Fig. 10). As the porosity and pore connectivity increased during shale pyrolysis, the development of microfractures played a more significant role than minerals and organic matter. The contribution of pyrite cannot be ignored since it has the highest Young's modulus of minerals. At 20 • C the volume fractions of organic matter, calcite and pyrite are 44.0%, 51.4% and 4.6% respectively in cube 1, and 45.0%, 55.0% and 0% respectively in cube 3. They have a similar volume percentage of dark phase and bright minerals. However, cube 1 displays a lower elastic modulus (47.0 GPa) even with more pyrite than cube 3 (53.1 GPa). The increase in total connectivity of calcite and organic matter also has a positive effect on the overall stiffness of sample E as temperature increases (Fig. 11). Pyrites in sample E are predominantly euhedral crystals instead of framboids. The organic matter is enriched near the pyrite and the stiffness of sample tends to decrease with a high content of organic matter since organic matter usually exhibits lower Young's modulus compared with other minerals (Huo et al., 2018). It can be inferred that the positive influence of mineral content on the elastic modulus at high temperature is mainly controlled by mineral connectivity. The apparent connectivity varies depending on the image resolution and segmentation, and so this is important to account for when comparing elastic moduli inferred from studies using differing resolutions. In addition, porosity displays a negative impact on shale stiffness of sample E (Fig. 10). The configuration of mineral compositions, porosity, fracture development, organic matter and mineral connectivity all influence elastic properties.
Pyrite has a significantly higher gray value than other components, and X-ray microscopy is therefore extremely well suited for pyrite identification. Discrete pyrite crystals in sample E display a slightly Fig. 10. Relationship between elastic modulus and phase content and porosity of sample E for each experimental temperatures (Y axis is elastic modulus (GPa) and X axis is volume fraction of each phase; ones without trend lines mean that the trend lines are difficult to predict due to low correlation coefficient where R 2 is lower than 0.2). negative influence on the elastic modulus for subvolumes where the pyrite content exceeds 1%. When the pyrite content is less than 1%, its impact on shale stiffness is limited regardless of maturity (Fig. 10). To investigate the influence of pyrite content further, we selected several subvolumes at each temperature and replaced the pyrite phase with calcite phase to remesh and perform FEM again. This replacement process simply involves incorporating the grayscale of pyrite into the mineral phase. The results demonstrate that the inclusion of pyrite can increase the overall stiffness to a certain extent in a given subvolume, but this increase is quite limited in most cases (Fig. 12). Although pyrite has a high Young's modulus, it does not significantly enhance the stiffness of sample E in our models. Stiffness is not only content-depended, but shape or connectivity also contribute. According to the simulation results (Table S2), the "soft phases" (organic matter and pore) are a pronounced control of stiffness in sample E. High content of "soft phases" reduces the stiffness of shale significantly even with the presence of high-modulus pyrite. Pyrite may increase shale brittleness due to its high elastic modulus, but its discrete distribution may also damage the connectivity of mineral matrix. Furthermore, the development of pyrite can indirectly promote the fracturing of shale reservoirs. Pyrite is primarily developed in a reducing environment, which is favourable for the accumulation and preservation of organic matter (Demaison and Moore, 1980;Wilkin et al., 1996;Wignall and Newton, 1998). Pyrite can accelerate the catalytic hydrocarbon generation, and the development of pores can be further facilitated after hydrocarbon expulsion (Mango et al., 1994;Wang et al., 2014;Ma et al., 2016b). A complex natural pore network is conducive to shale gas production. It should be stressed that the maximum pyrite content is 5.07% in our simulations, with most being only ~1%, as determined from the experimental results. We are unable to test the impact of higher pyrite content owing to the actual pyrite distribution in the samples. Furthermore, framboidal pyrite is a common phase in many shales, whereas pyrites in these sample are not framboids, i.e., they are crystalline without pores. Further research is needed to obtain a more comprehensive understanding of the influence of different pyrite structures on the elasticity of shales.
It is difficult to discuss the impact of organic matter alone in sample K dues to its incorporation into the "clay phase". The elastic modulus displays a relatively weak downward trend with increasing porosity, whereas increasing bright minerals potentially enhance the stiffness (Fig. 13). The relationship between phases and modulus is less clear than in sample E. The presence of pyrite is also does not appear to cause increased stiffness (Fig. 12). It can be concluded that the interactions between different phases in sample K are more complex, and cannot be accurately assessed from a single factor. Higher resolution imaging and segmentation of more phases is necessary to reveal the controlling factors of stiffness evolution in such a sample. Combined with the uneven distribution of strain in the same phase (Fig. 8 A3-C4), the strain variation also not only depends on the content of different phases but also on the dynamic evolution of shape, connectivity and coupled relationship between each component and its neighbours.
Since it is difficult to distinguish clay minerals and organic matter, we further assigned the dark gray voxels as clay phase and organic matter for sample E and K respectively to find uncertainties in the elastic properties (Table 3). The uncertainty was calculated by comparing with the original modelling results. It can be seen that the elastic properties of sample E at 20 • C are heavily influenced by the stiffness of dark phase due to the very low porosity. As pores developed with increasing temperature, the uncertainty decreases which shows that stiffness disparity due to the phase change is weakened. Elastic deformation therefore appears more sensitive to how much pore space is created during heating.

Effects of organic matter on elastic deformation
Organic matter dominates hydrocarbon generation, accumulation and transportation processes (Jarvie et al., 2007). The enrichment of organic matter is controlled by a variety of factors, such as phytoplankton blooms, nutrient supply and recycling, and anoxic events (MacIlvaine and Ross, 1973;Dimberline et al., 1990;Sageman et al., 2003;Schieber et al., 2010). Bitumen, Tasmanites, inertinite, telalginite, Fig. 11. Relationship between elastic modulus and total volume percentage of connected compositions of sample E for each experimental temperature. etc., are common organic macerals in shale, which indicate different origins. Abarghani et al. (2020) found that two adjacent bitumen particles within the same shale specimen showed distinct chemical structures, demonstrating separate parent macerals. Yang et al. (2017) proposed that a relatively high mechanical stiffness is often displayed in organic macerals with abundant aromatic carbon. In addition to differences in the chemical structure of organic matter caused by different sources of parent material, organic components would also change after undergoing the transformation of hydrocarbon generation during burial (Okiongbo et al., 2005). Emmanuel et al. (2017) investigated modulus changes of kerogen and bitumen in shales during experimental heating by in-situ AFM measurements. They found that the kerogen remains at an almost stable stiffness, while bitumen becomes significantly softer with a modulus of 14.1 GPa falling to 3.9 GPa with heating from 25 • C to 150 • C. Although it is hard to distinguish different organic macerals in the X-ray microscopy image, it can be observed that some, but not all organic matter generated large pores and fractures, which may indicate distinct organic macerals (Fig. S3). Particularly in sample K, the high TOC (15%) lead to large fractures developing inside the samples during heating from 350 • C to 400 • C, and at higher temperatures led to the sample breaking into pieces.
To research the influence of hardening or softening of organic matter on elastic deformation, one cube at each temperature (4 cubes in total) was selected for sample E to perform elastic simulation. Separate simulations were performed with 5 GPa, 10 GPa, 13 GPa, 15 GPa and 16 GPa as the input Young's modulus of organic matter. Here, low maturity shale is defined as the sample at 20 • C, while high maturity shale corresponds to the same sample heated to high temperatures. It can be observed that the presence of stiffer or softer organic matter has little impact on elastic deformation in high maturity shale (Fig. 14). However, low maturity shale is sensitive to the effects of organic matter stiffness. The source and maturity are two main factors influencing elastic properties of organic matter. Therefore, we infer that the source of organic matter-induced stiffness differences in high thermal maturity shale may not be significant when considering elastic deformation since the developed fractures contribute more prominently to the material stiffness. Conversely, organic-matter induced stiffness differences should be taken into consideration for low thermal shale. This result implies that when performing reservoir fracturing and reconstruction, stiffness changes induced by differences in organic chemical structure should be considered to maximise prediction accuracy.

Implications for energy applications
For the development of shale gas, high elastic modulus usually indicates a high stiffness. Brittle shale responds well to hydraulic fracturing for enhanced gas recovery (Rickman et al., 2008). In addition, heat generated by radioactive materials can affect sealing behaviour in shales and this is one of the major impacts that needs to be considered. Furthermore, geothermal or heat flows in a sedimentary basin may impact the mechanical behaviour of shales (Tester et al., 2006;Ghassemi, 2012). The high-yield Barnett Shale in the USA broadly benefits from high Young's modulus and low Poisson's ratio (Grigg, 2004). However, shale is often highly heterogeneous, even within the Barnett Shale, involving carbonate-rich, quartz-rich and clay-rich lithotypes. The results presented here suggest that the connectivity of organic matter and mineral phases in shale tends to dominate variations in stiffness. The content of brittle minerals and connectivity should therefore be considered as important factors when designing fracturing programs. For high maturity shales, pores, microfractures, connectivity of mineral matrix and organic matter are each important for determining stiffness. As maturity increases, hydrocarbon generation and expulsion lead to rapid pore development. The preferred fracture growth direction lies parallel to the bedding plane, which is in agreement with our previous studies (Ma et al., 2016a). Based on the simulation results under varying axial stresses, we speculate that if organic matter or minerals are well connected, it is easier to form large and continuous cracks at their interfaces during hydraulic fracturing. Additionally, the enhanced connectivity of pore networks in high maturity shale can provide larger pathways for gas migration and diffusion to the mineral matrix. During diagenesis, these fractures are likely to be filled with mineral cements which lead to a limited improvement of reservoir physical properties. However, such healed fractures can also act as zones of weakness, enhancing the effectiveness of the induced fractures during hydraulic fracturing (Bowker, 2007).

Limitations and outlook
It should be underlined that this work was performed in an open system and the images were taken at room temperature after cooling. Therefore, the microstructural response is determined at a steady state after heating and cooling cycles. The immediate reaction of rapid temperature changes are not considered in this study. In reality, the increasing of temperature from room temperature to 400 • C takes millions of years therefore the thermal response of rocks normally occurs in a relative steady confined subsurface system. Also, the grayscale boundary of clay minerals is difficult to be distinguished with organic matter especially in sample K. For the clay-rich shale, more realistic results that differ from this work may be drawn if different clay minerals could be reasonably identified from organic matter in a closed system with confining pressure. Smectite tends to be transformed to illite with increasing burial temperature in the subsurface (Eberl and Hower, 1976;Lynch, 1997). Especially in a closed system, the free silica yielded by clay mineral transformation will be removed, thus this process can create potential conditions for shale silicification (Thyberg et al., 2010;Metwally and Chesnokov, 2012). Changes in mineral composition and pore network affected by swelling and transformation of clay during heating should also be expected to influence shale stiffness. In addition, there is a large well-connected region of organic matter in our sample E. The shape and connectivity of the original organic matter will also influence the artificially matured experiment. Comparison experiments are still needed to obtain more comprehensive understanding of shale pyrolysis processes, such as shale with different distributions of organic matter and minerals, or in fluid systems with varied pressures. Characterisation results in this work largely reveal reasonable trends of microstructural and elastic variation during shale pyrolysis which may not represent the accurate values, because quantification and simulation depend strongly on how many phases can be reasonably identified and segmented in the current resolution. It is worth re-iterating that confining pressure is crucial in the shale heating process and some kinds of confined heating experiments will be necessary to fully understand these processes. Our ongoing works on various shale compositions and structures based on high resolution and multi-scale imaging techniques aim to explore these systems further. Multiscale mechano-chemical effects on rock stiffness caused by the geochemical reaction of fluid and rock under high temperature and pressure environments will also be investigated in future research. Widespread applications can be expected with the development of advanced high-resolution 3D imaging and modelling work, for example, we can monitor the in-situ behaviour of catalysis and track microscale characteristics of fuel cells (Sheppard et al., 2017;Singh et al., 2019). Finally, our approach can be applied to other materials and scales. Further studies into other materials should carefully consider the nature of their material's temperature-sensitivity, since the process of scanning after cooling may be of varying significance in different materials.

Conclusions
(1) This study presents the deformation of microstructure and elastic properties during experimental shale pyrolysis based on ex-situ 3D imaging at a resolution of 34 nm voxel size. Synchrotron based Transmission X-ray Microscopy is performed to investigate the thermal response of shales at the nm-scale for the first time, revealing the coupled relationships between compositions, microstructure and elastic deformation of shale.
(2) Nano-pores initially develop within organic matter and the mineral matrix (below 300 • C) and further develop into microfractures within the organic matter and at the interface of organic matter and minerals at higher temperatures (300-350 • C). The primary direction of pore-fracture propagation is along the kerogen-rich bedding with a maximum thickness of 21.6 μm. The connectivity of organic particles was reduced during heating. A temporary closure of some microfractures was observed owing to the effect of local stresses. (3) The elastic properties are affected by not only the compositions and pores, but also the connectivity of organic matter and mineral phases. The major reduction of elastic properties occurred below 350 • C for both samples. The high content of "soft phases" reduces the stiffness of shale significantly, even with the presence of high modulus blocky pyrite. The overall stiffness of shale is weakened during heating even though the organic matter becomes stiffer. We have shown that subtle changes in mineralogy and microstructure have an impact on mechanical properties of shale. (4) Under axial stress, the highest strain is concentrated in the newly developed pore-fractures, followed by organic matter zones. These pores and fractures may be closed under higher confining pressures in the subsurface and pore collapse may occur in the mineral or organic matter zones which have high strains. (5) We have demonstrated the importance of the spatial distribution of different phases at nm-scale on the thermal response of shales. The combination of high-resolution 3D imaging and deformation simulations can be applied to various energy systems to research time-dependent dynamic changes, such as geothermal reservoirs, catalysis and batteries.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.