4D Synchrotron X-ray Imaging of Grain Scale Deformation Mechanisms in a Seismogenic Gas Reservoir Sandstone During Axial Compaction

Abstract

Understanding the grain-scale processes leading to reservoir compaction during hydrocarbons production is crucial for enabling physics-based predictions of induced surface subsidence and seismicity hazards. However, typical laboratory experiments only allow for pre- and post-experimental microstructural investigation of deformation mechanisms. Using high-resolution time-lapse X-ray micro-tomography imaging (4D µCT) during triaxial deformation, the controlling grain-scale processes can be visualized through time and space at realistic subsurface conditions. We deformed a sample of Slochteren sandstone, the reservoir rock from the seismogenic Groningen gas field in the Netherlands. The sample was deformed beyond its yield point (axial strain > 15%) in triaxial compression at reservoir P–T-stress conditions (100 °C, 10 MPa pore pressure, 40 MPa effective confining pressure). A total of 50 three-dimensional µCT scans were obtained during deformation, at a spatial resolution of 6.5 µm. Time lapse imaging plus digital volume correlation (DVC) enabled identification of the grain-scale deformation mechanisms operating throughout the experiment, for the first time, both at small, reservoir-relevant strains (< 1%), and in the approach to brittle failure at strains > 10%. During small-strain deformation, the sample showed compaction through grain rearrangement accommodated by inter-granular slip and normal displacements across grain boundaries, in particular, by closure of open grain boundaries or compaction of inter-granular clay films. At intermediate and large strains (> 4%), grain fracturing and pore collapse were observed, leading to sample-scale brittle failure. These observations provide key input for developing microphysical models describing compaction of the Groningen and other producing (gas) reservoirs.

Highlights

• Time-lapse synchrotron micro-CT imaging reveals grain-scale deformation processes for compaction of Slochteren sandstone, from the Groningen gas field

• At small axial strains, digital volume correlation shows local strains by small displacements along grain boundaries, without intra-granular cracks

• After yielding, the sample deformed through pervasive grain failure and pore collapse

1 Introduction

Hydrocarbon production from sandstone reservoirs leads to pore pressure reduction, while the overburden stress remains constant. Consequently, production causes an increase in effective stress, which results in elastic and inelastic compaction of the reservoir with little or no lateral strain (Hol et al. 2018; Pijnenburg et al. 2018). At the surface, reservoir compaction may be expressed as surface subsidence, while gradients in deformation, normal, and shear stress across faults may lead to induced seismicity (Buijze et al. 2017; Evans et al. 2012; Grünthal 2014; Van Wees et al. 2014). Surface subsidence can reach values on the order of several tens of centimeters (Ketelaar 2009; Mallman and Zoback 2007), which means that compaction at the reservoir level typically corresponds to small vertical strains (< 1%). Since subsidence and induced seismicity arise from reservoir compaction, it is important to understand the deformation mechanisms operating in the reservoir during such deformation. Elastic deformation is well understood and can be quantified in terms of poro-elasticity (Wang 2000). However, inelastic deformation is less well understood, especially at reservoir-relevant strains (< 1%). To accurately model reservoir behavior during hydrocarbon production, it is therefore important to understand and quantify the grain-scale mechanisms controlling compaction: the 3D stress–strain behavior of the sandstone is dependent on the partitioning between elastic and inelastic deformation, and the precise mechanisms leading to inelastic deformation (Pijnenburg et al. 2018, 2019a; Pijnenburg and Spiers 2020). Indeed, these mechanisms can strongly influence lateral stress development during reservoir compaction, which in turn affects in-situ stress development on faults in the reservoir and their potential for seismic activation.

Prolonged production of gas from the Groningen gas field in The Netherlands has led to approximately 34 cm of subsidence, in response to a 25 MPa decrease in pore pressure, since production started in 1963. This subsidence has been accompanied by earthquakes in the central part of the field (Fig. 1), increasing in number and magnitude since the 1990s, with the largest magnitude event (M_w = 3.6) occurring near the village of Huizinge in 2012 (Cannon and Kole 2017; Dost and Kraaijpoel 2013; NAM 2015). The main reservoir rock in the Groningen gas field is the Slochteren sandstone, a formation within the Permian Upper Rotliegend Group. This is overlain by claystone of the Ten Boer Member of the Silverpit Formation followed by a 50 m thick anhydrite-carbonate and evaporite sequence which acts as cap rock (Buijze et al. 2017). The Slochteren sandstone reservoir is mostly between 100 and 200 m thick and located at about 3 km depth (Geluk 2007). The largest vertical strains in this reservoir are of the order of 0.2–0.3% (Cannon and Kole 2017; NAM 2016; Pijnenburg et al. 2019a), corresponding to a maximum subsidence of 34 cm at the center of the subsidence bowl (Fig. 1a), where the Slochteren sandstone is most porous (18%—22%) (Fig. 1b). Note that the overall compaction of the reservoir of at most 0.2–0.3%, since subsidence initiated has not decreased this porosity significantly.

Fig. 1

reproduced with permission from Pijnenburg et al. (2019a)

a The Groningen gas field, with the extent of the gas-bearing reservoir at a depth of approximately 3 km. Contour lines show the surface subsidence in centimeters as measured in 2013 and induced earthquakes are shown with their magnitudes and locations in the period of 1991–2017. b Map showing the mean porosity in the Slochteren reservoir sandstone, in combination with the main faults cutting the top of the reservoir. Figure

