Introduction

Compaction bands are defined as tabular zones of purely localized compressive deformation that are developed perpendicular to the direction of maximum compressive principal stress (Issen and Rudnicki 2001; Eichhubl et al. 2010; Wong and Baud 2012). Compaction bands play a critical role in controlling fluid flow. Compaction bands can have a contrasting impact on fluid flow depending on their geometry (Baud et al. 2012). The majority of studies emphasize the local decrease in porosity and consider compaction bands as flow barriers for fluid flow in reservoirs (Vajdova et al. 2004; Sternlof et al. 2006; Deng et al. 2015). Detrimental effects of the compaction bands have also been reported in the wellbore region during injection or extraction of fluid for storage or hydrocarbon production (Haimson 2001; Haimson and Kovacich 2003; Rahmati et al. 2014). Some authors also attribute this phenomenon to the mechanism of sand production in hydrocarbon reservoirs (Li et al. 2006; Detournay 2009). It is therefore important to understand the role of heterogeneity for the nucleation of compaction bands which is still controversially discussed (Fortin et al. 2006; Cilona et al. 2014; Baud et al. 2015). In this contribution, we present a time-lapse micro-CT experiment of compaction of a limestone (Mt Gambier limestone) of high porosity of ~ 50% to investigate this problem using digital imaging technology.

Time-lapse Micro-CT imaging of triaxial experiments provides new insights into the 4D-evolution of compaction bands in porous media including their effect on permeability. A recent study by Huang et al. (2019) used this technique to document the evolution of compaction bands in porous limestone (Leitha limestone) featuring an initial porosity of 22.5% and a macropore size of 165 μm. The study clearly showed the effect of the strong heterogeneity of local porosity near the top of the sample on the formation of distributed and localized compaction. A first region of collapsing highest porosity clusters in the top part of the sample was observed followed by the second stage of collapse of clusters of macropores in the lower part. The interesting role of local heterogeneity on the dynamics of pore collapse also has strong implications on the permeability evolution of the sample.

Compaction bands were first reported through field observation in Jurassic Aztec sandstone of south-eastern Nevada by Hill (1993) and later studied by others (Sternlof et al. 2005; Eichhubl et al. 2010; Sun et al. 2011). Compaction bands have been also observed in naturally deformed carbonates (Tondi et al. 2006; Rustichelli et al. 2012).

Understanding the control of compaction bands on fluid flow and possible sand production is crucial for oil and gas production, and therefore, numerous triaxial experiments have been conducted in the past supported by thin section analyses post deformation. These laboratory triaxial experiments, mimicking geological stress conditions, have been performed using different rocks, e.g., sandstone (Baud et al. 2004; Fortin et al. 2005; Townend et al. 2008; Charalampidou et al. 2011), carbonates (Ji et al. 2012; Baud et al. 2017), soil (Arroyo et al. 2005), clay (Hashimoto et al. 2019), and mudstone (Oka et al. 2011). Interestingly, compaction bands have been observed in applications outside of geoscience areas such as foam (Bastawros et al. 2000), honeycomb (Papka and Kyriakides 1998), snow (Barraclough et al. 2016), and puffed rice (Guillard et al. 2015).

The controlled laboratory experiments have shown that the porosity and permeability of the compaction bands can be significantly reduced by internal pore structure changes inside the band with potential pore collapse (Holcomb and Olsson 2003; Baxevanis et al. 2006; Lenoir et al. 2010; Chen et al. 2017). The accompanied permeability reduction can be empirically linked to axial strain recorded during the formation of compaction bands in laboratory experiments (Vajdova et al. 2004; Fortin et al. 2005). Tembe et al. (2008) noted that the thickness and length of compaction bands obey a quadratic scaling relation. In a similar study, Townend (2008) showed that the compaction bands nucleate at the edge and propagate across the specimen at 0.08 mm/s by collecting acoustic emission responses for Diemelstaedt sandstone samples. In other experiments, localized compaction bands were triggered in a pre-determined location of the specimen by carving a geometric imperfection (e.g., V-shaped circumferential notch) around the middle of the sample (Vajdova and Wong 2003; Tembe et al. 2006; Charalampidou et al. 2011; Charalampidou et al. 2014). Baud et al. (2004) and Fortin et al. (2006) identified three different patterns of bands in triaxial testing including discrete compaction bands, diffuse compaction bands, and high-angle shear bands. More recently, McBeck et al. (2018, 2020) reported that the proximity to macroscopic failure can be predicted using local strain information from dynamic in situ X-ray tomography. The proposed technique is to combine triaxial compression experiments on different rock types with the machine learning analyses. The results suggested that similar strain accumulation processes may precede macroscopic failure. However, the technique did not allow direct measurement of changes in the local stress fields of the samples.

Pioneering works with the aim of measuring the local stress variations were conducted to investigate the fracturing and crushing of grains and their effect on the macroscopic behavior of granular media (Hurley et al. 2018). The method relied on the use of analog granular materials (e.g., sapphire) with known mechanical properties so that the local stress field can be inverted from the X-ray CT images. Unfortunately, this technique does not lend itself to our experiments with highly porous limestones. The local material properties are unknown as it is a heterogeneous material that has complex microstructures, microfossils, and morphology. Our method is therefore similar to existing techniques for assessing meso-scale localized strain fields down to the strain of individual grains. This is accessible using X-ray tomography and DIC analysis at the grain scale (Hall et al. 2012; Takano et al. 2015).

