New approaches in experimental research on rock and fault behaviour in the Groningen gas field

Abstract

This paper describes a research programme recently initiated at Utrecht University that aims to contribute new, fundamental physical understanding and quantitative descriptions of rock and fault behaviour needed to advance understanding of reservoir compaction and fault behaviour in the context of induced seismicity and subsidence in the Groningen gas field. The NAM-funded programme involves experimental rock and fault mechanics work, microscale observational studies to determine the processes that control reservoir rock deformation and fault slip, modelling and experimental work aimed at establishing upscaling rules between laboratory and field scales, and geomechanical modelling of fault rupture and earthquake generation at the reservoir scale. Here, we focus on describing the programme and its intended contribution to understanding the response of the Groningen field to gas production. The key knowledge gaps that drive the programme are discussed and the approaches employed to address them are highlighted. Some of the first results emerging from the work in progress are also reported briefly and are providing important new insights.

Introduction

As described in numerous previous reports, and in several contributions to the present special issue, the Groningen gas field is one of the largest onshore gas fields in the world. This giant field occupies a 30km by 30km region of the northeast Netherlands. Discovered in the 1950s, it has been producing gas since 1963 and since that time has contributed a vast revenue to the Dutch economy. The gas-bearing reservoir formation is the Slochteren Sandstone, part of the Permian Rotliegend series (Fig. 1A). The formation consists of aeolian, alluvial and fan deposits, some 100–250m in total thickness, located at a depth of c. 3km below the surface (Geluk, 2007). It is overlain by the Ten Boer Claystone unit (10–50m thick) and by the basal anhydrite–carbonate unit of the Zechstein evaporite series (also c. 50m thick), which forms a thick caprock sequence. It rests unconformably on organic-rich shales, siltstones and coals of Upper Carboniferous age, which constitute the source rocks for the gas. For details of the stratigraphy of the Groningen field, the reader is referred to the reports published by the field operator, NAM (Nederlandse Aardolie Maatschappij), and to other papers in this volume. The vertical stress within the reservoir is c. 65 MPa, with the initial gas pressure of c. 35 MPa now having fallen to c. 8 MPa after some 55 years of gas production. The reservoir temperature is c. 100°C and its lower portion is charged with highly saline brine. The in situ conditions are therefore relatively hot, reactive and characterised by high effective stress.

Fig. 1. (A) Schematic illustration of the stratigraphy of the Groningen gas field, with enlargement of a faulted and internally fractured section of the Slochteren Sandstone reservoir and immediate over- and underburden units. (B) Schematic diagram illustrating vertical compaction of the reservoir sandstone due to an increase in vertical effective stress caused by depletion of the pore fluid pressure accompanying gas production. Compaction of the reservoir sandstone by poro-elastic and permanent inelastic processes leads to deformation of the surrounding rock mass and to surface subsidence. Differential compaction occurring across faults in the reservoir, notably where more rigid over- and underburden formations are juxtaposed against the compacting reservoir, can lead to shear stress concentrations across the faults, which may in turn lead to fault reactivation and induced seismic rupture.

After an initial delay of 5–6 years, gas production from the Groningen field has been accompanied by more or less steady subsidence at the surface, reaching a little over 30cm in the centre of the field in 2016 (NAM, 2016). This reflects compaction occurring in the reservoir system as the vertical effective stress increases. Since the 1990s, an increasing frequency and magnitude of induced seismic events has also been recorded, notably in the central region of the field where subsidence and porosity is greatest, culminating in 2012 in the largest event to date, the M _L=3.6 Huizinge earthquake (Dost & Kraaijpoel, 2013). This caused significant public concern. Since that time, the operator, NAM, has invested heavily in extending seismological instrumentation of the field, and it is becoming increasingly clear that the events recorded tend to occur within or close to the reservoir interval (Spetzler & Dost, 2017) and near, or on, the numerous pre-existing (mainly normal) faults that cross-cut the reservoir (Van Eijs et al., 2006; Dost & Haak, 2007). There has been much speculation as to whether increased production rates lead, after some delay, to increased seismicity caused by time-dependent processes operating in the reservoir and surrounding system (van Thienen-Visser & Breunese, 2015; van Thienen-Visser et al., 2015; Paleja & Bierman, 2016). In the last few years, production has accordingly been reduced in the central part of the field, where larger events have occurred, in an attempt at mitigation. The basis for this strategy is questionable, however, as data are limited and the physics of how the subsurface responds to production is poorly understood.

In broad terms, it is widely believed that the principal processes occurring in depleting gas reservoirs involve compaction of the reservoir rock, as the vertical effective stress increases. This leads to non-uniform deformation of the generally faulted and mechanically heterogeneous reservoir and over/underburden system, and to an associated non-uniform stress field (e.g. Zoback, 2010; Sone & Zoback, 2014;). Gradients in compaction and associated deviations in shear and normal stress in the neighbourhood of pre-existing faults and other heterogeneities can then trigger fault slip or other rock failure phenomena, and potentially induced seismicity (Roest & Kuilman, 1994; Mulders, 2003). Particularly relevant are differential compaction and associated shear stresses generated across faults that juxtapose compacting reservoir rock against relatively rigid over- or underburden units (Mulders, 2003).

