Europe’s largest gas field, the Groningen field (the Netherlands), is widely known for induced subsidence and seismicity caused by gas pressure depletion and associated compaction of the sandstone reservoir. Whether compaction is elastic or partly inelastic, as implied by recent experiments, is a key factor in forecasting system behavior and seismic hazard. We sought evidence for inelastic deformation through comparative microstructural analysis of unique drill core recovered from the seismogenic center of the field in 2015, 50 yr after gas production started, versus core recovered before production (1965). Quartz grain fracturing, crack healing, and stress-induced Dauphiné twinning are equally developed in the 2015 and 1965 cores, with the only measurable effect of gas production being enhanced microcracking of sparse K-feldspar grains in the 2015 core. Interpreting these grains as strain markers, we suggest that reservoir compaction involves elastic strain plus inelastic compression of weak clay films within grain contacts.

Hydrocarbon production reduces reservoir fluid pressure and increases the in situ (effective) stress, frequently leading to surface subsidence and induced seismicity (Geertsma, 1973; Yerkes and Castle, 1976). Subsidence is caused by one-dimensional (1-D) vertical compaction of the reservoir rock, whereas seismicity reflects coupled stress buildup at preexisting faults (Segall, 1989; Zoback, 2007). Reservoir compaction is widely assumed to be poro-elastic and, hence, reversible (Segall, 1992; Zoback, 2007). However, experiments on sandstones show that irreversible deformation contributes continuously during loading in the assumed poro-elastic range of 0.1%–0.2% strain (Bernabé et al., 1994; Shalev et al., 2014; Hol et al., 2018; Pijnenburg et al., 2018, 2019a), with sharp “plastic” yield occurring only at strains >1%–2%; i.e., far beyond those relevant for production-induced compaction (<0.5%). Including constitutive relations describing small-strain inelastic deformation may substantially modify forecasts of production-induced stress evolution (Pijnenburg et al., 2019a), subsidence, and induced seismicity. However, direct observations of the mechanisms responsible for irreversible, production-induced deformation are lacking.

One of the world’s best-known gas fields showing production-induced compaction, subsidence, and seismicity is the vast Groningen field, the Netherlands (Fig. 1A). The reservoir is the 100–350-m-thick Slochteren Sandstone (Gaupp and Okkerman, 2011), located at ∼3 km depth, where the overburden pressure is ∼65 MPa (Fig. 1B; Van Eijs, 2015). Conglomeratic and mudstone facies dominate the south and north of the field, whereas the most porous, central part is typified by fluvial and eolian sands (de Jager and Visser, 2017). Since gas production started in the 1960s, gas pressure has fallen from 35 to ∼8 MPa, and subsidence, driven by the increase in vertical effective stress, has reached ∼33 cm in the center of the field (Fig. 1A). This implies vertical reservoir compaction strains of 0.1%–0.3%, consistent with in situ fiber-optic strain measurements (Cannon and Kole, 2017). Since ca. 1990, production has been accompanied by seismic events within the reservoir (Spetzler and Dost, 2017), causing property damage and public concern, notably in 2012, when the largest-magnitude earthquake (local magnitude ML 3.6) occurred near Huizinge village. In 2015, field operator Nederlandse Aardolie Maatschappij BV (NAM) launched an unprecedented drilling operation at the Zeerijp ZRP-3a well (Figs. 1A and 1B) to retrieve reservoir core from the highly depleted, most seismogenic, central part of the field.

We report an investigation of Slochteren Sandstone samples from the ZRP-3a well, aimed at searching for evidence of small-strain (0.1%–0.3%), inelastic deformation processes. We used Slochteren Sandstone core recovered before depletion (1965), from the nearby Stedum SDM-1 well, as an undeformed benchmark. Potential mechanisms of inelastic reservoir compaction include (1) fracture of grains forming the load-supporting rock framework, (2) grain-scale dissolution-precipitation, and (3) processes operating within grain contacts, like asperity crushing/dissolution, clay film deformation, and slip (Spiers et al., 2017). We tested which of these may have operated by investigating the microstructure of reservoir core recovered from multiple depths in both wells (Fig. S1 in the Supplemental Material1), thus deliberately sampling the lithological variations that characterize the central part of the field.

We studied sandstones samples from the SDM-1 and ZRP-3a cores using techniques ranging from visual inspection to optical and backscattered electron (BSE) microscopy, X-ray mapping, cathodoluminescence (CL) imaging, electron backscatter diffraction (EBSD) mapping, and image analysis.

Qualitative Analysis

