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, Reference Mossop2012), and creep effects in the rock-salt caprock driven by reservoir deformation (Marketos et al., Reference Marketos, Spiers and Govers2016).

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, Reference Wang2000), the increasing effective stresses may also cause permanent, generally time- or loading-rate-dependent compaction (Heap et al., Reference Heap, Baud, Meredith, Bell and Main2009; Brantut et al., Reference Brantut, Heap, Baud and Meredith2014; Brzesowsky et al., Reference Brzesowsky, Spiers, Peach and Hangx2014b). 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, Reference Atkinson1984), 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. (Reference Verberne, Plümper, Matthijs de Winter and Spiers2014)).

As noted by Brzesowsky et al. (Reference Brzesowsky, Spiers, Peach and Hangx2011), 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., Reference Schutjens, Fens, Smits, Barends, Brouwer and Schröder1995; Hettema et al., Reference Hettema, Schutjens, Verboom and Gussinklo2000; Hol et al., Reference Hol, Mossop, van der Linden, Zuiderwijk and Makurat2015), 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, Reference Schutjens1991; Chuhan et al., Reference Chuhan, Kjeldstad, Bjørlykke and Høeg2002, Reference Chuhan, Kjeldstad, Bjørlykke and Høeg2003; Karner et al., Reference Karner, Chester, Chester, Kronenberg and Hajash2005; Hangx et al., Reference Hangx, Spiers and Peach2010; Brzesowsky et al., Reference Brzesowsky, Spiers, Peach and Hangx2014b) and of deviatoric creep of dense sandstones (Baud & Meredith, Reference Baud and Meredith1997; Ngwenya et al., Reference Ngwenya, Main, Elphick, Crawford and Smart2001; Heap et al., Reference Heap, Baud, Meredith, Bell and Main2009; Brantut et al., Reference Brantut, Heap and Meredith2013, Reference Brantut, Heap, Baud and Meredith2014), 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, Reference Wong and Baud2012; see also Fig. 3A). In reality, however, such envelopes are expected to show an expanding compaction cap (Karner et al., Reference Karner, Chester, Chester, Kronenberg and Hajash2005), 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, Reference Atkinson1984; Brantut et al., Reference Brantut, Heap, Baud and Meredith2014). 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, Reference Wong and Baud2012). (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., Reference Fischer, Gaupp, Dimke and Sill2007) 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., Reference Hangx, Spiers and Peach2010; Brzesowsky et al., Reference Brzesowsky, Hangx, Brantut and Spiers2014a)) 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., Reference Quintanilla-Terminel, Zimmerman, Evans and Kohlstedt2017) (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., Reference Renard, Cordonnier, Dysthe, Boller, Tafforeau and Rack2016). 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, Reference Guéguen and Fortin2013; Brzesowsky et al., Reference Brzesowsky, Hangx, Brantut and Spiers2014a, b), sandstones (Brantut et al., Reference Brantut, Heap, Baud and Meredith2014) and granular fault rocks (Niemeijer & Spiers, Reference Niemeijer and Spiers2007) 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, Reference Sayers and Kachanov1995; Sayers, Reference Sayers2006; Kern, Reference Kern2011) 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, Reference Orlic and Wassing2012). 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, Reference Orlic and Wassing2012).

Dynamic fault friction laws have been widely used in studies of natural earthquake rupture (Noda & Lapusta, Reference Noda and Lapusta2013; Noda, Reference Noda2016), and in a few recent studies of induced seismicity (Rutqvist et al., Reference Rutqvist, Cappa, Rinaldi and Godano2014; Urpi et al., Reference Urpi, Rinaldi, Rutqvist, Cappa and Spiers2016). These friction laws are referred to as Rate and State Dependent Friction (or RSF) laws (Dieterich, Reference Dieterich1979; Ruina, Reference Ruina1983), 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, Reference Bos and Spiers2002; Niemeijer & Spiers, Reference Niemeijer and Spiers2007; Carpenter et al., Reference Carpenter, Ikari and Marone2016). 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, Reference Dieterich1972; Karner et al., Reference Karner, Marone and Evans1997; Nakatani & Scholz, Reference Nakatani and Scholz2004).

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, Reference Mase and Smith1987; Andrews, Reference Andrews2002), 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., Reference Mizoguchi, Hirose, Shimamoto and Fukuyama2007; Brantut et al., Reference Brantut, Han, Shimamoto, Findling and Schubnel2011; De Paola et al., Reference De Paola, Hirose, Mitchell, Di Toro, Viti and Shimamoto2011; Di Toro et al., Reference Di Toro, Han, Hirose, De Paola, Nielsen, Mizoguchi, Ferri, Cocco and Shimamoto2011; Niemeijer et al., Reference Niemeijer, Di Toro, Griffith, Bistacchi, Smith and Nielsen2012). 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., Reference Brantut, Passelègue, Deldicque, Rouzaud and Schubnel2016).

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, Reference Chen and Spiers2016). 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, Reference Kemeny and Hagaman1992). 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., Reference Candela, Renard, Klinger, Mair, Schmittbuhl and Brodsky2012; Renard et al., Reference Renard, Candela and Bouchaud2013). 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, Reference Dieterich1979; Ruina, Reference Ruina1983; Marone, Reference Marone1998; Scholz, Reference Scholz1998) 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 (Reference Niemeijer and Spiers2007), den Hartog & Spiers (Reference den Hartog and Spiers2014) and Chen & Spiers (Reference Chen and Spiers2016). 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, Reference Chen and Spiers2016) 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., Reference Logan, Friedman, Higgs, Dengo and Shimamoto1979), 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, Reference Kemeny and Hagaman1992). 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., Reference Jiang, Ma, Zhang, Cao and Hou2003; Ma et al., Reference Ma, Ma, Liu and Liu2010), 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., Reference Togo, Shimamoto, Yamashita, Fukuyama, Mizoguchi and Urata2015; Yamashita et al., Reference Yamashita, Fukuyama, Mizoguchi, Takizawa, Xu and Kawakata2015), 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.