The typical macroscopic strain behavior of highly porous sandstones, such as those found in Groningen (porosity > 20%), includes roughly three stages, expressed in terms of mean effective stress and porosity reduction (Fig. 2). In the initial stage (stage 1) at low mean stress, the relationship between the porosity reduction and mean effective stress is generally non-linear concave-up. At higher confining pressure, this relationship transitions to a stage of near-linear behavior (stage 2). During the final stage (stage 3), non-linear behavior arises, leading to compaction and potentially localization of deformation, and ultimately to brittle failure. Although most triaxial experiments to date focus on understanding the failure and post-failure behavior seen at strains >  > 1% (stage 3) (Cała et al. 2016; Shi et al. 2021; Wong and Baud 2012; Yang 2018), in-situ strains relevant for depleting sandstone reservoirs typically fall within stage 2 behavior (< 1% strain). Generally, the linear behavior seen in stage 2 is interpreted as purely elastic deformation (Baud et al. 2006; Baud and Zhu 2000; Heap et al. 2015; Rutter and Glover 2012; Wong and Baud 2012). However, several experimental studies have shown that the near-linear relationship results from a combination of elastic and inelastic deformation (Bernabé et al. 1994; Hol et al. 2015, 2018; Pijnenburg et al. 2019a; Shalev et al. 2014), with the inelastic part constituting 20—70% of the total deformation, depending on the applied P–T-stress conditions and sample porosity (Bernabé et al. 1994; Hol et al. 2015, 2018; Pijnenburg et al. 2018).

Fig. 2

modified from Pijnenburg et al. (2019a)

Schematic representation of the typical mechanical behavior seen in triaxial compression tests on highly porous sandstones, expressed in terms of mean effective stress versus total porosity reduction. Figure

From “post-mortem” microstructural analysis of experimentally deformed samples, it has been inferred that the mechanisms operating at the small strains relevant for depleted reservoirs involve some combination of inter-granular (grain boundary) cracking, inter-granular slip, and intra-granular/trans-granular cracking (Bernabé et al. 1994; Brantut et al. 2013; Heap et al. 2015; Hol et al. 2018; Pijnenburg et al. 2019a; Shalev et al. 2014; Tengattini et al. 2014). However, given the micrometer-scale displacements between individual grains required to produce such small strains, it remains challenging to find direct evidence for the strain accumulating mechanisms. Furthermore, conventional analytical methods such as optical and electron microscopy only allow for comparison before and after experimentation, and only in 2D. With the development of high-resolution X-ray computed tomography scanning (µCT), it is possible to obtain a 3D representation of the rock sample (Cnudde and Boone 2013; Hanna and Ketcham 2017). This technique has been used to analyze undeformed versus deformed samples before and after the completion of an experiment (Baud et al. 2015; Liu and Zhao 2021). Most recently, developments in 4D, or time-lapse, µCT imaging have provided new opportunities to monitor grain-to-grain interactions during deformation, under realistic in-situ stress and temperature conditions (Andrade et al. 2011; Hurley et al. 2018; McBeck et al. 2018, 2019; Renard et al. 2004, 2009; Saadatfar et al. 2012; Takano et al. 2015; Van Stappen et al. 2018). In this manuscript, we illustrate the potential of using time-lapse (4D) µCT imaging to identify and quantify the deformation mechanisms that operate in producing sandstone gas reservoirs, such as Groningen, and in other sandstone reservoir systems, at the relevant in-situ conditions. Our results provide the first direct in-situ visualization of the precise deformation mechanisms operating in the center of the seismogenic Groningen reservoir at relevant P–T-stress conditions, or indeed in any depleting gas or other geo-energy reservoir.

2 Material, Experimental Methods, and Analytical Techniques

2.1 Material

For the single 4D experiment reported here, several cylindrical samples of Slochteren sandstone measuring 5 mm diameter and 10 mm length were prepared. Two samples were severely damaged when installed in the triaxial deformation apparatus. A third sample was successfully installed, with little damage, and the results are described in the present study. The sample material was taken from Slochteren sandstone core retrieved from the Zeerijp (ZRP)-3a well, located in the seismogenic center of the Groningen field (Fig. 1a), by the field operator (NAM, Nederlandse Aardolie Maatschappij) in 2015, i.e., after more than 50 years of production. The small cylindrical samples (5 mm diameter) were drilled sub-normal to bedding, from 25 mm-diameter plugs cored from the main NAM core. The main NAM core section used was taken at a true vertical depth of 3002 m and had a porosity of 22.9%. The small (5 mm diameter) cylindrical sample ultimately selected from this core section was similar to the intermediate porosity sandstone samples (porosity of 21.5%) tested in conventional triaxial compression by Pijnenburg et al. (2019a). More specifically, our sample replicates sample z24c of Pijnenburg et al. (2019a) in a qualitative manner, as described further below. For meaningful comparison, we aimed to follow the same hydrostatic and deviatoric loading path applied to z24c by Pijnenburg et al. (2019a).

Throughout the reservoir, the Slochteren sandstone is typically quartz-rich (72–90%), with varying amounts of feldspar (8–25%), clay (0.5–5.5%), and lithic fragments (3–10%) (Verberne et al. 2017; Waldmann et al. 2014). The sandstone has a predominant grain size distribution between 150 and 250 µm (Verberne et al. 2017). Some of the feldspar grains are chemically weathered to clay minerals, leading to micro-porosity. The clays are also present as clay films, with a thickness between 1 and 10 µm, coating quartz and feldspar grains, as well as many grain contacts (Gaupp et al. 1993; Waldmann and Gaupp, 2016). Figure 3A and B shows scanning electron microscopy images of a Slochteren sandstone sample taken from the same NAM core section as the plug used in this study. The composition is very similar to previous descriptions of the Groningen reservoir sandstone (Verberne et al. 2017; Waldmann et al. 2014), and closely resembles sample z24c tested by Pijnenburg et al. (2019a). The material is quartz-rich, with chemically weathered feldspar minerals and minor traces of a calcitic cement. From Electron-Dispersive X-ray spectroscopy (EDS; Fig. 3C), it is clear that microporous clay material fills in the pores and coats the quartz and feldspar grains.

Fig. 3

A scanning electron microscopy image of a Slochteren sandstone sample from the NAM core section sampled in this study. In the inset (B), minerals are defined through energy-dispersive X-ray spectroscopy (C). Clay minerals fill some pore space

2.2 Triaxial Experiment with 4D µCT Imaging