Depleted and undepleted sandstone samples were indistinguishable at hand-specimen scale. A single SDM-1 core plug, laboratory-deformed to ∼0.2% elastic plus ∼0.2% inelastic strain, under conditions simulating depletion at well ZRP-3a, was also visually indistinguishable from its SDM-1 precursor (section S1 in the Supplemental Material; Fig. S2). To investigate grain-scale microstructures, 28 polished sections were prepared for optical and BSE imaging, using core fragments taken from 13 near-evenly distributed depth intervals (15–20 m) in each well, and using the laboratory-deformed sample and its undeformed counterpart for comparison (Fig. S1; Tables S1 and S2). All samples showed millimeter- to centimeter-scale sedimentary layering, enabling the effects of varying porosity, grain size/shape, and mineralogy on deformation signatures to be investigated. All were dominated by rounded to subangular quartz grains (≥70%), with minor quantities of K-feldspar (KFs; ≤15%), Ca-Mg-Fe carbonates (3–6%), lithic fragments (3–10%), and clay (illite, kaolinite; 0.5–5.5%; Fig. 2; Fig. S3), with a clastic grain size of 100–200 µm (Fig. 2; Tables S1 and S2). Carbonates were mainly present as authigenic cement in specific sedimentary layers. The scattered feldspar and lithic grains were frequently “corroded” and porous. Clay platelets were widespread, forming vermicules in pores and thin layers coating grains and grain contacts (average thickness 2–5 µm; Waldmann and Gaupp, 2016). Intergranular films mainly consisted of preburial, tangential clay, with postburial, radial clay present on pore walls. Indented and sutured grain contacts indicating diagenetic pressure solution were widespread. Intergranular and transgranular cracks were rarely observed, but intragranular cracks were widespread, regardless of mineral type, depth interval, or source well. Both fresh and reacted feldspars frequently showed multiple intragranular cleavage cracks, made visible by dissolution. From qualitative inspection of optical or BSE micrographs, no distinction could be made between the depleted, undepleted, or laboratory-deformed samples.

Quantitative Microstructural Analysis

To assess the role of grain fracturing during reservoir depletion, porosity (ϕ), grain-size distribution (GSD), and intragranular crack density (crack count per unit grain area, ρcr) were quantified within microstructurally distinctive sedimentary domains measuring ∼10–190 mm2, as observed in section-scale BSE mosaics (section S2; Fig. S4; Tables S1 and S2). We identified tens of thousands of intragranular cracks within ∼50 domains in each of the SDM-1 and ZRP-3a core samples, plus seven in the laboratory-deformed sample. No correlations were found between ρcr and mean/median grain size, or with the standard deviation, skewness, or kurtosis of the GSD, in any of the materials studied (Fig. S5). Crack densities ranged from ∼1 to 30 mm–2 for porosities of 4–29%, with the ρcr-ϕ data for the entire suite of samples studied falling in a band (Fig. 3A), suggesting a tendency for crack density to increase with increasing porosity. However, the data for the depleted, undepleted, and laboratory-deformed samples were indistinguishable (Fig. 3A, inset). This suggests that grain failure did not play a significant role in producing inelastic strain in the ZRP-3a core versus the undepleted SDM-1 core, or in the laboratory test on SDM-1 material.

Crack density analysis combined with electron dispersive X-ray mapping revealed that KFs grains frequently showed multiple intragranular (cleavage) fractures (Figs. 2A and 3B; Fig. S3). To avoid bias due to varying KFs grain area forumla per sample, for each domain we determined the fraction of cracks cutting KFs grains forumla versus the KFs area fraction forumla (Fig. 3C). For both the depleted and undepleted samples, we found μ1 > μ2, implying that KFs grains were preferentially cracked compared with quartz and other minerals. Moreover, two-sample hypothesis testing of the ZRP-3a versus SDM-1 sample population means showed statistical discernibility at a 94% confidence level (section S4). This suggests that KFs grains in samples from the depleted core are preferentially cracked compared with samples from the undepleted core (Fig. 3C). Coupled with local intergranular displacements, this could “nucleate” bulk inelastic deformation. However, in view of the low abundance of KFs compared with quartz (∼15%), and their spatially scattered distribution in the reservoir rock, this is unlikely. More likely, feldspar grains act as passive markers indicative of deformation of the surrounding grain framework. With the grain strength of quartz being higher than feldspar (Hangx et al., 2010), especially when the latter is corroded by reaction, even small, elastic strains within the surrounding quartz-grain framework should lead to fracturing of the feldspar grains.

The possibility that enhanced stresses at quartz grain contacts due to gas depletion caused reactivation of healed intragranular cracks, or of pressure-dissolved contacts, was also investigated. Using high-resolution panchromatic CL imaging, we scrutinized ZRP-3a and SDM-1 samples for quartz overgrowths and (reactivated) healed fractures (Figs. 2C and 2D; Figs. S6 and S7). We found no evidence for more, or more recent, diffusive mass transport in depleted (ZRP-3a) versus undepleted (SDM-1) samples, nor for statistically significant differences in the reactivation of healed fractures (Table S3). In quartz, high compressive stresses can trigger Dauphiné twinning (Tullis, 1970), offering a potential indicator of enhanced grain-contact stresses (Mørk and Moen, 2007). EBSD mapping of Dauphiné twin (DT) boundaries showed that DTs were present in all samples analyzed, at grain contacts and within grain interiors (Figs. 2E and 2F; Fig. S8). The DT density (ρDT), defined as the ratio of the total DT boundary length to the total quartz grain area, measured 10.6 mm–1 in depleted core, 9.44 and 16.8 mm–1 in undepleted core, and 14.5 mm–1 in the laboratory-deformed sample (Table S4). We infer that (simulated) depletion was insufficient to measurably enhance Dauphiné twinning.