In the past decades, the majority of studies of compaction bands were conducted on sandstone. Cheung et al. (2012) concluded that a homogeneous grain size distribution is the main microstructural attribute promoting the development of compaction bands in sandstone. However, compaction localization has not been investigated to the same level in a more complex microstructure such as a carbonate (e.g., mineralogy, deposition, and diagenesis) (Rashid et al. 2017). New findings may be expected for a highly porous and extremely heterogeneous carbonate rock with a dual-porosity profile like Mount Gambier limestone as chosen in this study. Due to the macro-pore heterogeneity caused by fossil shells and micro-pore heterogeneity from bryozoae (Fig. 1a), the limestone provides an ideal testbed to investigate the heterogeneity problem. Previous studies suggested that rock heterogeneity may enhance compaction localization (Han et al. 2013; Cilona et al. 2014). Nucleation of compaction bands may, however, also be affected by experimental stress boundary conditions (such as boundary friction) (Fortin et al. 2006). Due to the lack of dynamic imaging, the causal relationship between the collapse of the high porosity clusters and nucleation of the compaction bands is still debated (Baud et al. 2015). In this study, we provide direct imaging insights and investigate whether the aforementioned studies can be generalized.

Fig. 1
figure 1

a Scanning electron microscope (SEM) image and b micro-CT tomograph image of the Mount Gambier limestone showing strong micro-structural heterogeneity. c Normalized incremental intruded volume, and cumulative porosity versus pore diameter in semi-logarithmic scale, obtained by mercury intrusion capillary pressure porosimetry (MICP). The three tests show representative results of a total of 9 measurements. The shadow areas correspond pore-size ranges of the two samples to be resolved by X-ray micro-tomography

Full size image

We have therefore conducted two sets of triaxial experiments on carbonate with similar rock composition. The main difference between the two samples was the location of the high porosity and low porosity regions. The first sample features the highest porosity slightly below the center of the sample, and the second features two areas of high porosity exactly in the middle of the sample and slightly above the center. Both samples have dual-porosity profiles with macropores ranging from 10 to 200 microns and nano-micropores ranging from 0.1 to 5 microns. The particular focus of this study is to use digital rock techniques to relate the initial heterogeneity of the sample to the compaction band evolution and record the microstructural changes during the formation of compaction bands.

The samples were imaged using micro-CT at several stages of deformation. The results have been analyzed by digital image correlation (DIC) technique to map the 2 cross-sections at right angles of the axially symmetric strain field.

Methods/Experimental

Sample preparation

The carbonate samples used in this study are from the Gambier Embayment of Otway, the basin of Paleogene-Neogene age in southern South Australia (Bourman et al. 2016). Mount Gambier limestone is a highly porous, fossiliferous carbonate which is originated from extensive colonies of Bryozoa (or lace coral) that flourished on an open-marine shelf, formed over 30 million years ago. Its composition consists of Bryozoa, Foraminifera, and Echinoid spines and plates (Love et al. 1993). X-ray diffraction (XRD) analyses showed that the sample essentially consists of calcite (96%) with a small contribution of quartz, magnesite, and magnetite (accounting for the remaining 4%). An X-ray CT tomograph image of the studied carbonate is shown in Fig. 1b. The image shows material microstructure mainly consisting of carbonate bioclasts and shell fragments. Following the standard triaxial testing protocol (ISRM 2007), specimens were cored with a 2:1 length-diameter ratio from a single carbonate block with approximately 12.7 mm in diameter and 25.4 mm in length. In order to reduce the effect of moisture, all samples were placed in a vacuum oven at 110 °C for 48 h and the temperature was then reduced to room temperature for another 24 h prior to experiments. All experiments were conducted at room temperature.

In order to collect information on the pore size distribution of the sample, mercury injection experiments were conducted using a micrometrics POREMASTER PM-33 (Roshan et al. 2016) (Fig. 1c). The mercury intrusion capillary pressure measurement (MICP) employs a pressurized chamber to force mercury into the porous materials, where the volume of that mercury enters the pores is related to pore volume. Based on the assumption of cylindrical pores, the pore size distribution can be described according to Washburn (1921) by:

$$ P=\frac{-4\gamma \cos \theta }{d} $$
(1)

where P is measured pressure, γ is the surface tension of mercury, θ is the contact angle of mercury, and d is the pore diameter when mercury enters at pressure P. In this study, the mercury surface tension and the contact angle value of 48 N/m2 and 140° were used, respectively (Adamson and Gast 1967). Nine samples were investigated by mercury intrusion experiments (three representative measurements are shown in Fig. 1c) leading to an average interconnected porosity of 51.5%. In order to gain confidence in the porosity measurement by MICP, independent porosity was measured on three core samples using helium porosimetry (TPI-219 Helium Porosimeter) and an average total porosity of 53.1% was obtained, i.e., as the sample is highly porous, the effective and total porosity are assumed identical.