We performed a single triaxial compression experiment using the HADES deformation rig (Renard et al. 2016) at the micro-tomography beamline ID19 of the European Synchrotron Radiation Facility (ESRF) in Grenoble, France. This deformation rig is specifically designed to submit solid samples to stress and temperature conditions representative for depths up to 8–10 km (Renard et al. 2016). The rig fits into the ID19 micro-tomography beamline, allowing the sample inside to be fully imaged upon rotation between 0° and 180°.

In conducting the experiment, the 5 mm-diameter × 10 mm-long sandstone sample was first placed in a Viton jacket (50 mm long, 1.5 mm thick), sealed between the two end pistons of the HADES rig using O-rings and inserted into the HADES machine, which was subsequently emplaced in the ID19 micro-tomography beamline. The jacket isolates the sample from the confining medium (at pressure P_c = σ₃ = σ₂). Axial stress and confining pressure were controlled via two syringe pumps, using silicone oil as the pressure medium (0.1–100 MPa range; uncertainty of 0.05 MPa). Loading was under load control. Pressure-compensation balancing of the axial loading system means that the applied axial stress represents the differential stress (σ₁—σ₃). Two stigma pumps control the pore fluid pressure (P_p; distilled water) and inlet and outlet flow (300 ml capacity; 0.1–100 MPa range; uncertainty of 0.01 MPa) (Renard et al. 2016). Heating of the sample is achieved using two sheathed resistance elements at the top and bottom of the cell (maximum 250 °C; uncertainty 1 °C), minimizing any thermal gradient in the sample and confining fluid.

The Slochteren sandstone sample was brought to a temperature of 100 °C, followed by application of isostatic confining pressure (P_c = σ₁ = σ₂ = σ₃) to 10 MPa, then pore pressure P_f to 10 MPa, and then P_c = σ₁ = σ₂ = σ₃ to 50 MPa in 4 steps. During this procedure, seven 3D scans were acquired, one at each step in temperature, P_f or P_c. Subsequently, the axial stress was increased in 43 steps to 151 MPa (σ₁—σ₃ = 101 MPa) in increments of 1, 2, and eventually 5 MPa. After applying each axial stress increment, all stresses were held constant to take a µCT scan (2.5 min). In total, 50 µCT scans were acquired during the course of the experiment. With reference to Fig. 2, the loading conditions employed were such that the sample would deform in stage 2 and then approach a peak strength in stage 3 by compaction localization (Fig. 2), as described for sample z24c by Pijnenburg et al. (2019a). We emphasize that we could not reproduce the in-situ field boundary conditions or the evolution of these in the unique triaxial/CT setup used in the present study. Instead, we focused on stress conditions believed to be close to those in the reservoir in its later (current) depleted state (Pijnenburg et al. 2019a; Van Eijs 2015) and studied the microscale processes associated with (ongoing) Stage 2 mechanical behavior relevant for the Groningen field, as characterized by Pijnenburg et al. (2019a).

The HADES deformation rig, combined with two 2.8 mm-thick aluminum plates, acted as a filter limiting the full pink X-ray spectrum of the synchrotron beam, such that an average energy of 93 keV with peak energy of 80 keV and a bandwidth between 61 and 110 keV was transmitted through the rig and sample assembly. This setup and parameters have been shown to provide maximum image contrast between different phases typically present in geomaterials. Details on the setup and image acquisition methods are provided by Renard et al. (2016).

During the experiment, axial stress, confining pressure, and pore pressure signals were measured and logged during every scan. All µCT scans were obtained with a voxel size of 6.5 µm³ and a spatial resolution close to 6.5 µm. Sample shortening was determined in each scan using vertical cross-sections through the µCT data. The sample length at the different stages was compared to the initial sample length, determined when the pore pressure and confining pressure were, respectively, at 10 MPa and 50 MPa (i.e., at the start of the deviatoric loading). Macroscopic compressive axial strain, or shortening, is taken as positive. During the scanning time, it is possible that time-dependent deformation (i.e., creep) took place within the sample. However, significant creep would have resulted in motion artifacts, or blurring, in the µCT images, which were not observed. Indeed, it can be said that time-dependent deformation, if present, has to be limited to a fraction of the imaging resolution and thus smaller than at least half of the pixel size (or 3.25 µm). In the µCT images, the local grey value is related to the local attenuation of X-rays. Dense phases with high atomic numbers attenuate X-rays significantly and appear bright, while low-density materials with a low atomic number are darker. In this way, we can distinguish between water-filled pores and minerals and we can follow the evolution of porosity. To do so, the same spatial volume within the sample was taken in each scan (i.e., the reference frame and voxels do not deform, but the sample within this frame does). The pore space of the rock was segmented with a grey value threshold < 16,000. This relatively low threshold results in minimum estimates of porosity. Using a higher threshold resulted in erroneous segmentation of grain edges and microporous feldspar grains as pore space. For this reason, the low threshold was chosen. Choosing this lower threshold provides a minimum estimate of the actual porosity in the sample. We must note that all µCT scans were obtained and reconstructed with identical parameters, so that using a fixed threshold (< 16,000) is meaningful, even though deformation takes place during the experiment.

Figure 4 illustrates a vertical cross-section through typical µCT data (obtained at a differential stress of 5 MPa, i.e., after the first increment of deviatoric loading). On the basis of the grey values and grain morphology, we can differentiate the following constituents in the sandstone sample, from dark (low grey value) to bright (high grey value): brine-filled pores, quartz and feldspar, zircon and barite (bright spots), and finally the steel pistons above and below the sample. Because of their similarity in grey value, it is very difficult to differentiate between quartz and feldspar. Some of the feldspars are altered and therefore microporous. These can be differentiated from the quartz grains on this basis, as illustrated in Fig. 4B. In terms of microstructure, both the sample grain size and pore size are spatially heterogeneous, as seen especially in the bottom half of the sample (inset B versus inset C in Fig. 4).

Fig. 4