Our results imply that neither fracturing nor dissolution of the quartz framework grains played a significant role in causing production-induced compaction at the center of the Groningen field. The only measurable microstructural difference between samples from the undepleted SDM-1 and depleted ZRP-3a wells is the breakage of KFs. We suggest that preferred fracturing of these scattered, corroded, cleavage-prone, and hence weaker, KFs grains occurred as a result of the vertical poro-elastic compaction of the quartz-grain framework that occurred during reservoir depletion. Our interpretation is that the preferentially fractured KFs grains persisted as a qualitative record of the elastic strain accommodated by the rock framework, without significantly contributing to permanent deformation. At the same time, it is important to consider the effects of the weak clay films that coat most detrital grain contacts (Fig. 2; Fig. S3). Because the Groningen reservoir lies close to its maximum burial depth (Gaupp and Okkerman, 2011), and because long-term clay consolidation depends mainly on effective normal stress (Brown et al., 2017), intergranular clay films compacted during burial would almost certainly deform further as the effective vertical stress increased from ∼30 to 57 MPa during depletion of the Groningen reservoir. As demonstrated by Pijnenburg et al. (2019b), such stress changes are sufficient to inelastically compact individual 2–5-µm-thick intergranular clay films in the Slochteren sandstone by ∼4%, i.e., by 0.1–0.2 µm. Coupled with the intergranular sliding (clay film shearing) required by the in situ condition of 1-D vertical strain, this implies reservoir compaction strains of 0.1–0.2% for typical detrital grain sizes of 100–200 µm (Pijnenburg et al., 2019b). While these compaction strains match field subsidence and in situ strain data well, they are virtually impossible to detect in the postmortem microstructure, as there are no internal strain markers either at the grain or clay film scales.

Combining (1) the irreversible deformation observed in small-strain experiments simulating depletion (Hol et al., 2018; Pijnenburg et al., 2018, 2019a) and attributed to clay film deformation (Pijnenburg et al., 2019b), (2) the widespread occurrence of intergranular clay coatings in the Groningen reservoir (Waldmann and Gaupp, 2016), and (3) the burial history of the reservoir (Gaupp and Okkerman, 2011), we infer that elastic deformation of the Slochteren Sandstone plus inelastic deformation by intergranular clay film compression and shear offer the most viable explanation for the strains recorded by cracked feldspars in the depleted core (Fig. 4).

To date, studies addressing inelastic compaction of porous rock have focused on the microphysical mechanisms controlling macroscopic plastic yield at strains far in excess (>1–2%) of the strains occurring in hydrocarbon reservoirs (Brantut et al., 2013; Menéndez et al., 1996; Wong and Baud, 2012). These studies frequently show strain localization accommodated by grain fracturing and associated grain-scale rearrangement. For the case of small-strain, production-induced, seismogenic compaction of the Groningen and similar reservoirs, our study suggests that deformation at grain contacts, involving micrometer-scale clay film compaction and shear, may be a key factor (Fig. 4). The compressive strength of a wide range of reservoir sandstones was recently shown to decrease with increasing clay content (Heap et al., 2019), supporting the proposed mechanism. Physics-based models aiming to forecast induced subsidence and seismicity (Bourne et al., 2014; Zbinden et al., 2017; Buijze et al., 2019), as in the case of the Groningen field, should aim to incorporate constitutive models describing inelastic processes, such as the intergranular clay film deformation inferred here, because the impact on horizontal stress evolution is substantial (Pijnenburg et al., 2019a, 2019b).

We thank J. van Elk, D. Doornhof, C. Visser, D. Doran, F. Marcellis, T. Wolterbeek, R. Keijzer, E. Gutiérrez, S. Matveev, S. Waldmann, and L. Bik. D. Healy, R. Gaupp, and an anonymous reviewer provided helpful comments. The work was funded by Groningen field operator Nederlandse Aardolie Maatschappij BV (NAM) (agreement UI49358). European Research Council starting grant SEISMIC 335915 supported Hamers.

1Supplemental Material. Detailed description of the materials and methods employed in this study, Figures S1–S8, and Tables S1–S4. Please visit https://doi.org/10.1130/GEOL.S.13262861 to access the supplemental material, and contact [email protected] with any questions.