Figure 1a, b shows that the Mount Gambier limestone has a geometrically complex percolative backbone through the macropores (10–200 μm) supported by a solid matrix with a population of micro-pores sizes ranging from 0.1 to 5 μm. Similar results have been reported by Shahin et al. (2019) using Tuffeau de Maastricht limestone (53% porosity) and Ji et al. (2014) using Majella limestone (30% porosity). Two distinct pore throat clusters are evident from Fig. 1c, i.e., incremental intrusion volume divided by the overall intrusion volume showing the pore size distribution. The area under the curves represents the normalized volume of the pores (total normalized volume is equal to 1). MICP measurements indicate that a small volume fraction of smaller clusters, here called micro-pores (< 5 μm), cannot be resolved by micro-CT. This comprises around 10% of the overall porosity.

Experimental set-up

Before running the experiments, a series of triaxial experiments on larger samples were conducted on the specimens at different confining pressure ranging from 1 to 5 MPa (see Appendix, Fig. A1). A confining pressure of 5 MPa was found to trigger pure compaction bands in Mount Gambier Limestone.

In our study, a new X-ray transparent triaxial cell was used (Roshan et al. 2019a) to conduct time-lapse X-ray CT imaging (schematic diagram shown in Fig. 2). The body of the cell was made of PEEK (polyether ether ketone), an X-ray transparent material with enough strength and stiffness to hold required confining pressures. The overall design is similar to a previously designed larger scale triaxial cell (Roshan et al. 2019b) used for the macro-scale experiments. Carbonate specimens with 0.5-inch (12.7 mm) diameter and 1-inch (25.4 mm) length were cored and placed in the micro-CT cell.

Fig. 2
figure 2

Schematic representation of the experimental setup for triaxial compression test for dynamic X-ray micro-CT imaging

Full size image

Axial loads were applied using a 100-ton servo-controlled loading frame with a displacements rate of 3 × 10-3 mm/s (corresponding to a nominal strain rate 1.18 × 10-4/s) at room temperature. The schematics of the experimental setup are presented in Fig. 2. The confining pressure is supplied by hydraulic oil with a Teledyne ISCO Pump (500D) with control up to ± 0.007 MPa. The axial load was logged externally by the servo-controlled loading frame system with ± 0.05 kN accuracy. An external Compression Disk Load cell (LPX-1000, maximum capacity of 1000 kg) was placed underneath the X-ray Triaxial Cell to register a higher sampling rate of the axial load of sample MG-5b. The load cell was connected to a data logger providing also higher accuracy (± 0.001 kN) readings of the axial load compared to the accuracy (± 0.05 kN) registration of the loading frame. The axial strain was measured with reference to the piston movement of the loading frame. The radial strain and volumetric strain were calculated from the obtained 3D tomography images. Axial and radial displacements have an accuracy of ± 0.001 mm and ± 0.02 mm, respectively. All the experiments and scans were conducted at room temperature (~ 25 °C).

In this work, compressive stress and strain are considered positive and tensile stress and strain negative. The volumetric strain is described as εv = εa + 2εr where εa is axial strain and εr is radial strain. The differential stress is defined by q = σ1 − σ3 and the effective mean stress is defined as p = (σ1 + σ2 + σ3)/3 − pp, where σ1 is axial compression (maximum) stress, σ3 = σ2 is the confining pressure (minimum) and pp is pore pressure.

Stress-strain analysis and failure patterns offer a macroscopic understanding of overall deformation in carbonate samples. However, the microscopic processes leading to macroscopic observations of localization also need to be investigated. Therefore, 3D X-ray computed tomography and 2D digital image correlation were used to study the development and evolution of compaction bands from micro-to-macro scale (Ji et al. 2014; Lv et al. 2019).

The specimens were initially scanned in a micro-CT facility (before compression test) and then loaded isotropically to 5 MPa stress while acquiring new series of tomographic images. Triaxial loading was then initiated and at several axial strain levels, the sample was scanned. We conducted two sets of experiments: In the first set of experiments, the sample (MG-5a) was scanned at 0%, 1%, 3%, 6%, 9%, 13%, and 17% of axial strain. Each scan took 1 h to have sufficient time for tomographic imaging where the compression of the sample was kept during 1 h scans under static load. A resolution of 23.8 μm/voxel was obtained. In order to verify the repeatability of our first experiment, a second set of experiments, the sample (MG-5b) was conducted at the same experimental condition but the sample was scanned at 0%, 1%, 3.5%, 6%, and 10% of axial strain; each scan took 2 h, the resolution of 12.5 μm/voxel was achieved. Figure 1c indicates the resolvable pore size (shadow area). Possible creep during the 1–2 h pauses was analyzed using the reconstruction code of the X-ray CT tomography images, which allows a maximum of 5 voxel deformation corresponding to an axial displacement of 100 μm. Within these uncertainty limits, no creep deformation was recorded.

X-ray computed tomography and image analysis

X-ray computed tomography is a non-destructive scanning technique which can map the internal structure of the rock sample with high resolution (Viggiani et al. 2004; Ji et al. 2012; Baud et al. 2015; Pirzada et al. 2018; Viggiani and Tengattini 2019). The technique has been used to get insights into material instabilities (Sheppard et al. 2014) such as observation of the parallel compaction bands in carbonate (Baxevanis et al. 2006; Arzilli et al. 2016) or to obtain localized 3D strain maps (Lenoir et al. 2007; Tudisco et al. 2015; Papazoglou et al. 2017; McBeck et al. 2018; Renard et al. 2019).