Representative axial cross-section through the µCT scan of the sandstone sample (left panel). Inset A shows brine-filled pores (related to a horizontal fracture) and bright barite minerals; inset B shows the differentiation between quartz and feldspar minerals based on the feldspar micro-porosity; inset C illustrates a section of the rock with a grain size typically lower than 150 µm, as opposed to the grain size observed in inset B

It can furthermore be observed from Fig. 4 that the sandstone sample was partly damaged: material was lost from the top part of the sample and a horizontal fracture can be observed just above its center (inset A in Fig. 4). The sample material was of low cohesive/tensile strength, which resulted in this damage upon emplacement in the HADES triaxial machine. Other samples were rejected during placement into the machine as damage was more pervasive.

2.3 Image Analysis: Digital Volume Correlation

A time-lapse record of the relative displacement fields of sample sub-volumes, i.e., volumes of 10-by-10-by-10 voxels, in the µCT images was obtained using the digital volume correlation (DVC) analysis code TomoWarp2 (Tudisco et al. 2017). This code has previously been used in geo-mechanical studies to understand strain localization in sand packs, shale, and sandstones subjected to triaxial testing (Andò et al. 2017; McBeck et al. 2018, 2019; Renard et al. 2019; Tudisco et al. 2015). The TomoWarp2 code is based on a local DVC approach, in which the volume of a µCT scan is subdivided in small sub-volumes (or nodes), which are then independently registered in a subsequent scan (Buljac et al. 2018; Tudisco et al. 2017). The DVC analysis by TomoWarp2 results in a list of spatial coordinates for each of these sub-volumes, their displacements and rotation components to reach their new position in the subsequent scan, and the correlation coefficient, which determines the quality of the fit (Tudisco et al. 2017). In our analyses, a node spacing (i.e., spacing of the different sub-volumes) of 20 voxels (i.e., 130 µm) was used in combination with a correlation window (i.e., volume taken into account per sub-volume) of 10 voxels (i.e., 65 µm). The search window (i.e., an estimate of the expected displacement of the sub-volumes) was altered according to the investigated deformation step. All voxels representing the rock sample were taken into account in the DVC analysis. Based on the size of the largest particles in the sample, combined with the spatial resolution of the images, selecting a higher spatial resolution discretization would have resulted in sub-volumes completely inscribing a single grain, without encompassing grain–grain or grain–pore contacts. Such sub-volumes are prone to introducing errors during the DVC analysis, due to a lack in strong differences (i.e., phase transitions) in grey value, making recognition in subsequent scans by the TomoWarp2 code more difficult. To understand the detection limit of the method for small local displacements, the TomoWarp2 code was run with one micro-CT volume as input. Any detected displacement between the identical volume should thus be considered as measurement errors. Figure 5 gives the measured displacements in the X-, Y-, and Z-direction obtained from this error measurement. Maximal displacements of 0.05 voxels (= 0.325 µm) were obtained. Therefore, in this work, we will consider displacements down to 0.1 voxel (0.65 µm) to be detectable in the DVC measurements between different micro-CT volumes.

Fig. 5

Histogram of the measured displacements for the DVC analysis on the same micro-CT volume. Maximal displacements of 0.05 voxels (= 0.325 µm) were observed in each Z-, Y- and X-direction

The output of the TomoWarp2 code can further be processed following the procedures described by McBeck et al. (2018). In these procedures, the displacement fields are used to calculate local strains within the sample. DVC allows sub-volumes with local dilation or contraction to be identified, as well as quantification of the degree of rotation of the sub-volume, based on the 3D displacement fields obtained from the TomoWarp2 analysis (McBeck et al. 2018; Tudisco et al. 2017). In the post-processing, we focus on the largest values (top 3%) and the mean values of local dilation, contraction, and rotation, because we investigate deformation between scans separated by both small axial strains (< < 1%) and larger values (> 1%). To determine the magnitude of local dilation or contraction, the divergence (div) of the 3D displacement field is calculated, yielding the first invariant of the strain tensor. Contrary to the macroscopic strain, which is considered positive for macroscopic compaction, a negative divergence corresponds to local contraction, while a positive divergence indicates local dilation (McBeck et al. 2018). Local contraction and dilation correspond to volumetric strain changes. To determine the magnitude of local rotation, the curl of the 3D displacement field is calculated. The curl reported here is a measure of angular rotation increment, hence angular velocity, within the horizontal XY-plane of the μCT images. A positive curl represents right-lateral rotation in the XY-plane, while a negative curl represents left-lateral rotation, in the reference frame of the sample (McBeck et al. 2018). Figure 6 illustrates this concept with an artificial displacement field, given by the gradient of the function $f\left(x,y\right)=2x {e}^{-{x}^{2}-{y}^{2}}$. This illustration is meant to visually show the concept of divergence and curl, but does not necessarily represent a displacement field one can expect in a triaxial compression experiment.

Fig. 6

Divergence and curl calculated for an artificial displacement field. Negative and positive divergence are caused by contraction and dilation in the displacement field, respectively. Left-lateral and right-lateral rotation result in a negative and positive curl, respectively

The divergence and curl measurements are subsequently divided by the change in macroscopic axial strain (or macroscopic contraction), Δε^M_zz measured between the two μCT images used in each DVC calculation. In this way, the local measures of divergence and curl are normalized with respect to the axial displacement observed in these two μCT images (McBeck et al. 2018). This step reduces the influence of changes in the Δε^M_zz between scan pairs on the reported magnitude of the local strain components.

The post-processing of the TomoWarp2 output allows us to identify volumes of interest where strain is locally accommodated. These volumes of interest can then be visually inspected to help determine the grain-scale deformation mechanisms at play. For clarity, it must be noted that strains reported in the DVC analysis are incremental and not cumulative. They represent the strain accumulated between two scan acquisitions, and not the total cumulative strain experienced by the rock. As the considered macroscopic and local strains are small, we focus on the top 3% of maximum local strain to determine strain accumulation. This yields manageable and displayable amounts of data while showing the trends qualitatively visible in larger percentiles of the data.