Geomechanical modelling of the evolution of reservoir systems, and of the initiation and propagation of quasi-static and dynamic fault rupture and seismic energy release, is progressing rapidly (Ida, 1972; Andrews, 1976a,b; Madariaga, 1976; Tang & Kaiser, 1998; Settari & Walters, 2001; Mulders, 2003; Rice et al., 2005; Cappa & Rutqvist, 2011; van den Bogert, 2015; Lele et al., 2016; Urpi et al., 2016). However, the physical basis for such work is as yet poorly constrained and uncertainties remain large and often unknown. The reliability of geomechanical modelling and of our understanding of systems such as the Groningen reservoir, depends to a large extent on having a physically based (i.e. mechanistic) understanding and quantitative description of the mechanical behaviour of the reservoir rock, the surrounding rock types and of faults. Obtaining such descriptions, in particular of reservoir compaction and fault friction behaviour, is crucial for modelling stress–strain gradients developing across faults in the reservoir and the coupling with induced fault slip. Only then can trends in system behaviour and seismicity be understood and modelled with any predictive capacity, in the manner needed to steer production strategies. At present, it has to be admitted that there is insufficient scientific basis to evaluate how the field will respond to changes in production rates, or even to stopping production.

Much work has been done on rock physical properties (e.g. Paterson & Wong, 2005; Wong & Baud, 2012; Brantut et al., 2013) and fault mechanical properties (Byerlee, 1978; Tembe et al., 2010; Niemeijer & Collettini, 2014; Carpenter et al., 2016) that can and has been applied to the Groningen case (van den Bogert, 2015; Lele et al., 2016). However, much of the currently available laboratory data has been obtained on rock materials with a different composition, under room temperature and other non in situ conditions, and neglecting chemical effects of pore fluid or effects of loading rate and time. Moreover, little physical or mechanistic basis exists to support extrapolation of laboratory data to the in situ conditions, and to the temporal and spatial scales relevant for the Groningen field. Addressing the issue of induced seismic hazard in Groningen with the rigour required by society, regulator, government and operator requires a new level of rock physics understanding and quantification. A physics-based understanding is needed so that results can reliably be applied to interpret field observations on subsidence and seismicity and improve forward-modelling capability. The required laboratory research must address true in situ conditions and must reproduce and characterise the processes that actually occur in the field.

Since 2015, a team of senior scientists, postdoctoral researchers and PhD students at Utrecht University (UU) has been engaged in a 5-year laboratory-centred research programme sponsored by NAM. The team is led by the Experimental Rock Deformation Group (High Pressure and Temperature or HPT Laboratory) which specialises in rock and fault mechanical behaviour under the in situ temperature, pressure and chemical conditions pertaining throughout the Earth's crust. The team also incorporates specialists in electron microscopy, microstructural analysis, geophysics, hydrogeology and theoretical and numerical (geo)mechanics. The project is aimed at providing the fundamental physical understanding and quantitative description of rock and fault behaviour needed to advance understanding of reservoir compaction and fault behaviour, and the relationship between compaction and fault slip, in the context of induced seismicity and subsidence in the Groningen field. In this paper, the key knowledge gaps that drive the programme are discussed and the approaches employed to address these gaps are highlighted. First results emerging from the work done by the various PhD students and postdoctoral researchers involved are indicated briefly. Detailed results are not presented as these are under consideration by specialist journals. Rather, we focus on describing the programme and its intended contribution to understanding the response of the Groningen field to gas production. Project results will consist of multiscale models for the constitutive behaviour of sandstones (deformation/compaction) and faults (friction), which will provide input for geomechanical computation of reservoir–fault system behaviour, i.e. the coupling between reservoir deformation and fault slip, such as that conducted by Orlic & Wassing (2012), Bourne et al. (2014) and Urpi et al. (2016).

Reservoir compaction – the problem

Production of gas from reservoir rocks decreases the gas pressure (P _f) in the pores. Since the overburden stress (σ _v) is constant, gas production causes an increase in compressive vertical effective stress, σ _e=(σ _v–P _f), which generally leads to compaction. The associated changes in the stress–strain and elastic energy fields that occur within the reservoir formation, and in the surrounding rock mass, lead to surface subsidence and provide the driver for induced fracturing and fault slip, and hence induced seismicity. Understanding the compaction behaviour of the Slochteren reservoir rock, under in situ conditions, is therefore essential for understanding induced seismicity in the Groningen field and for modelling field behaviour. Other effects that can potentially influence deformation of the subsurface system, and hence subsidence and seismicity, include drainage of water and gas from low-permeability formations, such as the Ten Boer Claystone, into the depleted reservoir (Mossop, 2012), and creep effects in the rock-salt caprock driven by reservoir deformation (Marketos et al., 2016).

From a rock physics perspective, compaction of a sandstone reservoir rock occurs because decreasing fluid pressure during depletion leads to an increase in effective stresses transmitted through, and across the contacts of, the mineral grains that constitute the rock (Fig. 1B). In addition to conventional poro-elastic deformation of the grain framework, which is reversible, time-independent and relatively easily quantified (Wang, 2000), the increasing effective stresses may also cause permanent, generally time- or loading-rate-dependent compaction (Heap et al., 2009; Brantut et al., 2014; Brzesowsky et al., 2014b). Permanent deformation of this type is believed to be predominantly controlled by grain-scale yield, grain and grain contact fracturing by equilibrium or by subcritical crack growth (Atkinson, 1984), cement fracturing, micro-dissolution at contacts and/or grain rearrangement processes (Fig. 2), many of which are thermally activated (i.e. governed by stress-dependent, Arrhenius-type kinetic behaviour).

Fig. 2. Illustrations of grain-scale mechanisms that can potentially accommodate permanent (inelastic) deformation during compaction of reservoir sandstones. These operate in addition to conventional poro-elastic deformation, which is reversible. Permanent compaction and creep of sandstone may be controlled by (A) cement fracture, cement dissolution or clay mineral dewatering/desorption within grain contacts, (B) dissolution at/within quartz grain contacts, (C) grain- or grain-contact fracture or crushing, and (d) intergranular frictional slip (illustrated here by nano-slickensides observed in calcite by Verberne et al. (2014)).