In our experiments, all images were acquired at the Tyree X-ray CT facility at the University of New South Wales, Sydney, Australia. The Tyree X-ray CT facility uses a HeliScan micro-CT scanner which was designed to allow easy access and ability to integrate complex flow experiments within the imaging system. The equipment has a GE Phoenix NanoFocus X-ray Tube and high-quality flatbed detector. In our study, X-ray energy and tube current were set to 100 kV and 120 μA, respectively. A double helix trajectory was used for all the scans. We used 0.5 mm stainless steel sample filter notes, and the exposure time was 0.6 s. The resolution difference between sample MG-5a (23.8 μm/voxel, 1 h scan) and sample MG-5b (12.5 μm/voxel, 2 h scan) is because for sample MG-5a, a camera binning mode was used to reduce the scanning time. The X-ray CT images were reconstructed from raw cone-beam X-ray projection data obtained from helix trajectories based on Katsevich inversion formulation (Varslot et al. 2010; Kingston et al. 2016). The workflow of the image analysis includes the following steps: (a) intensity calibration, (b) noise reduction, (c) 3D image registration, (d) segmentation, and (e) 2D digital image correlation.

Intensity calibration

The equilibration of the grayscale intensity was conducted across all images acquired from each sample to obtain higher accurate registration and thresholding segmentation. Because of the effect of filament lifetime, each scan may have a slight X-ray flux difference (less than 100 intensity in our study). In order to calibrate the grayscale intensity, the grayscale of selected regions (large pores, shell fragments, cell body, rubber sleeve, and air) of images of all samples was correlated with each other, and then, differences were corrected for images of each sample using the obtained linear correlation.

Noise reduction

Prior to registration and segmentation, a nonlinear anisotropic diffusion filter (Sheppard et al. 2004) was used for edge-preserving image denoising, using a method similar to that of Frangakis and Hegerl (2001). A key parameter of the anisotropic diffusion filter is the contrast which was set as the intensity difference between the solid phase and pore phase.

3D image registration

The tomography images required registration, because the distance between X-ray source and the X-ray CT transparent cell can be affected by small variations between each scan. The obtained images were first cropped to a cylindrical shape and then registered using a 3D image registration technique developed by Latham et al. (2008). The technique brings two or more images into geometric alignment for further digital image correlation (DIC) and porosity profile analysis. In our study, each scan was registered to the previous scan, once all the images are geometrically aligned. The information from all sources can be correlated to resolve ambiguities which cannot be resolved through the examination of the individual images alone.

Segmentation

Image segmentation is the process of converting a multiphase dataset into two or more phases. In our study, the converging active contours (CAC) method was used for image segmentation (Sheppard et al. 2004; Sheppard et al. 2014). This method uses a combination of the watershed and active contour methods for the segmentation of the grayscale data. The watershed algorithm is figuratively a divide between two adjacent drainage patterns. It is used in image processing to identify changes in grayscale to a thresholding process for segmentation. The active contour method assists in this segmentation process as it allows to identify object outlines from noisy images. The two threshold values for the solid phase and pore phase were chosen based on a layer-by-layer inspection of the 3D tomogram images. In detail, once the segmentation with the initial guesses is conducted, the overall segmented images were compared to that of the original radiographs by the user. This process required a trial and error procedure where convergence was accepted, e.g., when the segmented pores were visually aligned with the original image at one-pixel resolution identified by the user. The true boundary between the solid phase and pore phase in three dimensions then converges by a function which depends on the gradient, intensity, interfacial curvature, and the geodesic distance.

2D digital image correlation

Digital image correlation (DIC) is a modern tool which enables the full-field measurement of displacement and strain fields to obtain detailed information of shear and compaction deformations in geomaterials (Viggiani and Hall 2008; Hall et al. 2010). To investigate the localization and propagation of compaction bands, we used an open-source 2D DIC software package: Ncorr (Blaber et al. 2015), which is implemented in MATLAB. DIC uses image registration algorithms to track relative displacement of material points between a reference image and a deformed image. Ncorr uses a nonlinear correlation optimizer for image alignment. The image alignment algorithm consists of an inversion of the image difference between an identified template to assess the movement and deformation of the template. The template can be either the original undeformed image, or it can be the previous deformation step; thus, total or incremental changes can be detected from the X-ray tomographic slices. The compressive strain εyy and shear strains εxy were extracted from micro-CT images using the Green-Lagrangian strain tensor in an incremental form. Both εyy and εxy are obtained by four displacement gradients.

Results

Characteristics of stress-strain response

Figure 3 summarizes the results of axial stress versus axial strain at 5 MPa confining pressure of two samples. Experiments with a similar high porosity (52%) carbonate, Tuffeau de Maastricht calcarenite, (Baxevanis et al. 2006; Shahin et al. 2019) reported that at 4 MPa confining pressure, a transition from shear bands to compaction bands occurs. At lower confining pressure (< 4 MPa), the specimens developed high-angle shear bands while above 4 MPa confining pressure, compaction bands are formed. This transition of failure mode has been also reported in previous studies (Wong et al. 2001; Bedford et al. 2018; Bouissou et al. 2019).