3 Results and Discussion

3.1 Macroscopic Stress–Strain Data and Time-Lapse µCT Imaging

The stress–strain evolution of the sample during deviatoric loading is shown in Fig. 7, in terms of differential stress (σ₁—σ₃) versus axial shortening. The sample was able to support a differential stress of 92 MPa, attained at 4.37% axial strain, before sharp yielding occurred, leading to > 10% additional axial strain as the differential stress was increased to the maximum applied value of 101 MPa. In this interval, two µCT scans were acquired, at differential stresses of 97 MPa and 101 MPa and axial strains of 14.54% and 17.71%, respectively. Then, the sample was axially unloaded, depressurized, cooled, and recovered from the HADES rig.

Fig. 7

Stress–strain evolution of the sample, after the initial hydrostatic loading of 50 MPa. Prior to pervasive yielding/failure beyond 92 MPa differential stress (left panel), the experiment can be subdivided into three different phases, separated by sample-scale deformation events of significant magnitude, as illustrated by the more detailed plots of the individual phases, A, B and C. Each black data point indicates the stress–strain coordinates corresponding to acquisition of a 3D scan; and the 11 red data points correspond to scans used for digital volume correlation (DVC) analyses. Ten DVC calculations were calculated from pairs of 3D scans corresponding to two successive red points (see Table 1). 0% axial shortening is considered after hydrostatic loading of 50 MPa of the sample

In the present experiment, stage 1 behavior, i.e., microcrack closure (Fig. 2), occurred during hydrostatic loading and was influenced by the damage induced during the experimental setup. The hydrostatic loading is not shown in Fig. 7, which only shows the deviatoric stress loading path of the experiment. Prior to sharp yield at 92 MPa differential stress, the deviatoric stress loading path is characterized by three active loading plus axial deformation phases (A, B, and C) featuring shortening strains up to 0.1% (see Fig. 7). These are bounded by the scans acquired at fixed differential stresses of 69 MPa, 86 MPa, and 92 MPa. The larger axial deformation events that occurred during scanning at 69 and 86 MPa took place rather abruptly, during scans 32 and 43, respectively, leading to motion artifacts within these scans. Therefore, additional scans (numbers 33 and 44, respectively) were retaken immediately afterward at the same stress conditions.

At the start of deviatoric loading, stage 2 behavior is observed (refer Figs. 2 and 5, cf. Pijnenburg et al. (2019b)). The sample shows a quasi-linear stress–strain path in phase A, up to a cumulative macroscopic axial strain (ε^M_zz) of 0.05%, when non-linearity begins. At 69 MPa, larger strain (ε^M_zz > 1%) is accumulated suddenly. It must be noted that 0% axial shortening is considered after hydrostatic loading of 50 MPa of the sample. The sudden increase in ε^M_zz separating phases B and C, without concurrent increases in differential stress, indicates that the sample experienced significant internal deformation events producing sample-scale yield. This yielding occurred at similar stress conditions as in the reference experiment on sample z24c described by Pijnenburg et al. (2019a) (~ 50—55 MPa differential stress—cf. Figs. 2 and 7), and in various other sandstones of similar porosity (see Fig. 5b in Wong and Baud 2012). This provides confidence that the deformation behavior shown by the presently studied sample is representative for a wider range of sandstones of similar porosity, tested at similar stress conditions. We consider phases B and C to fall within the so-called stage 3 behavior (Fig. 2).

We illustrate the time evolution of a µCT cross-section taken along the vertical axis and through the middle of the deforming sample in Fig. 8 and complementary in video S1, presented in the Supplementary Information. This evolution shows that the sub-horizontal fracture present in the undeformed sample, which developed while setting up the experiment (e.g., Fig. 4), closes in the axial strain event between phases A and B (Fig. 7), and that subsequent compaction of the sample occurs in three consecutive events, consistent with the mechanical data shown in Fig. 7 (phases B, C and beyond C).

Fig. 8

The evolution of a single vertical cross-section through the µCT in function of the applied differential stress. The pre-existing fracture and the upper piston are indicated in the cross-section at 25 MPa differential stress. The fracture is seen to close throughout the experiment, which ends with an increase in total axial strain up to 17.71% at a differential stress of 101 MPa

Consistent with Fig. 8, the porosity profile, determined in horizontal slices at various heights above the sample base (Fig. 9), shows the presence of a highly porous, horizontal fracture through the middle of the sample, in the early stages of differential loading, i.e., in scans 8 to 31 (i.e., phase A). Comparing scans 8 and 31 (the first and last scan of phase A, respectively), we observe a slight porosity decrease in the middle of the sample, where the fracture is located. During the axial strain event between phase A and phase B, the fracture closed (cf. Figs. 7, 8 and 9). Furthermore, it is clear from Fig. 9 that the sample consisted of a more porous (probably damaged) portion above the initial fracture, and a less porous portion below the fracture. With ongoing differential loading, the porosity in the upper zone decreased more than the porosity in the lower zone up to scan 47, up to an effective differential stress of 92 MPa, that is the end of phase C. Only beyond phase C, i.e., in scans 48 and 49 at stresses > 92 MPa, do we observe a distinct decrease in porosity (i.e., compaction) in the lower zone of the sample, coinciding with an almost total loss of imaged porosity in the upper zone. Figure 9 shows that the initial porosity and the degree of initial damage in the rock played an important role controlling the subsequent shortening and change in porosity (i.e., compaction) upon differential loading of the sample, i.e., in controlling macroscopic stress–strain behavior.

Fig. 9

Evolution of the relative amount of pore space in a horizontal rock slice (thin disc) of fixed volume translated along the sample axis, as a function of strain or scan number in the experiment. Each slice has a thickness equal to the resolution of the scan (i.e., 6.5 µm). The relative amount of pore space is shown for each horizontal slice through the rock volume. The top and bottom of the sample is indicated as well as the location of the fracture within the scans. The total axial strain experienced by the sample is given below each plot (in black: at start of scanning interval; in red: at end of scanning interval)