As noted by Brzesowsky et al. (2011), a theoretically and experimentally based understanding of these mechanisms is needed in order to develop micromechanical models that can eventually be used in macroscale modelling of the compaction of clastic hydrocarbon reservoirs in the (post-)production phase. Laboratory measurements have been performed on samples from the Groningen reservoir, focusing on empirical quantification of compaction in terms of a compaction coefficient (linearised compaction strain per unity change in effective vertical stress) in uniaxial compaction mode (Schutjens et al., 1995; Hettema et al., 2000; Hol et al., 2015), but with limited attention for partitioning of deformation between poro-elastic, inelastic and time-dependent components, or for active deformation mechanisms such as subcritical cracking, dissolution or intergranular slip. Moreover, despite numerous more general studies of compaction of loose or poorly consolidated sands (Schutjens, 1991; Chuhan et al., 2002, 2003; Karner et al., 2005; Hangx et al., 2010; Brzesowsky et al., 2014b) and of deviatoric creep of dense sandstones (Baud & Meredith, 1997; Ngwenya et al., 2001; Heap et al., 2009; Brantut et al., 2013, 2014), relatively few data or microphysically based models exist that address the role of these processes in controlling compaction of reservoir sandstones with intermediate porosities of 15–25%, such as the Slochteren sandstone.

Detailed data that are available on permanent deformation of sands and sandstones generally assume that poro-elastic behaviour is limited by time or rate-independent inelastic yield that can be described for a rock of given porosity using a single yield envelope, consisting of shear failure and compaction cap components, drawn in deviatoric versus mean stress space (Wong & Baud, 2012; see also Fig. 3A). In reality, however, such envelopes are expected to show an expanding compaction cap (Karner et al., 2005), with ongoing yield-related damage, and to exhibit a dependence on loading rate, chemical environment and temperature, due to the role of thermally activated processes such as stress corrosion cracking (Atkinson, 1984; Brantut et al., 2014). Understanding the form and evolution of such envelopes, the associated stress–strain behaviour, and the dependence of these on porosity, deformation rate and in situ conditions, expressed through appropriate constitutive equations, is crucial for evaluating reservoir rock response to gas production. For example, exceeding the vertical effective stress at which significant inelastic compaction strain begins to occur in a porous sandstone, under the more or less uniaxial compaction conditions expected in the centre of a large, unfaulted reservoir compartment, can lead to a marked acceleration of reservoir deformation at a given rate of pressure depletion (Fig. 3B). In addition, the partitioning of reservoir rock deformation between poro-elastic versus permanent inelastic strain determines the amount of elastic energy that is stored in the reservoir and is potentially available to be released seismically, versus that dissipated as heat by irreversible ‘plastic’ yielding or creep. Conventional compaction coefficients that express the amount of compaction resulting from depletion or loading under uniaxial strain conditions are inadequate to quantify deformation behaviour occurring under more general boundary conditions. More general relationships are needed to describe reservoir deformation behaviour, in particular in the neighbourhood of faults, where steep gradients are expected in the stress–strain field and may drive seismogenic fault slip. Obtaining a detailed, mechanism-based constitutive description of the 3D mechanical response of reservoir sandstones that can be applied to any boundary conditions is therefore key to addressing and modelling how reservoir compaction is linked to fault slip and induced seismicity. Moreover, the equations obtained must embody the behaviour characterising the processes going on during depletion of the real reservoir, as well as those that might continue after production has stopped.

Fig. 3. Schematic illustrations of inelastic deformation behaviour of porous sandstones. (A) Illustration of a typical yield envelope drawn for a reservoir sandstone, consisting of shear failure and compaction cap components, drawn in deviatoric stress (Q) versus effective mean stress (P) space. The inner compaction cap (solid red curve) illustrates behaviour expected for a weak, highly porous sandstone. The outer cap (dashed red curve) illustrates behaviour expected for a denser, stronger sandstone. Note that the data delineating such envelopes are sometimes described by fitting a vertically orientated ellipse (e.g. Wong & Baud, 2012). (B) Illustration of the change in vertical effective stress (σ_e) versus axial strain (ε) curve exhibited by a porous sandstone sample showing inelastic yield (permanent deformation) during uniaxial compaction. Note that for a given rate of increase of differential (or effective vertical) stress due to steady pore pressure depletion (green arrows), loading beyond yield leads to an increase in the rate of axial deformation, compared to the same increase in stress prior to yielding.

New directions in compaction research

The approach adopted to advancing understanding of the compaction behaviour of the Slochteren reservoir rock, in the research being conducted under NAM sponsorship at UU, involves the following key aims:

1. I. To identify and quantify the inelastic micromechanical processes that operate in the Slochteren sandstone during deformation under in situ conditions (i.e. during depletion);

2. II. To determine the relative importance of elastically stored versus inelastically dissipated energy;

3. III. To develop microphysically based constitutive laws that describe compaction by the observed inelastic deformation mechanisms, including rate-dependent or creep effects, alongside poro-elastic response;

4. IV. To upscale the lab-derived compaction/dissipation models obtained to provide input for modelling at the reservoir scale.

State-of-the-art experimental rock deformation and microstructural/microbeam methods are being used to investigate the deformation behaviour and mechanisms controlling deformation of the Slochteren sandstone (aim I). The strategy followed involves studying core material taken from a strongly depleted portion of the Groningen reservoir, namely from the Zeerijp-3A well, drilled by NAM in 2014/15 (Fig. 4), as well as undepleted core taken from the nearby Stedum-1 well in the early days of field exploration and exploitation.