Fig. 3
figure 3

Axial stress versus axial strain at 5 MPa confining pressure and 2D digital image correlation analysis results of a sample MG-5a and b sample MG-5b (the green color curve of sample MG-5b is the axial strain which was recorded by an externally high-resolution disk load cell). A detailed description of the incremental compressive strain shown in the DIC images will be further discussed in this section

Full size image

In Fig. 3, arrows indicate the points when the axial load was held constant, and the specimens were scanned. The results of 2D digital image analysis were also shown in this figure, but they will be discussed later in detail. The stress-strain curves after yielding showed relatively large and periodic stress drops which are similar to those reported for low strength and porous rocks (Han et al. 2013; Shahin et al. 2019). Such large stress drops are not usually observed in rocks with relatively lower porosity e.g., low porosity carbonate rock (Nicolas et al. 2016) and sandstones (Chen et al. 2017). The results showed that the absolute radial strain is much lower than the axial strain (Fig. 4) showing that inelastic compaction occurred with almost no radial strain (Baud et al. 2017). This is because the highly porous structure of the sample is prone to collapse along the axial direction rather than the radial direction. A significant radial expansion can be seen at 17% axial strain (sample MG-5a).

Fig. 4
figure 4

Effective mean stress (MPa) versus axial strain, radial strain, and volumetric strain of sample MG-5a and sample MG-5b, respectively

Full size image

Microstructure observation

In order to acquire a detailed map of the internal structure, the microstructural analysis was performed on two specimens using full-length X-ray computed tomography. Digital image correlation analysis was performed on both triaxial tests with 5 MPa confining pressure to track the initiation and evolution of localized compaction bands during the experiments. Sample MG-5a was scanned seven times at 0%, 1%, 3%, 6%, 9%, 13%, and 17% axial strain, respectively, and sample MG-5b was scanned five times at 0%, 1%, 3%, 6.5%, and 10% axial strain, respectively.

The aforementioned intensity calibration was performed in a sequential manner for each subsequent test in order to be able to correlate the results for all scans. A linear relationship (R2 ranging from 0.9994 to 0.9998) shows the perfect consistency of our intensity of each scan (Fig. A2 in the Appendix). Furthermore, 3D image registration was conducted based on the calibrated X-ray CT images. The registration process aligns all the scans at one voxel accuracy such that the images overlay with maximum similarity (Latham et al. 2008).

The segmentation process of carbonate rocks presents difficulties due to a large amount of submicron pores (Gray et al. 2016). The segmentation process and results can be subjective because of over-segmentation of pores or grains. From our MICP results, there are some micropores that cannot be segmented because of the limitation of the X-ray CT resolution. Those pores account for 15% porosity of sample MG-5a and 10% porosity of sample MG-5b. This does not negatively affect our analysis as we are not focused on obtaining accurate absolute porosity values but on the evolution of the porosity during the formation of compaction bands. The overall porosity can be estimated from X-ray microtomography (Renard et al. 2009), and the porosity profile can be calculated as the porosity of each layer (1 voxel thickness) along the sample from top to bottom. Figure A3 (in the Appendix) shows the porosity profile of undeformed sample MG-5b when choosing three different intensities ranges.

Figure 5 shows the porosity profile along the sample from the top (0 mm) to the bottom (approx. 26 mm) of sample MG-5a and sample MG-5b. A strong heterogeneity of the porosity can be observed in both samples, ranging from 43 to 58% in sample MG-5a and 38 to 59% in sample MG-5b. In the first scan (scan 1), both samples are undeformed, and therefore, the porosity profile along the sample has an average porosity of 50% (Fig. 5). Significantly larger porosity in both samples MG-5a and MG-5b, with local porosity values up to approximately 60% is seen at different locations near the middle of the sample.

Fig. 5
figure 5

Porosity profile along the z-direction of the specimen extracted from images. a Sample MG-5a and b sample MG-5b. The zones of lowest initial porosity are used as reference points for the registration

Full size image

At 1% axial strain (scan 2), deformation occurs around the top and bottom boundaries where porosity is slightly reduced (Fig. 5). At 3% axial strain (scan 3), the porosity of the top and bottom of the sample significantly reduces from 53 to 33%. Compared to Fig. 3, the stress-strain curve of sample MG-5a shows a relatively larger drop before scan 3 and sample MG-5b shows yielding and slight stress fluctuation before scan 3 followed by two large stress drops between S3 and S4. Comparing with Fig. 3, we observe that the first compaction band corresponds to the yielding or first drop of the axial stress in the stress-strain response. The first compaction bands hence form on top and bottom boundaries of the sample. From scan 3 to scan 4, the axial stress of sample MG-5a increases without any further drop (Fig. 3), while the reduced porosity zones grow towards the middle of the sample and the top and bottom boundary porosity reduces to 28%. However, the axial stress of sample MG-5b shows two large stress drops corresponding to two drops in the porosity profile. These are initiated at the zone of highest local porosity, indicating that the mechanism for compaction band formation is macropore collapse of large porosity zones.

Between scan 4 and scan 5, a large drop in axial stress of sample MG-5a shows similar stress drop of sample MG-5b which occurred between scan 3 and scan 4, the porosity profile of sample MG-5a also shows a significant drop in the middle of the sample having the highest initial porosity of 58%, i.e., where the porosity reduces from 58 to 35% and a 3 mm thickness compaction band is formed.