3.2 Strain-Accommodating Mechanisms Observed Through µCT Imaging

A summary of the stress–strain values associated with the scans used in the DVC analysis is shown in Table 1 and in Fig. 7. We performed the DVC analysis with consecutive scan pairs (i.e., 08–11, 11–17, and so on until 48–49, see Table 1). These pairs show the successive deformation within phase A, constitute the beginning and end of each phase (A, B and C) and extend to the end of the triaxial test, as described in the previous section. We also performed DVC using scans 08 and 31, which are the first and last scan of phase A, respectively. The TomoWarp2 code was run across the entire µCT volume representing the rock sample. Post-processing of the DVC results enabled differentiation between zones present within the sample, including the fractured middle zone, the undamaged lower zone, and the damaged upper zone. The results of the DVC analysis are considered in the sections that follow.

Table 1 Selected scans, and the corresponding stress–strain levels used for DVC analysis. Each set of two consecutive scans was taken as a scan-pair for input to the analysis. In addition, one extra DVC analysis was carried out between scans 08 and 31

3.2.1 Micro-displacements During Small Strain Deformation

Reservoir depletion is typically associated with small-strain deformation (< < 1%). In the present experiment, three small-strain deformation phases (A, B, and C), each with < 1% macroscopic strain, occurred (Fig. 7). These phases provide the opportunity to quantify the microscale displacements and identify the mechanisms that produce this macroscopic behavior. However, due to the sample damage that developed prior to the experiment, a significant portion of the deformation observed during loading is controlled by the induced horizontal fracture and its associated closure (Fig. 8). Furthermore, the sample volume located above the fracture is likely not representative of the intact sample material, also due to pre-loading damage. Therefore, in this section, we focus on the rock volume below the horizontal fracture, i.e., the bottom part of the sample. By selecting this part of the sample for post-processing of the DVC analyses, we are able to isolate the deformation mechanisms operating in intact sample material at small macroscopic strains, as relevant for deformation of the sandstone under in-situ reservoir conditions. The final result of the DVC analysis comparing scans 08 and 31 (first and last scan of phase A, respectively) is shown in Fig. 10.

Fig. 10

Result of the post-processing of the DVC analysis comparing scans 08 and 31 (differential stresses of 10 MPa and 69 MPa, respectively) for the undamaged bottom portion of the sample tested. Only the locations of the top 3% values of local dilation (in blue), contraction (black) and right-lateral and left-lateral rotations of the sub-volumes (in red and purple, respectively) are shown for clarity

Figure 10 shows the location of the sub-volumes used in the DVC analysis in which the top 3% values for dilation, contraction, and rotations were measured when comparing scan 08 (start of phase A) with scan 31 (end of phase A). These locations represent where in the sample most local deformation took place. It is interpreted that there is no clear localization of deformation, or link to sedimentary structure to be determined from these images. Deformation occurs over the entire section of the sample below the fracture. For the DVC analyses performed at sub-stages of phase A (i.e., DVC analysis between scan pairs 07–11, 11–17, 17–24 and 24–31), the magnitude of the top 3% most significant deformation increases slightly with axial strain (cf. Fig. 11A). When taking into account all deformation in the lower subsection of the sample, average values for the local dilatational divergence, contractional divergence, and positive and negative curl can be calculated. These average values are smaller for the sub-stages of phase A than for the complete phase A deformation (cf. Fig. 11B compared to Fig. 11D). Furthermore, the magnitudes of positive and negative curl (Fig. 11A–C), are roughly equal, but opposite in sign, suggesting that these motions cancel each other out on the macroscopic level. Overall, more local contraction is observed than dilation.

Fig. 11

Magnitudes of div/Δε^M_zz and curl/Δε^M_zz for a detailed breakdown of phase A (DVC comparing scan pairs 08–11, 1–17, 17–24, and 24–31) measured in the sub-volumes used for DVC where the top 3% most significant local deformation was measured compared to the average deformation measured in all sub-volumes in figures A and B, respectively; Figures C and D show the magnitudes for the 3% most significant deformation compared to the average values for the magnitudes of div/Δε^M_zz and curl/Δε^M_zz for phases A, B, and C. These values indicate very small local deformation during small-strain deformation of the sample (< < 1%)

Comparing the magnitudes of local deformation for the three small-strain deformation phases, i.e., A (scan-pair 08–31), B (scan-pair 33–39), and C (scan-pair 44–47), reveals that the incremental local strain values remain fairly constant throughout the experiment, with slightly more compaction observed during phase A (Fig. 11C and D). On average, there is more right-lateral rotation than left-lateral rotation during phases A and B. Visual inspection of the locations of the motions indicates that deformation is not localized or linked to sedimentary structures (the DVC results for phases B and C are shown in the supplementary information, Figure S2).

For the locations where the 3% highest local deformation is measured, the absolute values of div/Δε^M_zz and curl/Δε^M_zz measured across all of phase A (shown in Fig. 11C) are four-to-five times lower than the sub-stages investigated in phase A (shown in Fig. 11A). The mean values, however, are more similar (Fig. 11B and D). There is no significant difference between the behavior observed in scans prior to and after the yield point identified from the macroscopic behavior (cf. phase A vs. phases B and C).