Fig. 4. (A) NAM's Zeerijp-3A drilling operation, (B) sections of recovered core and (C) map showing location of the Zeerijp-3A (ZRP-3A) well and Stedum-1 (SDM-1) wells (taken from NAM (2016)). The Zeerijp-3A well was drilled in 2014/15 recovering samples from the depleted reservoir. Samples recovered from Stedum-1 represent the (near) undepleted state (photographs kindly provided by Jan van Elk, NAM).

A central component of the research is experimental work conducted using triaxial testing machines (Fig. 5A and B) to investigate the deformation behaviour of cylindrical sandstone samples (Fig. 5C) taken from both undepleted and depleted sections of the reservoir, i.e. sampled from the Stedum-1 and Zeerijp-3A wells mentioned above. The effects of variations in porosity and grain size (distribution) within the reservoir rock will be investigated by testing samples showing systematically chosen values. Samples from quarry outcrops of lateral equivalents of the Slochteren sandstone outcropping in Germany (Flechtinger sandstone; see Fischer et al., 2007) are also being used where core material is limited. Systematic sets of experiments are being performed at true in situ P–T conditions, in situ chemical/pore-fluid conditions and at a variety of loading/unloading and stress relaxation conditions, addressing loading-rate-dependence, creep and time-independent effects. The energy stored elastically versus that dissipated during compaction is addressed by direct measurement of elastic strain, hence strain energy, reversibly stored in experimentally tested samples, alongside calculations of the total mechanical work done during compaction from overall stress–strain behaviour (aim II). Acoustic emission methods, pore fluid additives (crack-growth and dissolution inhibitors (e.g. see Hangx et al., 2010; Brzesowsky et al., 2014a)) and isotopic tracers are being employed to help identify active deformation processes at the grain scale, e.g. to assess if subcritical crack growth, reaction or solution–precipitation processes control deformation. The latest optical, electron microscopy (Focused Ion Beam Scanning Electron Microscopy (FIB-SEM), Electron Backscatter Diffraction (EBSD), Transmission Electron Microscopy (TEM)) and microbeam analysis methods (Secondary Ion Mass Spectroscopy (SIMS), Fourier Transform Infrared (FTIR) and Raman spectroscopy, Laser Ablation Microprobe Inductively Coupled Plasma Mass Spectrometry (LAM-ICP-MS)) are also being employed to identify the microscale physical processes that accommodate deformation under different loading conditions (Fig. 6A), along with cutting-edge techniques such as grain-scale displacement/deformation analysis (Quintanilla-Terminel et al., 2017) (see Fig. 6B). One of the most exciting developments being pursued is that of 4D X-ray microtomography (μCT), which allows imaging of samples during active deformation under in situ P–T conditions (Renard et al., 2016). This is now possible at modest P–T conditions using small-scale deformation cells that enable deformation of 5mm by 10mm cylindrical samples. The Utrecht team is involved in a collaborative project with Prof. F. Renard (Oslo University, Norway) and the European Synchrotron Radiation Facility (ESRF) in Grenoble, France, to image active sandstone deformation and to quantify how strain is accumulated and propagates as individual grains deform. The first ever μCT ‘see-through’ time-lapse experiment on Slochteren sandstone was conducted by the Utrecht team in October 2016, with success, at in situ reservoir conditions (Fig. 7). Similar μCT and in situ triaxial deformation experiments are also being investigated in collaboration with Ghent University, Belgium (Prof. V. Cnudde, Department of Physics).

Fig. 5. (A, B) Triaxial testing machines and (C) typical 25mm diameter sample of the Slochteren Sandstone used in the experimental programme on reservoir compaction at Utrecht University.

Fig. 6. (A) Back-scattered electron image illustrating the microstructure visible in a polished section of an experimentally deformed sandstone sample. Irreversible, inelastic deformation occurred by microcracking at the grain scale. (B) Secondary electron image of grains in the polished surface of a sandstone sample deformed in the laboratory under in situ reservoir conditions. Note the displacement markers (micro-trenches) that were ‘machined’ into the polished sample surface prior to deformation. These initially straight micro-trenches were offset during deformation demonstrating strain accommodation by intergranular sliding, alongside grain fracturing visible in the central sand grain.

Fig. 7. (A) Triaxial testing machine enabling time-lapse X-ray μCT imaging of actively deforming sandstone samples under in situ conditions at the European Synchrotron Radiation Facility in Grenoble (collaboration with ESRF, Université Grenoble Alpes and Oslo University). (B, C) Two-dimensional X-ray μCT image of a sample of Slochteren Sandstone prior to deformation (B) and during deformation (C) at in situ Groningen reservoir P–T conditions in the apparatus shown in (A).

Though the mechanisms causing permanent deformation of sandstones have been studied and imaged (see Fig. 2), quantitative descriptions of the grain-scale strain response to stress are still lacking, particularly in terms of the time dependence. The above-mentioned systematic experimental deformation work and micro-mechanistic studies will identify the main deformation mechanism responsible for compaction under different environmental and boundary/loading conditions. On this basis, work is in progress to develop mechanism-based models describing the deformation and compaction behaviour of Slochteren sandstone samples as a function of sample porosity and grain size, under specific and general loading conditions, i.e. to obtain constitutive equations relevant for reservoir conditions (aim III). Microphysical modelling methods will extend those recently used by the UU team and other groups to obtain constitutive models for deformation, compaction and creep of sands (Guéguen & Fortin, 2013; Brzesowsky et al., 2014a, b), sandstones (Brantut et al., 2014) and granular fault rocks (Niemeijer & Spiers, 2007) via a micro-mechanistic description. The models obtained will be treated as hypothetical descriptions of deformation behaviour that will be tested against experiments and calibrated accordingly. Simplified upscaling from lab to modelling mesh and reservoir scales will be achieved by volume averaging with respect to porosity and compositional variations across the full vertical extent of the reservoir, at different surface locations, making use of well- and core-log data (aim IV). This will provide descriptions of reservoir compaction behaviour that vary laterally with lithological changes in the reservoir, notably porosity. A further approach being considered for modelling and upscaling sandstone deformation behaviour is to incorporate microphysical models, developed to describe grain-scale processes such as grain crushing or dissolution, into Discrete Element Modelling (DEM) codes in the form of new particle interaction laws.