An interesting observation can be made between scans 4 and 5 in Fig. 5b showing that in the sample MG-5b, the region of highest porosity continues to collapse although the local porosity reaches lower values than the bottom of the sample between 15.5 mm and 22.8 mm. We attribute this to the loss of coherency in the crushed zone. The remainder of the sample still bears structural integrity and can support the applied load. A similar observation can be made for scan 5 to scan 6 in sample MG-5a (Fig. 5a), where the axial stress increases gradually (Fig. 3) and the region from 15.4 to 20.7 mm of the sample loses structural coherence while the top of the sample from 3.4 to 11.1 mm remains comparatively intact. Note that the gradual increase in axial stress corresponds to a diffuse growth of compaction in this section. A sudden collapse of the entire sample is observed in scan 7 of sample MG-5a through a dramatic decrease in the porosity throughout the entire sample where porosity reduces to an average of 20% and the majority of macropores have collapsed. At this stage, the load-bearing framework of the entire sample is lost and the macroscale behavior transition into that of a granular fluid which is not the subject of this paper.

Microscopic view of localization and propagation of compaction bands

As previously mentioned, an open-source 2D digital image correlation (DIC) code: Ncorr (Blaber et al. 2015) was used to analyze the localization and propagation of compaction bands. Figures 6 and 7 present the strain along z-direction on XZ and YZ planes, i.e., representative tomograph slices are shown at right angles to each other. As the mechanical load is monotonic and compression sign is positive, only the compressive strain is developed in the sample along z-axis on both XZ and YZ planes.

Fig. 6
figure 6

a Incremental shear strain εxy and b incremental compressive strain εyy along the z-direction (intersecting XZ plane and YZ plane) from DIC analysis. c Porosity reduction profile derived from segmentation images analysis of sample MG-5a

Full size image
Fig. 7
figure 7

a Incremental shear strain εxy and b incremental compressive strain εyyalong the z-direction (intersecting XZ plane and YZ plane) from DIC analysis. c Porosity reduction profile derived from segmentation images analysis of sample MG-5b

Full size image

No shear or compressive strain localization appears at 1% axial strain for both specimens (Figs. 6 and 7). However, strain localization is seen for all other scans immediately following the yield point, illustrating that the yield phenomenon is promoted by localized pore collapse. The localization appears at the top and bottom of the specimen and then migrates towards the middle. The nucleation of compaction on the bottom and top platens needs to be treated with caution as this is affected by the interaction of the sample with the platens. With respect to the shear strain profile in Figs. 6b and 7b, only a small boundary effect is noted. Overall shear strains are negligible compared to compressive strain suggesting that the observations are true compaction bands.

When focusing on sample MG-5a, the porosity profile (Fig. 5a) supports the DIC observations where porosity across the sample starts changing considerably after yield. When progressing towards larger strains (17%), the sample porosity drops to 24.8%, where DIC results show significant complex failure patterns at this stage. The DIC analysis for such high strain regions breaks down as it resolves only the large-scale pore structure. When the structure collapses, an error is reported and the areas are marked by showing the background tomographic images. These zones define regions of strain larger or smaller than the scale bar due to errors in the DIC measurements caused by the local collapse of reference points.

Using the stress-strain curve from Fig. 3 with the porosity profile in Fig. 5a, we infer that each drop in the stress profile after yielding is likely to correspond to the nucleation of a new compaction band while the strain hardening between stress drops represents the growth of the already existing compaction bands. A similar conclusion is also reported by Baud et al. (2004) inferred from experimental observations and Das et al. (2013) through numerical analysis. The stress drop is therefore chosen as an indicator of compaction localization.

A series of 2D tomography slices superposed with porosity profile on the YZ plane showing the structural changes inside and outside the compaction bands are presented in Fig. 8 i.e. scan 3 to scan 7 of sample MG-5a and scan 1 to scan 5 of sample MG-5b were shown. In order to facilitate the correlation of the images, the arbitrary kinematic reference system was chosen to be the undeformed section of the sample. Scan 3 to scan 7 of sample MG-5a and scan 1 to scan 5 of sample MG-5b were shown here. The red boxes indicate the compaction zones. Almost all of the collapsed pores are located in the compaction band area, and there is no observable collapse outside of the compaction zone. This is additionally supported by the porosity profile (Fig. 5) and strain map (Figs. 6 and 7). The growth of the existing compaction bands leads to a reduction in porosity also seen from the porosity profile (Fig. 5).

Fig. 8
figure 8

2D slices of X-ray tomographs of scan 3 to scan 7 on YZ plane shown on the porosity profile of each scan of a sample MG-5a and b sample-5b, respectively, black color represents pores and white color represents grains. Red boxes indicate compaction zones. Note scan 7 of sample MG-5a shows some residual large-scale porosity at the top of the sample (compare with Fig. 5a)

Full size image