Based on the x-, y-, z- coordinates of the local deformations represented in Fig. 10, the µCT images can be evaluated to identify the strain-accommodating grain-scale mechanisms. Consider, for example, Fig. 12 for scan-pair 08–31. In this figure, a part of the sample approximately 600 µm below the horizontal fracture is shown. Two of the most significant local contraction events in the intact part of the sample, determined with the DVC analysis, are located here (x-, y-, z- coordinates 3.77, 4.55, 5.65 and 4.42, 4.55, 5.65, respectively, in Fig. 10). These locations are indicated with circles in Fig. 12. It is difficult to observe changes directly in static, 2D representations of the microstructure of the sample (Fig. 12). However, we added the x-, y-, and z- displacements determined by DVC. These displacements are larger than the determined error of 0.05 voxels (Fig. 4) and can thus be considered as true displacements. Figure 12 shows that all grains experience a downward component of displacement. This occurs, while also macroscopic shortening is measured, through the displacement of the upper piston, while the lower piston remains stationary. The local volumetric contraction measured at the site of Fig. 12 is inferred to be governed by small inter-granular displacements, though it is difficult to see whether they are normal and/or shear displacements. In the direct neighborhood of the locations identified having the most volumetric contraction, no fractures or intra-granular cracks can be observed. The inter-granular displacements are inferred to result in a minor local decrease in pore volume, due to closure of slightly open grain boundaries or compaction of inter-granular clay layers (Pijnenburg et al. 2019b; Verberne et al. 2021). However, this closure cannot be directly observed at the current spatial resolution of 6.5 µm, as the displacements associated with these movements would be on the order of a few hundred nanometers or, at most, a few micrometers, as also indicated by the x-, y-, z- displacements in Fig. 12. Visual inspection of the µCT scans of the local contraction and dilation locations identified in DVC (locations based on the 3% highest values of deformation, Fig. 10) provides no clear evidence that contraction or dilation was linked to local grain failure, as no new cracks can be observed (µCT images in Fig. 12). Although fractures might facilitate small grain displacements during phase A, they are inferred not to be the main mechanism of deformation in the intact part of the sample, below the sub-horizontal fracture.

Fig. 12

Vertical µCT slices from scans 08 and 31. The circles show the locations of maximum local contraction, as measured between scan 08 and 31. As can be seen, no new intra-granular cracks can be observed in Scan 31 over Scan 08, meaning that grain failure therefore likely did not contribute significantly to the local deformation. Below the µCT images, the displacement in the z-, y-, and x-direction is shown, necessary to go from scan 08 to scan 31. All displacements are expressed in voxel size (= 6.5 µm). The displacement is limited in all directions; however, our interpretation is that local contraction must therefore have been accommodated by inter-granular normal and/or shear displacements

3.2.2 Fracture Closure and Deformation at Intermediate Strains

At the fixed stresses maintained between the small-strain deformation phases A and B, and B and C, two “interphase” deformation events occurred that led to a substantial amount of sample shortening strain, namely 3% and 1.5% contraction, respectively. This macroscopic contraction is associated with closure of the horizontal fracture within the sample (Figs. 13, 14). Figure 13 shows the DVC results for these interphase events. Similar results for the local dilation and left-lateral rotation are included in the supplementary information (Figure S3, supplementary information). High local contractions, significant grain rotations, and grain failure occur in the material directly surrounding the fracture (top row in Fig. 14 for scan-pair 39–44). For both interphase events, localized deformation was also observed in the intact, bottom portion of the sample, i.e., below the fracture (bottom row in Fig. 14). This deformation was more localized than the deformation observed during the small-strain deformation phases (Fig. 10).

Fig. 13

Result of the post-processing of the DVC analyses comparing A scans 31 and 33, and B scans 39 and 44. At the full sample-scale, strain is localized around the early horizontal fracture, while the deformation in the lower part of the sample is less localized (see also Fig. 9). Only the locations of the top 3% of the local contraction (black) and right-lateral rotation (red) values are shown for clarity

Fig. 14

Vertical µCT slices from scans 39 and 44, through the rock volume near the pre-existing horizontal fracture (top row) and in the bottom portion of the sample (bottom row). The slices illustrate local deformation at the intermediate strains between phase B and C of loading. In the direct vicinity of the early fracture, strain is accommodated through grain failure and displacement leading to local pore collapse and pore filling with grain fragments, while in the intact, lower sample area, only grain displacement can be observed resulting in small porosity changes

Furthermore, local contraction in the intact, bottom part of the sample is also larger in magnitude compared to the small-strain behavior observed during phases A, B and C. However, deformation in the bottom part of the sample is of smaller magnitude than the deformation around the fracture. High local contractions are measured close to the horizontal fracture between phases B and C (see top row in Fig. 14 for scan-pair 39–44). The magnitude of the vertical deformation in the rock volume surrounding the fracture (i.e., div/Δε^M_zz) is up to ten times higher than for the deformations in the intact, bottom portion of the sample (– 43.39 vs. – 3.6, respectively, Fig. 14). Local contractions near the fracture, leading to increased local stresses exerted on individual grains, resulted in grain breakage (top row, Fig. 14). Local porosity reduction occurs as fragments of broken grains fill in space between grains, and/or due to small grain rotation and translation without breakage. In this regard, we believe also the interplay between competent (e.g., quartz) and weak (e.g., feldspar) grains can be important: the more competent grains can displace through small grain rotations and translations, exerting local stress on surrounding weak minerals which then typically deform in a brittle manner (Hangx et al. 2010).

In the lower half of the sample, local contraction is predominantly accommodated by grain rearrangement without significant grain breakage or cataclasis. Although the spatial resolution is limited to 6.5 µm, we infer that the rearrangement of the stronger quartz grains involves inter-granular slip and closure of open grain boundaries, thereby stressing the weaker feldspar minerals (Verberne et al. 2021). The feldspar minerals are chemically weathered and therefore contain micro-porosity (Fig. 14, bottom right), which could be closed upon loading. Again, there is no clear evidence in the µCT scans for grain breakage being an important strain-accommodating mechanism in the intact, bottom part of the sample, implying that Stage 3 behavior represented in Fig. 2 was not achieved in this part of the sample.

3.2.3 Sample Deformation at High Stresses and Strains

Increasing the differential stress from 92 MPa (scan 47) to 97 MPa (scan 48) resulted in pervasive grain failure (Fig. 15), leading to an increase in axial strain of 10%. Increasing the differential stress to 101 MPa (scan 49) resulted in an additional 3% axial strain.

Fig. 15