An important part of the work on reservoir compaction is the need to confirm that the deformation mechanisms observed and characterised in the laboratory experiments, in terms of constitutive descriptions, are similar to those occurring within the reservoir during gas production. This is being attempted by applying the same state-of-the-art optical, electron microscopy (FIB-SEM, EBSD, TEM) and microbeam analysis methods, employed for studying the deformation mechanisms operating in lab-deformed samples, to (a) undeformed control samples taken from the near-virgin reservoir, i.e. from the Stedum-1 well, and (b) samples taken directly from the depleted reservoir, i.e. ‘production-deformed’ core material recovered from the Zeerijp-3A well. The aim is to assess what inelastic deformation processes may have operated in situ during depletion of the reservoir, through quantitative comparison of microstructures observed in samples taken from the depleted and undepleted reservoir intervals, and whether similar processes are active in the laboratory experiments. Possible microstructural differences that are being investigated between the undepleted and depleted core samples include differences in microcrack density and orientation distribution. Work is also in progress in collaboration with Christian-Albrechts University of Kiel (Dr H. Bahadur Motra, Institute of Applied Geoscience) to assess if such differences, and possibly even the in situ stress state in depleted reservoirs, can be sensed by means of 3D acoustic velocity measurements (see also Sayers & Kachanov, 1995; Sayers, 2006; Kern, 2011) combined with microstructural work on crack orientation distribution.

The challenge of quantifying fault friction and its role in controlling induced seismicity

To date, most numerical studies investigating fault (re)activation due to subsurface activities, such as gas production or hydraulic stimulation of oil, gas or geothermal fields, evaluate the potential for fault slip using a standard static Mohr–Coulomb failure or Byerlee friction criterion (e.g. Orlic & Wassing, 2012). This approach is adequate for assessing whether faults will be activated, but not for evaluating whether fault motion will be aseismic or seismic, what the areal extent of the rupture and the displacement field associated with it will be, what changes in stress will occur, or how much strain energy will be released. However, this detail is essential if we are to evaluate the potential for compaction-driven fault slip and induced seismicity in gas fields such as Groningen. To capture such detail, constitutive relations are needed that robustly describe how fault strength depends on slip rate, slip distance and time, at a length scale commensurate with mesh scales used in geomechanical modelling codes (1–10–100m; see Orlic & Wassing, 2012).

Dynamic fault friction laws have been widely used in studies of natural earthquake rupture (Noda & Lapusta, 2013; Noda, 2016), and in a few recent studies of induced seismicity (Rutqvist et al., 2014; Urpi et al., 2016). These friction laws are referred to as Rate and State Dependent Friction (or RSF) laws (Dieterich, 1979; Ruina, 1983), and have been successful in reproducing the seismic versus aseismic slip behaviour of faults, as well as fault healing and reactivation, in a qualitative sense. However, these relations are largely empirically derived, are often based on friction experiments done at room temperature only, and have limited physical basis to allow extrapolation to in situ P–T and chemical conditions (Bos & Spiers, 2002; Niemeijer & Spiers, 2007; Carpenter et al., 2016). The key parameters in them also vary strongly from material to material, as well as with temperature, effective stress and slip rate, with limited understanding of why. Moreover, they are determined under fixed displacement rate boundary conditions, by introducing step changes in displacement rate, and their validity for other boundary conditions, for length scales beyond the lab scale, and for slip rates too slow or too fast to study easily in the lab is unknown. Finally, static fault rock failure plus displacement weakening effects are only partly accounted for through the incorporation of fault strength recovery or healing effects that evolve with static fault-healing duration (Dieterich, 1972; Karner et al., 1997; Nakatani & Scholz, 2004).

As a result, though very useful for investigating qualitative trends, standard RSF laws have not been fully successful in quantitatively modelling earthquake fault behaviour, except by tuning the embedded (poorly constrained) parameters, in particular the length scale or critical slip distance associated with re-establishing steady-state sliding friction after a step change in sliding velocity. An additional important point is that most RSF descriptions to date have been obtained from experiments addressing low slip velocities (0.01–100μms⁻¹), which are the most relevant to rupture nucleation. If slip velocities accelerate unstably towards higher seismic velocities, other dynamic processes may be activated, such as frictional heating and associated pore fluid pressurisation (Mase & Smith, 1987; Andrews, 2002), which can drastically reduce effective normal stress and hence fault friction. These high-velocity friction effects become important at slip velocities in excess of c. 1–10 cms⁻¹, depending on normal stress, and have been extensively investigated in high-speed (1–3ms⁻¹), large-displacement (1–10m) sliding experiments, in the context of natural earthquake fault behaviour (Mizoguchi et al., 2007; Brantut et al., 2011; De Paola et al., 2011; Di Toro et al., 2011; Niemeijer et al., 2012). Whether such effects play a role in induced seismic slip, which is generally limited in displacement, is still poorly known, though evidence has been reported for significant frictional heating in laboratory stick-slip events (Brantut et al., 2016).