A compaction band can be described as discrete or diffuse. A “discrete band” is defined as a few grains thick tabular structure with cumulative strain accommodated by the initiation of additional discrete bands (Wong et al. 2001; Baud et al. 2004). A “diffuse band” is much thicker and accommodates the cumulative strain by lateral propagation of damage (Holcomb and Olsson 2003). Figure 8a shows a single discrete compaction band developed in the middle of the sample MG-5a between scan 4 and scan 5 and few diffuse compacting features in the lower part of the specimen picked up by the DIC analysis (Fig. 7). The discrete compaction band forming in this stage can be compared to experiments on notched sandstone specimens. However, the single discrete compaction band has a thickness of approximately 1180 μm compared to the large pore size 60~80 μm which is thicker than the sandstone where a thickness of 2 or 3 times sand grain diameter was reported (Vajdova and Wong 2003; Tembe et al. 2006). Note that in the tomograph shown in Fig. 8a, this observation is not clear and we assume that the diffuse compaction zones picked up by the DIC may not propagate through the entire sample at this early stage. A similar observation, where both discrete and diffuse compaction bands co-exist in Berea sandstone was previously reported by Baud et al. (2004).

A representative image of the collapsing macropores obtained from the micro-CT images is shown in Fig. 9 in two selected windows: (a) 310 voxels × 232 voxels and 3.90 mm × 2.92 mm, and (b) 300 voxels × 193 voxels and 3.77 mm × 2.43 mm from XY plane of sample MG-5b at scan 3 (left) and scan 4 (right). Figure 9a shows large pores which are undeformed at scan 3, when increasing the axial strain from 3 to 6%, the zone suffered significant deformation and the structure of the fossil shell was broken into fragments while surrounding grains were crushed into powder and pushed into the macropores at scan 4. Figure 9b shows that some mesopores which have stronger structure still retain their mechanical integrity. Thus, macropore collapse leads to a reduction of the pore throats and acts as a flow barrier.

Fig. 9
figure 9

Two selected 2D tomography image windows a 310 voxels × 230 voxels and 3.90 mm × 2.92 mm, and b 300 voxels × 193 voxels and 3.77 mm × 2.43 mm from XY plane of sample MG-5b at scan 3 and scan 4 shows the collapse of macropores. The figure on the left shows the 3D volume rendering of the tomography image and the red line indicates the layer of the chosen window. The weaker structure of the fossil shells which is shown by the dashed boxes is broken into several angular fragments and crushed calcite powder filled the pore space

Full size image

Discussion

We have shown that the dynamic micro-CT analyses can provide direct insights into the mechanism of formation and growth of compaction bands with special emphasis on quantifying the role of heterogeneity to trigger the initiation and evolution of compaction bands. Digital rock techniques such as micro-CT images and digital image analysis (DIC) can therefore provide insight into deformation and its localisation in porous rocks. Our main findings are summarized below.

Effect of heterogeneity of local porosity on nucleation and growth

The nucleation of compaction bands in such a high heterogeneous carbonate rock is still an open question. Previous studies showed the initiation and propagation of compaction bands depend on grain crushing, pore collapse, and breakage of the cementation (Vajdova and Wong 2003; Charalampidou et al. 2011; Das et al. 2011; Han et al. 2013). Previous laboratory experimental studies (Olsson 2001; Baud et al. 2004) observed that the nucleation of compaction bands from the top and bottom boundaries which can be explained as the stiffness mismatch between the sample and the testing device. This causes a local stress heterogeneity on the platens triggering the formation of the first compaction band. Fortin (2006) also pointed out that the material heterogeneities (e.g., high local porosity cluster) can collapse first and trigger compaction bands in Bleurswiller sandstone. Numerical breakage mechanics simulations studies also showed that the boundary stiffness mismatch and weak points in a relative homogeneous sample can trigger compaction bands (Das et al. 2013; Das and Buscarnera 2014).

We have confirmed insights on the role of local porosity variation for the nucleation of compaction bands reported in dynamic CT analyses of a lower porosity (26%) limestone (Huang et al. 2019) which has local porosity heterogeneity near the top of the sample. Their sample was much smaller (5 mm diameter and 10 mm in height) and local heterogeneities therefore played a more significant role. In particular, the effect of the top and bottom platen is more prominent for smaller geometries. We conclude that the propagation of localized compaction bands initiated from both ends of the sample first can be considered a boundary effect attributed to the aforementioned stiffness contrast of sample and testing machine. We have, however, clearly shown that the region of highest porosity in the sample is the nucleation site of the first major compaction band in two repeated experiments (Cilona et al. 2014). This together with the previous work provides strong evidence for the role of local microstructural heterogeneity controlling this mechanical instability. However, we have also shown that the stochastic heterogeneity control reverts back to a continuum formulation as further compaction grows from this initial seed. This results in a diffuse growth of the compaction band from the initial seed.

Link between diffuse and discrete compaction bands