A Vertical cross-sections through µCT scans 47 (92 MPa), 48 (97 MPa) and 49 (101 MPa) taken at the center of the sample; B results of the DVC analyses (local dilation in blue; contraction in black; right-lateral rotations in red; left-lateral rotation in purple); C insets a–a′, b–b′, and c–c′ showing the local evolution of the microstructure, displaying pore collapse, brittle grain failure of the quartz grains, and feldspar crushing throughout pervasive failure of the sample

Vertical cross-sections made through the center of the sample after the large compaction events at 97 and 101 MPa (Fig. 15A, right two images) illustrate ongoing grain crushing and related local contractions and rotations, notably in the top part of the sample (scans 47 and 48 in Fig. 15B). Strain in the top part of the sample is accommodated through pore collapse, brittle grain failure of the quartz grains, and extensive feldspar crushing (Fig. 15C, images a-a′ and b-b′). Larger quartz grains appear to be less prone to failure (Fig. 15C, image b’). Significant pore volume reduction during grain failure can be observed in the upper portion of the sample (Fig. 9, scans 47 and 48), while porosity appears to be fairly well preserved in the lower part of the sample (Fig. 15C, images a′ and b′ vs. c′).

Although large-scale sample failure is initially concentrated in the top portion of the sample and led to a significant amount of compaction (scans 47 and 48 in Fig. 15A), reaching a total macroscopic strain of 15.4% by the end of phase C, subsequent local strain development concentrates in the lower half of the sample (scans 48 and 49 in Fig. 15A). At high differential stress, volumetric contraction in the lower portion of the sample localizes along a sub-planar surface (see Fig. 15B, as indicated by the DVC analysis between scans 48 and 49), with mainly right-lateral rotation above and left-lateral rotation below, suggesting the development of a right-lateral shear fracture. The orientation and location of the shear plane shows a marked correlation to the distribution of local contraction observed between scans 39 and 44 (78 and 86 MPa differential stress; Fig. 13B). Detailed investigation suggests that the apparent shear zone is located in a zone with a finer average grain size (Fig. 15A, scan 48 vs. 49, at 97 and 101 MPa differential stress). The formation of the shear plane appears to again be related to local pore collapse, brittle failure of quartz grains, and feldspar crushing (see Fig. 15C: insets c to c′). Its development is consistent with the previous work in that, at such large strains, conventional triaxial experiments on sandstones typically show strain localization in shear and/or compaction bands formed by micro-fractures nucleating at the scale of several grains followed by fracture coalescence (McBeck et al. 2018; Renard et al. 2019).

4 Conclusions

This study demonstrates the application of in-situ 4D µCT imaging (spatial resolution of 6.5 µm) to help understand the grain-scale deformation processes responsible for compaction of Slochteren sandstone, obtained from the seismogenic Groningen gas field, when subjected to axial compression in a small-scale triaxial apparatus (P_c = 50 MPa, P_p = 10 MPa, T = 100 °C). The experiment allowed us to observe and quantify strain-accommodating mechanisms at in-situ stress conditions relevant for producing sandstone reservoirs. Using Digital Volume Correlation (DVC) analyses, results pinpoint the areas undergoing the most significant deformation during compression of the material. In this way, substantial information was gained on deformation and the associated grain-scale deformation mechanisms operating in the sample, though the resolution of the scans limited what was possible. Small displacements along grain boundaries could only be observed via DVC analysis, and clay layers present between grains and potentially accommodating displacement were not distinguishable in the µCT images. Further developments are necessary to increase the level of detail in the µCT images. Nonetheless, we could rule out grain failure as inelastic deformation mechanism operating under the conditions of the Groningen field in this experiment. At differential stresses in the range of 97–101 MPa, and axial strains in the range 14.5–17.7%, which are well in excess of those expected in clastic gas fields, we also observed the onset of a shear zone associated with failure of the sample. In summary, our conclusions are the following:

• 1. At small macroscopic axial strains (< < 1%), grain-scale strains were not localized, but rather homogeneously spread throughout the sample. At intermediate macroscopic axial strains (in the order of 1%), the distribution of local strains throughout the sample was more heterogeneous. However, it was only after sample failure that we could clearly recognize strain localization, in this case in the form of a shear zone.

• 2. For small macroscopic axial strains (< < 1%), relevant for gas depletion in sandstone reservoirs, local strains are governed by small grain displacements occurring via inter-granular slip, closure of open grain boundaries or compaction of inter-granular clay layers. This last result is an interpretation, because it was below the spatial resolution of the imaging method. At these strains and spatial resolutions, we did not observe intra-granular cracking.

• 3. Although the 4D imaging technique employed in this study shows potential for studying grain-scale deformation processes at realistic in-situ conditions (i.e., small-strain conditions), at present, it is limited by the spatial resolution. For imaging of grain displacements relevant for the very small strains occurring during reservoir depletion, the spatial resolution should be higher, i.e., in the submicron range. Potentially, by focusing on smaller volumes of the sample during imaging, it may be possible to achieve the desired resolution, though this has not yet been verified and requires adaptations to the current experimental setups.

• 4. At intermediate strains (1–3%), local strain development was related to the closure of a pre-existing fracture, which developed during sample assembly by intra-granular cracking and subsequent grain displacement through inter-granular slip. In intact parts of the sample, only grain displacement through inter-granular slip was observed without significant intra-granular cracking.

• 5. At high strains (> 4%), local strains were accommodated by extensive pore collapse and the development of a shear zone. The grain-scale strain-accommodating mechanisms involve widespread intra-granular cracking of the quartz grains, coupled with crushing of weaker feldspar grains.

• 6. Despite limitations of the present work, our results do confirm the findings reported by Pijnenburg et al. (2019b), who performed incremental triaxial tests on split cylinders of Slochteren sandstone under simulated reservoir depletion conditions. From SEM observations performed at different stages of loading through Stage 2 deformation behavior relevant for the Groningen field, these authors also inferred that macroscopic strain and compaction is accommodated by inter-granular clay film compression and shear. Confirming these findings strengthens the basis for development of constitutive models for reservoir compaction that are based on the physical processes actually operating.