Clearly, then, to model induced seismicity in the Groningen field, better mechanistically based, RSF-type constitutive models that cover the entire velocity range for fault slip, and that incorporate fault healing and slip weakening effects, are needed (e.g. Chen & Spiers, 2016). These must apply to fault rocks (crushed rock wear products or fault ‘gouges’), that are representative for the Groningen field, being derived for example from the key lithologies present in the field, or from mixtures of these materials simulating gouge smearing effects. Little or no such friction data has been previously available for Groningen rock compositions. Also essential is a first attempt to develop upscaling rules, linking fault friction measured at lab length scales (1–10cm) to the geomechanical modelling mesh-scale (1–10m, above which geometric irregularities can be explicitly accounted for in the model mesh). Upscaling models for fault friction have been developed assuming multiscale asperity or roughness distributions (e.g. Kemeny & Hagaman, 1992). Moreover, recent studies on natural fault surfaces have shown that some degree of self-similar scaling indeed exists between roughness profiles at scales ranging from less than 0.1mm to more than 1000km in length (Candela et al., 2012; Renard et al., 2013). However, there are currently no multiscale friction data covering the required lab- to mesh-scale range for gouge-filled faults, and there are no experimentally validated theoretical or numerical models for predicting their friction behaviour at the 1–10m scale. Experiments at different scales are therefore needed to establish scaling relationships for key fault friction parameters across scales ranging from mm/cm to dm/m. In addition, the possible effects of frictional heating and fluid pressurisation phenomena need to be evaluated.

Addressing the fault friction issue in the context of Groningen

The work on fault reactivation and frictional behaviour, in progress at UU in the context of the NAM-funded project on Groningen, is directed at understanding the behaviour of pre-existing faults that cut the reservoir and over/underburden system. The broad aims are:

1. 1) To obtain constitutive laws describing the failure, slip weakening and rate- and state-dependent frictional behaviour of simulated fault rocks representing the reservoir and over/underburden rock compositions present in the Groningen field. For extrapolation beyond laboratory conditions, these must address the physical mechanisms that operate in faults in situ;

2. 2) To investigate upscaling relations enabling fault mechanical behaviour and friction measurements obtained in laboratory experiments to be upscaled across the 1cm–10cm–1m–10m scales, i.e. at least to the scale of finite element meshes used in rupture modelling;

3. 3) To incorporate the fault mechanics models obtained into (quasi-)dynamic rupture simulators;

4. 4) To test resulting models and rupture simulations against multiscale laboratory experiments;

5. 5) To apply the fault mechanics models developed to simulate fault reactivation, induced seismic events and earthquake statistics.

The experimental work on fault strength and rate- and state-dependent frictional behaviour (aim 1) is being conducted using a range of direct shear testing assemblies that allow simulated fault gouge layers c. 1mm in thickness and roughly 40mm by 35mm in area to be sheared in conventional triaxial testing equipment under true in situ P–T and pore fluid chemical conditions. Imposed slip rates from 10nms⁻¹ to 1mms⁻¹ can be investigated, and work is in progress to allow strength evolution to be monitored during unstable slip events at seismic slip velocities, as well as slip rate evolution under conditions of controlled shear stress. The controlling microphysical mechanisms are investigated using the same optical and electron optical methods employed to study deformation processes in the reservoir sandstone. The experimentally measured frictional behaviour is cast first into conventional RSF and healing rate descriptions (Dieterich, 1979; Ruina, 1983; Marone, 1998; Scholz, 1998) to obtain an immediately usable description for geomechanical and rupture modelling work in progress at Shell Global Solutions, Rijswijk, and the Netherlands Organisation for Applied Scientific Research (TNO), Utrecht. In addition, mechanism-based friction, slip weakening and healing models, that can be rigorously tested against the experimental results, are being developed following the microphysical theory developed by Niemeijer & Spiers (2007), den Hartog & Spiers (2014) and Chen & Spiers (2016). In contrast to standard RSF models, these theoretical models are based on competition between granular flow and thermally activated, grain-scale interaction and deformation processes relevant for fault gouge deformation as opposed to rock-on-rock sliding friction. They contain parameters with direct physical meaning, such as porosity, grain size and fault rock dilatation angle, while the parameters describing the rate dependence of friction can be directly related to kinetic parameters for processes such as microscale plasticity, crack growth or stress-driven dissolution. These new models (e.g. Chen & Spiers, 2016) provide a framework for extrapolating laboratory data to subsurface conditions under any boundary conditions and are able to reproduce all of the transient RSF and healing behaviours typically seen in gouge friction experiments.