We were able to show that the new triaxial cell in conjunction with in situ X-ray CT imaging analysis can directly quantify the location and thickness of the compaction bands during deformation, i.e., within the method resolution, confirming the results of the thin section method. We noted a transition from discrete to diffuse compaction band formation mechanism after the discrete band developed in scan 4 (sample MG-5a). This is similar to observations in some sandstones, for instance, where a combination of diffuse compaction zones as well as discrete compaction bands was found (e.g., Baud et al. (2004)). While the findings from sandstone seem to also apply to the extreme end-member carbonate, the analysis suggests that the switchover from one mode to the other can happen at any particular moment in the experiment. We suggest that a phenomenological observation for discrete compaction bands is a signature large stress drop (Figs. 3a, 5a, and 6b, scans 4–5; Figs. 3b, 5b, and 7b, scans 3–4) while the growth of continuous bands are supported by minor drops followed by strain hardening. This can be seen when comparing Fig. 3a with the Fig. 5a (scans 5–6) without large stress drops. Similar results have been observed by Klein et al. (2001) and Baud et al. (2004) who showed that discrete bands are associated with overall strain hardening punctuated by episodic stress drops. The magnitude of the stress drops has been reported to be related to the ratio between the thickness of the bands and the length of the sample. Small bands in large samples are not inferred to produce noticeable stress drops (e.g., Baud et al. (2004)).

The initial nucleation of discrete compaction bands may be attributed to boundary friction on the top and bottom platen; however, a second more prominent band nucleates in the local region of highest porosity near the center of the sample. Note the different locations of the highest porosity zone in both samples trigger the site of initial localization and subsequent diffusive growth. The discrete bands are argued here to act as the initial seed for a wave-like propagation of discrete micro-compaction bands away from the initial seed. These micro-compaction events are associated with minor stress drops, and we propose that they increase the thickness of the compaction bands in a linear fashion. In order to verify this hypothesis, the technique presented here would need to be carried out in the Synchrotron X-Ray beam to achieve the time resolution required for the detection of the speed of the proposed propagating micro-compaction bands.

The effect of boundary friction on the top and bottom platen is difficult to clearly identify. However, the incremental shear strain on the boundaries (see Figs. 5b and 6b) is not affected, and the compaction band that formed near the top and bottom of the sample is abandoned after scan 4 in sample MG-5a and after scan 3 in sample MG-5b (Fig. 5).

Conclusions

We have presented a comprehensive analysis of the formation of compaction bands in highly porous carbonates using advanced image analysis techniques comprising triaxial testing and analysis, micro-CT tomography, and DIC analysis to perform a 4D analysis of the principal mechanism and their evolution. The techniques were shown to mutually reinforce the main findings offering different perspectives of microscopic and macroscopic observations.

Given the comprehensive analyses of the formation of compaction bands that have been performed on sandstones and siliciclastic rocks and the affirmative conclusions of the similarity between earlier studies and the extreme end-member heterogeneous carbonate in this study, as well as the recently published micro-CT work (Huang et al. 2019; Shahin et al. 2019), we may attempt to interpret the phenomenon of compaction band formation using the continuum theories that have been proposed. Issen and Rudnicki (2001) summarized the incipient condition of compaction bands in porous rocks based on a rate-independent constitutive model (Rudnicki and Rice 1975; Rudnicki 1977; Olsson 1999; Issen and Rudnicki 2001). Yield cap failure criteria have been used to determine the onset of instabilities, both in shear and compression (Rudnicki 2004; Grueschow and Rudnicki 2005). Katsman and Aharonov (2006) presented a new damage elastoplastic constitutive model taking into account the formation of the compaction bands under a variety of boundary conditions. Das et al. (2011) proposed a constitutive model based on the concept of grain crushing, which uses a breakage measure as an additional internal variable to the continuum model to track the evolution of grain size distribution. Chemenda and his colleagues (2011; 2012) claimed that the onset of compaction bands is triggered by constitutive instabilities resulting from deformation bifurcation rather than anti-crack models of linear elastic fracture mechanics. A recent model emphasized the rate-sensitivity of the instability and proposed a new mechanism forming compaction bands and dilation bands (Veveakis and Regenauer-Lieb 2015; Regenauer-Lieb et al. 2016).

All of the abovementioned models fail to explain the observation of growth of diffuse compaction zones as well as the transition from discrete to diffuse compaction band formation reported in our experiments. This calls for an extension of the existing theories to capture the dynamics reported in our experiments. On the flipside, future experimental work needs to adequately address how the strain data can be used for the prediction of failure. This requires innovations for the direct measurement of the behavior of both stress and strain fields at the mesoscale of the localization band.

Triaxial tests show a transitional phase of nucleation of first compaction bands on the top and bottom platens which can be attributed to boundary effects and porosity heterogeneity (Huang et al. 2019). These boundary zones may either progress smoothly into the center of the sample or as in our experiments be aborted and be followed by nucleation of compaction bands on zones of highest porosity. This observation of nucleation of compaction bands on the highest porosity zones shows the strong effect of material heterogeneity to control the site of this mechanical instability also controlling the sites of continued diffusive growth. We therefore conclude that there are two components in the nucleation of compaction bands: (1) the macroscopic boundary conditions (size of the sample), friction on the platen, and its mechanical response as proposed in the incipient continuum mechanics models and (2) the local heterogeneity of the sample (large porosity zones or bedding planes) may become the dominant feature for initial localization (Cilona et al. 2014). This information may also help in identifying reasons for a switch in compaction band formation from discrete to continuous zone growth. Larger porosity values trigger the initial pore collapse followed by a diffuse growth of compaction bands due to loss of mechanical strength inside the collapsed zone. We have also identified in our experiment a second phase transition from macroscopic solid behavior to macroscopic fluid-like behavior when the majority of macropores have collapsed.