In the context of upscaling, field studies of natural faults that cut lateral equivalents of (and analogues for) the Slochteren sandstone are being conducted, where exposed at the surface in the UK and Germany. These investigations are confirming that natural (tectonic) shear on faults, with mechanically significant throw (10–100m), is accommodated in gouge-filled fault cores and is often highly localised internally into mm- or cm-wide principal slip zones (Fig. 8). This suggests that induced slip events will also occur within similar principal slip zones and lends confidence about the field-scale relevance of laboratory friction experiments performed on simulated gouge samples of just 1–2mm in thickness. However, such field studies also point to the ubiquitous development of a multiscale, self-similar, internal hierarchy of cross-cutting or anastomosing Riedel shears, P-shears and/or Y-shears (Fig. 8), remarkably similar to the structures observed in experimentally sheared samples of simulated gouges (Logan et al., 1979), and suggesting that shear bands may control fault zone strength. Motivated by this observation, Finite Element Methods (FEM) are being applied to explore the mechanical significance of such multiscale networks of localised, weak shear band structures, and to establish whether upscaling laws can be derived that link friction measured in direct shear tests at the cm length scale to the 10cm, 1m, 10m and 100m scales (aim 2). Classical asperity length-scaling methods are also being considered (Kemeny & Hagaman, 1992). The existence of scaling relations and the validity of the above scaling models is being simultaneously investigated by comparison with experiments performed on the 5–50cm length scale, systematically varying the simulated fault rock structure, its compositional components, fault geometry and asperity structure (aim 4). The smaller-scale experimental data (5–25cm) emerge directly from direct and ring shear experiments in progress at UU. Large-sample (25–50cm), laboratory-scale fault reactivation experiments, with acoustic emission, displacement/strain field and rupture propagation monitoring capability, are being conducted using the large biaxial testing machine at the Laboratory for Earthquake Dynamics, China Earthquake Administration (CEA), Beijing (see Jiang et al., 2003; Ma et al., 2010), employing rock and polymer forcing blocks to vary system stiffness.

Fig. 8. (A) Field exposure of a fault juxtaposing Rotliegend-equivalent sandstone against Upper Carboniferous shales (Whitley Bay, UK). Note internal structure of localised principal slip zones (PSZ) and shear bands visible at the 1–10cm scale. (B) Internal P, R1 and Y shear band structures at the 1–10mm scale.

In the context of aim 3, the fault friction laws emerging from our experiments, and ultimately from our upscaling work, will be implemented into quasi-dynamic and dynamic rupture models, such as the open-source, seismic cycle simulator for modelling natural fault slip/rupture, QDYN (https://github.com/ydluo/qdyn), but also fully dynamic finite element codes. This work is being performed in close collaboration with Prof. J.-P. Ampuero of CalTech, who developed the QDYN simulator, and with the geomechanics team at TNO using the FEM package DIANA. The resulting rupture simulations will in turn be tested by comparison with small and larger-scale rupture experiments performed at Utrecht and using the large biaxial machine at CEA (also aim 4). Under aim 5, the multiscale fault strength models obtained from experiments and upscaling relations will be applied in fully dynamic, reservoir-scale geomechanical models to evaluate (a) the conditions under which seismogenic slip can be initiated on faults at the reservoir scale, (b) the rupture area, slip magnitude, stress drop and energy release (potential earthquake magnitude) associated with seismic slip events, and (c) event frequency versus magnitude relationships. Again this will be done in collaboration with TNO using the DIANA FEM platform, with NAM/Shell Global Solutions using in-house software, and with Prof. Ampuero at CalTech.

More on upscaling: ‘big friction experiments’ at NIED (Japan)

As indicated already, in order to model fault rupture and associated seismic or aseismic fault slip, in a finite element or similar modelling framework addressing dynamic fault rupture at the reservoir scale, fault friction must be described at least at the mesh scale, i.e. at a scale of at least 0.5–2m in current codes. In order to extend the length scales that can be accessed for establishing fault friction scaling relations beyond the 5–50cm range, and for testing scaling models derived from FEM analysis and the other methods discussed above, collaboration has recently been established between the UU team and the group responsible for the Large Scale Earthquake Simulator facility at the Japanese National Research Institute for Earth Science and Disaster Resilience (NIED), in Tsukuba, Japan (Prof. E. Fukuyama, Dr F. Yamashita). In 2018, the two groups will collaborate, via NAM funding, to conduct the first-ever experiments on simulated gouge-bearing faults at the 1.0 and 1.5m fault length scale. The NIED group recently published results of rock-on-rock friction experiments at this scale (Togo et al., 2015; Yamashita et al., 2015), which show that the threshold velocity for dynamic weakening between (initially) bare sliding surfaces of gabbro is lower than that determined using cm-scale samples.

The facility in question (Fig. 9) consists of a so called ‘shaking table’ used to drive one large rectangular block of rock over another. The blocks are loaded normal to the shearing direction by roller-mounted flat-jacks, capable of supplying a normal load of up to 6 MPa, depending on the fault area. The shaking table is capable of driving the bottom block at shearing velocities ranging between 10⁻⁴ and 0.1ms⁻¹, but can only be run for a total duration of 15min. Forward shear motion is possible, but not reverse. Load cells are used to measure the macroscopic shear and normal loads. Displacement can be measured in three directions using three laser displacement transducers, attached to the edge of the blocks. Strain gauges can be installed on the blocks to measure local strains, which can be used to estimate local stress. We also plan to install one or more high-speed video cameras to track the evolution of the displacement and velocity field occurring throughout the block–gouge assembly, using digital image correlation to track surface marker displacements. Accelerometers will also be deployed at appropriate points on the sliding block assembly.

Fig. 9. Large Scale Earthquake Simulator facility at the Japanese National Research Institute for Earth Science and Disaster Resilience (NIED), in Tsukuba, Japan.

We will use relatively low-porosity (~10%) quartz-rich sandstone blocks of 0.5, 1.0 and 1.5m in length scale, sandwiching a simulated Slochteren sandstone fault gouge of varying thickness. The preliminary plan is to shear the gouge layer at 100μms⁻¹ to determine its macroscopic friction for direct comparison with conventional 5cm-scale friction experiments and 10–25cm ring shear friction experiments at UU. In these 5–25cm-scale experiments, we will reproduce the lower NIED machine stiffness (1.8GPam⁻¹ vs ~100GPam⁻¹) and the use of sandstone forcing blocks at low normal stress, to optimise cross-comparison of results between the two scales. The scaling interval between 25 and 50cm is covered through our experimental programme at the China Earthquake Administration. The internal fault gouge microstructure will be investigated at all scales. The results will be used to assess if empirical friction scaling relations can be established up to the 1.5m scale, for testing our FEM-based and other upscaling models, and for testing the rupture modelling approaches described above. The experiments planned at NIED offer the first opportunity to date to validate rupture models of the type needed to simulate induced seismicity in the Groningen field.

First results

The first results to emerge from the programme are presently being submitted to specialist journals so will not be reported here beyond a brief indication of some of the highlights. A key point to emerge from the compaction experiments on the reservoir sandstones under simulated in situ conditions is that up to 40 or 50% of the deformation occurring in samples with porosities of 19–20% and above is inelastic and irreversible (see Fig. 3B for a qualitative illustration of increased permanent or inelastic deformation beyond yield). This is accommodated by a (coupled) combination of intergranular sliding and grain fracturing. In such highly porous sandstones, the associated mechanical behaviour seems to be governed by a porosity-dependent yield envelope, featuring both shear failure and compaction cap portions (Fig. 3A). With ongoing compaction, or decreasing porosity, the compaction cap expands towards higher mean stresses (P), resulting in a continuously growing yield envelope, as schematically illustrated in Figure 3A (see also Karner et al., 2005). If such behaviour occurs in situ during reservoir depletion, then stress and strain development within the reservoir will differ markedly from expectations based on conventional poro-elastic compaction predictions. This requires investigation via future geomechanical modelling. Time-dependent inelastic deformation (creep) is also observed in the experiments and is clearly affected by pore fluid composition, implying that stress corrosion cracking processes may be involved. The rates are relatively rapid, though, producing mainly transient strains, and can be viewed as contributing to time-independent reservoir deformation on field depletion timescales. Microstructural studies of reservoir samples retrieved from the Zeerijp-3A well (depleted state) and the Stedum-1 well (pre-/early depletion) show only subtle differences. Quantitative comparison of crack densities in samples that contain sedimentary bands of similar porosity and grain size indicates higher crack densities in the Zeerijp-3A material than in the Stedum-1 core, suggestive of inelastic deformation by the same grain-scale microcracking mechanisms as seen in our experiments. More microstructural work is needed to confirm or refute this, along with results from the planned acoustic velocity experiments.

Work on fault friction has progressed rapidly, yielding quantitative descriptions of the rate- and state-dependent frictional properties of all the main rock units in the Groningen reservoir system, and of mixtures simulating fault gouge smearing effects. This has enabled a mechanical stratigraphy to be established for the strength, healing and rate- and state-dependent frictional properties of faults cutting the reservoir and over/underburden intervals. Major effects of reservoir temperature, pore fluid chemistry and gas/water saturation have been observed and quantified, specifically in simulated fault gouges derived from the basal Zechstein carbonate/anhydrite unit that overlies the Ten Boer Claystone and reservoir. Field studies of faults in the Rotliegend succession and laterally equivalent analogues exposed in the UK support the hypothesis that internal shear band structure plays a role in controlling the strength and frictional behaviour of faults at multiple scales. In addition, first efforts in numerical modelling of the effect of internal shear band development on fault strength and frictional behaviour confirm scale-dependent effects. Friction experiments performed at the 5 and 50cm length scale in Utrecht and at CEA have demonstrated complex differences in behaviour, probably related to differences in experimental details. Substantial progress has been made in explaining these differences, however, as well as in FE modelling of the stress–strain field and fault slip behaviour seen in the 50cm-scale experiments conducted at CEA. Geomechanical modelling of fault reactivation and rupture propagation at the reservoir scale has progressed assuming simple, slip-weakening frictional behaviour to date. The results show that the frictional stratigraphy of the units cut by an individual fault (cf. Fig. 1A) has a strong influence on the upward and downward extent of rupture propagation and hence on seismogenic potential, making inclusion of the advanced fault friction laws emerging from the programme a must.

Conclusions

A NAM-funded research programme has been initiated at Utrecht University with the aim of providing a fundamental, physically based understanding and quantitative description of reservoir deformation and fault mechanical behaviour in the Groningen gas field. The programme will run for the period 2015–19. It involves an integrated approach employing experimental rock and fault mechanics work conducted at in situ conditions, microscale observational studies to determine the physical processes that control reservoir rock deformation and fault slip in the field and in the laboratory, plus modelling and experimental work aimed at establishing upscaling rules between laboratory and field scales. The results obtained are providing constitutive laws describing reservoir deformation and fault friction behaviour. These are being included in geomechanical models aimed at coupling reservoir compaction and fault reactivation to investigate fault rupturing and seismogenic behaviour at the reservoir scale.

The work is being executed by a team of senior scientists, postdoctoral researchers and PhD students with expertise ranging from experimental rock and fault mechanics to electron microscopy, geophysics and numerical modelling, and involves important collaboration with earthquake science groups in the USA, Japan and China. The key basic science questions addressed by the programme are: what processes control deformation of the Slochteren Sandstone under in situ reservoir conditions? how can we describe these processes and the associated mechanical behaviour quantitatively? what processes control fault mechanical behaviour in the Groningen situation? how should fault reactivation and frictional slip behaviour be described as input for geomechanical and earthquake rupture modelling? how do we upscale lab data to the field-scale? and how do the processes that are found to control reservoir and fault behaviour impact induced earthquake generation? First results are providing important new insights, notably into (a) the contribution of inelastic deformation to reservoir deformation behaviour, (b) fault friction and healing behaviour in simulated fault rocks derived from the rock types present in, above and below the Slochteren Sandstone reservoir unit, and (c) the crucial role played by these aspects of reservoir and fault behaviour in determining the response of the Groningen field to gas production.