# Reservoir creep and induced seismicity: inferences from geomechanical modeling of gas depletion in the Groningen field

Reservoir creep and induced seismicity: inferences from geomechanical modeling of gas depletion... Summary The Groningen gas field in the Netherlands experienced an immediate reduction in seismic events in the year following a massive cut in production. This reduction is inconsistent with existing models of seismicity predictions adopting compaction strains as proxy, since reservoir creep would then result in a more gradual reduction of seismic events after a production stop. We argue that the discontinuity in seismic response relates to a physical discontinuity in stress loading rate on faults upon the arrest of pressure change. The stresses originate from a combination of the direct poroelastic effect through the pressure changes and the delayed effect of ongoing compaction after cessation of reservoir production. Both mechanisms need to be taken into account. To this end, we employed finite-element models in a workflow that couples Kelvin-Chain reservoir creep with a semi-analytical approach for the solution of slip and seismic moment from the predicted stress change. For ratios of final creep and elastic compaction up to 5, the model predicts that the cumulative seismic moment evolution after a production stop is subject to a very moderate increase, 2–10 times less than the values predicted by the alternative approaches using reservoir compaction strain as proxy. This is in agreement with the low seismicity in the central area of the Groningen field immediately after reduction in production. The geomechanical model findings support scope for mitigating induced seismicity through adjusting rates of pressure change by cutting down production. Geomechanics, Induced seismicity 1 INTRODUCTION Time-dependent compaction of hydrocarbon reservoirs is a well-known phenomenon. Laboratory experiments indicate that time-dependent compaction can be of the same order of magnitude as direct (elastic) compaction for clastic and shale reservoirs (De Waal 1986; Hettema et al.2000; Sone & Zoback 2014; Pruiksma et al.2015). Results from subsidence inversion demonstrate that reservoir compaction rates are delayed relative to pressure drop, in line with laboratory data (Hettema et al.2000; Mossop 2012; Fokker & Van Thienen-Visser 2016). Compaction of clastic reservoir rocks is predicted to continue for decades after a potential shut-in of production or end of field life, and creep may constitute a significant portion (up to 10 per cent) of the compaction prior to cessation of production (Mossop 2012; Pruiksma et al.2015). In addition, time-dependent strain of overburden and underburden rock can contribute to time-dependent compaction and subsidence (Orlic & Wassing 2013; Chang et al.2014; Marketos et al.2015). Bourne et al. (2014) introduced a predictive model for seismicity due to fluid extraction from hydrocarbon reservoirs, using compaction strains as a proxy for potential stress change and seismic moment, in analogy to the strain–stress seismic moment proxy introduced by Kostrov (1974) and used for induced seismicity by McGarr (2014). For time-dependent compaction following a production stop, this model implies that seismicity rates remain significant and reduces slowly. This outcome was unexpected to us, as it does not agree with the most recent observations in the Groningen field in the Netherlands: a significant cut in production rates in the central area that was applied in 2014 was immediately followed by a statistically significant reduction in seismic event rates in the next year (Nepveu et al. 2016). In this paper, we demonstrate that the unexpected outcome of the proxy model can be attributed to a break which arises in the compaction strain–stress relationship, and which is related to the cut in production. The fundamental difference between our study and that of Bourne et al. (2014) is that we consider the stresses acting on the activated faults as the driving parameters of induced seismicity, rather than using the reservoir compaction as a proxy (Fig. 1). These stresses originate from a combination of the direct poroelastic effect through the pressure changes and from the delayed effect of ongoing compaction after cessation of reservoir production. Both mechanisms need to be taken into account, and this can explain the observed break. To demonstrate this, a simplified geomechanical model for the central area in the Groningen field is employed. Section 2 introduces the relevant data of the Groningen gas field and their significance in unraveling the relative control of pressure change and time-dependent compaction on seismicity. In Section 3, the geomechanical modeling approach is introduced. It is based on a viscoelastic (Kelvin Chain) finite-element solution for a two-dimensional reservoir section. The model evaluates the stress response of the reservoir, including time-dependent reservoir compaction, and highlights the reduction of stress build-up upon termination of production. In Section 4, the results in terms of seismic moment increase are evaluated using a novel analytical approach to predict seismic moment from the viscoelastic stress solution. The analytical technique is required as the finite-element code we use does not allow for combining Kelvin-Chain creep and Mohr–Coulomb failure. It is demonstrated that the model predicts a significant reduction of seismic moment increase after production stop, qualitatively in agreement with the observed reduction of event rates. We argue that the findings of this study are expected to hold for depleting reservoirs in general, assuming similar reservoir conditions. Figure 1. View largeDownload slide Schematic diagram of the geomechanical workflow in this paper, marked by blue lines. The model predicts a physical discontinuity in stress loading rate on faults upon the arrest of pressure change after cessation of reservoir production, with some delayed effect on stress change due to creep. The orange line represents the alternative model using compaction strain as proxy to predict seismic moment (Bourne et al.2014). The numbers refer to sections in the paper. Figure 1. View largeDownload slide Schematic diagram of the geomechanical workflow in this paper, marked by blue lines. The model predicts a physical discontinuity in stress loading rate on faults upon the arrest of pressure change after cessation of reservoir production, with some delayed effect on stress change due to creep. The orange line represents the alternative model using compaction strain as proxy to predict seismic moment (Bourne et al.2014). The numbers refer to sections in the paper. 2 GRONINGEN FIELD DATA The Groningen gas field is a very large gas field located onshore in the Netherlands. Gas production from the Groningen field started in 1963. Induced seismicity of the field first was recorded in 1991 (ML 2.4). During the subsequent 10 yr, induced seismicity stayed at a rate of about five events (ML ≥ 1.5) per year, increasing afterwards (e.g. van Thienen-Visser & Breunese 2015). In 2012, the largest event (ML 3.6) occurred, which caused the most damage to date. As a consequence, many studies have been carried out to look at the relationship between induced seismicity and gas depletion. These include predictive approaches based on compaction strain (Bourne et al.2014), and 2-D and 3-D finite-element-based mechanical studies (e.g. Sanz et al.2015; Van den Bogert 2015). These studies focused on the period prior to an 80 per cent production reduction, which was implemented on the 2014 January 17 to the production in the clusters in the central area of the Groningen field (Fig. 2). This area, constituting over 25 per cent of total field production prior to reduction, had shown the majority of induced events over the lifetime of the field (van Thienen-Visser & Breunese 2015). The seismic data of the central area of the field are shown in Fig. 2. The seismic event rate is marked by a statistically significant reduction of seismic events according to a recent study of Nepveu et al. (2016). This strongly suggests that the rate of seismic moment increase is considerably reduced after the production cut. Figure 2. View largeDownload slide Outline of the Groningen field (red line denotes gas water contact), coloured with pressure change predicted from NAM’s reservoir model 1.5 yr prior (left) and after (right) reduction in production. Induced events of the central area within the green polygon is shown by red dots, for the remainder of the field by grey dots (source www.knmi.nl), production clusters by triangles. The hydraulic diffusivity of the reservoir is such that pressure gets equilibrated in between production clusters in the central area within months to half a year, whereas pressure diffusion from the edges of the field to the centre would take over a few years (van Thienen-Visser & Breunese 2015). Figure 2. View largeDownload slide Outline of the Groningen field (red line denotes gas water contact), coloured with pressure change predicted from NAM’s reservoir model 1.5 yr prior (left) and after (right) reduction in production. Induced events of the central area within the green polygon is shown by red dots, for the remainder of the field by grey dots (source www.knmi.nl), production clusters by triangles. The hydraulic diffusivity of the reservoir is such that pressure gets equilibrated in between production clusters in the central area within months to half a year, whereas pressure diffusion from the edges of the field to the centre would take over a few years (van Thienen-Visser & Breunese 2015). The reduction in production in the central area of the Groningen field, the effects of which are closely monitored via the recorded seismicity, provides a unique data set documenting the break in trend of seismicity which can be attributed to relative change in rates of time-dependent compaction and pressure change in the reservoir. Past studies and associated data sets on gas depletion in other fields do not provide such data for the following reasons: A mid-field-life production cut at the scale applied in the central area of the Groningen field, is exceptional. Induced seismicity is not always monitored by a dedicated local network prior to the occurrence of induced events leading to stop of production, such as in the case of the Rotenburg Field in Germany (Dahm et al.2007). Therefore, seismic catalogues generally contain an insufficient number of events for statistical analysis. Gas fields are typically marked by a gradual reduction in production and the associated pressure change decays prior to production stop (e.g. Van Wees et al.2008). This appears to be consistent with a gradual reduction of seismic event rates occurring in the tail-end of production, as observed in the Strachan field in Canada (Baranova et al.1999). In fields where production has stopped, seismicity appears to return to background rates, as observed in the Strachan Field and the Roswinkel and Bergermeer fields in the Netherlands, with 25 and 4 seismic events (M > 1.5) during production and no events for 9 and 8 yr after production stop, respectively, in accordance with absence of natural seismicity (Van Wees et al. 2014). In the aforementioned fields, there was considerable time to allow a gradual cessation of time-dependent compaction prior to production stop, because of the decay in pressure rate changes. This renders these cases not suitable for the analysis pursued here. A strong aquifer support can result in equilibration of reservoir pressure on short timescales and the abundance of natural seismicity can render it difficult to analyse post-production evolution in view of induced seismicity. In summary, the Groningen field is marked by a sharp and prolonged drop in production by 80 per cent and pressure change rates in the central area (Fig. 2), of which the effects in terms of a sharp reduction in seismic event rates are well documented in the seismic catalogue (Nepveu et al. 2016). The field is located in a tectonically stable area, and characterized by limited aquifer pressure support. Therefore, it qualifies as an excellent case study to disentangle the pressure- and time-dependent compaction control on induced seismicity. 3 GEOMECHANICAL MODEL Semi-analytical and finite-element-based geomechanical models demonstrate that a progressive increase in induced events relates to progressive stress perturbations due to depletion, and that the relative magnitude of stress perturbations is strongly dependent on initial stresses, pressure change, reservoir geometry (Segall 1989; Roest & Kuilman 1994; Mulders 2003), reservoir and overburden creep (Orlic & Wassing 2013). In finite-element geomechanical models predicting displacement on faults, the expected evolution of cumulative seismic moment evolution (M0) can be determined via (e.g. Aki & Richards 2002):   $$M0 = \int G\ u\ dS\$$ (1)where S is surface area of faults, u is displacement and G is the shear modulus. The onset of slip is determined by the Coulomb failure function (CFF > 0) (e.g. Zoback 2010):   $${\rm CFF}\, = {{\sigma}_s}\, - \, \mu {{\sigma}_n}\!\!^{'}$$ (2)where σs is shear stress, σn' effective normal stress on the fault and μ the friction coefficient. Slip causes the fault to return to the yield condition (CFF = 0). This can include slip weakening effects and/or rate and state friction effects during rupture events (Ruina 1983; Dieterich 1994; Rutqvist et al.2013; Wassing et al.2014). Part of the slip may occur aseismically. In this simplified approach, the dynamic effects of slip and slip weakening are discarded, and it is assumed that all incremental slip is released seismically. In elastic models, the change in the CFF (ΔCFF) scales linearly to pressure change ΔP through changes in the total stress tensor on the fault (Mulders 2003; Soltanzadeh & Hawkes 2008):   $$\Delta \ {\sigma _{ij}} = \ \alpha ({\gamma _{ij}} - {\delta _{ij}})\ \Delta P$$ (3)where σij are the total stress tensor components, and γij normalized stress path parameters, describing the stress path in terms of change in stress tensor as a function of pressure change. α corresponds to Biot's coefficient. The stress path parameters are dependent on poroelastic parameters, fault geometry and offset, and can be determined semi-analytically from linear poroelasticity or numerically in finite-element models (Mulders 2003; Soltanzadeh & Hawkes 2008). They vary over the fault surface, in case of fault offset, but are static through time under spatially uniform depletion. Prior to the production cut, the Groningen field production strategy targeted at maintaining a uniform pressure change over the field. The central area of the field is marked by high permeability capable to equilibrate pressures in between well clusters in months, rather than years (van Thienen-Visser & Breunese 2015). For this reason, we ignore effects of spatial changes in pressure depletion which can cause stress perturbations (e.g. Altmann et al.2010). In the Groningen field, the pressure change has been decreasing almost linearly during the last 30 yr, prior to the cut in production (e.g. Bourne et al.2014; Sanz et al.2015). The associated stress paths are well capable of predicting an accelerating growth of seismic moment as a function of pressure change (Van Wees et al.2014). Mulders (2003) and Van Eijs et al. (2006) demonstrate that the ΔCFF is amplified with increasing contrast in the elastic properties of reservoir and surrounding rock. 3.1 Adopting time-dependent compaction Insufficient consensus of the exact creep mechanism in the Groningen field is currently available to justify a particular choice for the time-dependent compaction constitutive law (Mossop 2012; Hol et al.2015; Pruiksma et al.2015). For this study, a Kelvin-Chain model has been adopted. This model is based on a specific form of the generalized creep formulation (cf. Bazant et al.1988; De Borst & Van den Boogaard 1994):   $${\boldsymbol{\varepsilon}}(t) = \mathop \int \nolimits_0^t J\left( {t,\tau } \right){\boldsymbol{C\dot{\sigma }}}\left( \tau \right)d\tau$$ (4)ε(t) is strain tensor at time t, $${\boldsymbol{\dot{\sigma }}}$$ is stress rate tensor, J(t, τ) is creep function; τ  is time of load application and the components of C are defined by (De Borst & Van den Boogaard 1994):   $${C_{ijkl}} = \frac{1}{2}\left( {1 + \nu } \right)( {{\delta _{ik}}{\delta _{jl}} + {\delta _{il}}{\delta _{jk}}} ) - \ \nu \ {\delta _{ij}}{\delta _{kl}}$$ (5) In the Kelvin-Chain model, a single creep element in series with the elastic response is adopted (Fig. 4). The creep function, J(t, τ) corresponds to:   $$J\left( {t,\tau } \right) = {E^{ - 1\ }}\ + {E_{1\ }}^{ - 1}\left[ {1 - {e^{ - \left( {\frac{{t - \tau }}{{{\lambda _1}}}} \right)}}} \right]$$ (6) The spring with stiffness E allows for an elastic stress response, which should be in accordance to the static Young's modulus. Lab measurements document E typically in the range of 8–25 GPa for rock samples which are representative for the Groningen gas field (Hol et al.2015; Pruiksma et al.2015). The Kelvin-Chain component combining the spring with constant E1 and the dashpot with Newtonian viscosity η1 represents the time-dependent compaction. For uniaxial compaction, the creep component εcz(t) is constructed by:   $$\varepsilon _z^c(t) = \int_{0}^{t}{{\left[ {1 - {e^{ - \left( {\frac{{t - \tau }}{{{\lambda _1}}}} \right)}}} \right]}}{C_{{\rm{tinf}}}}\dot{\Delta }P\left( \tau \right)d\tau$$ (7)where $${C_{{\rm{tinf}}}} = \ \frac{{( {1 - v - 2{v^2}} )/( {1 - v} )}}{{{E_{1\ }}}}$$. It follows that the uniaxial Kelvin-Chain model corresponds to the decay model for time-dependent compaction introduced by Mossop (2012):   $$\ {\varepsilon ^c}\!\!_z\ (t) = f'\ *g\ = \ {{\cal F}^{ - 1}}\left[ {{\cal F}\left[ {f'\ } \right]\ {\cal F}\left[ g \right]} \right]$$ (8)with   \begin{equation*}f'(t) = \frac{d}{{dt}}\ \left[ {1 - {e^{ - \left( {\frac{{t - \tau }}{{{\lambda _1}}}} \right)}}} \right]\ = \frac{1}{{{\lambda _1}}}{e^{ - t/{\lambda _1}}}\end{equation*}  \begin{equation*}g\ (t) = {C_{{\rm{tinf}}}}{\rm{\ \Delta }}P(t)\end{equation*} where $${\cal F}$$ denotes the Fourier Transform in the time domain and $${{\cal F}^{ - 1}}$$ its inverse, ΔP represents pressure change and $${\lambda _1} = \frac{{{\eta _1}}}{{{E_1}}}$$ is the relaxation time of time-dependent compaction. The Kelvin-Chain model cannot account for irreversible deformation, and stress paths in alternative constitutive models will deviate from this paper's model results. However, the relative change in stress response upon cessation of depletion is not likely to change significantly, as this is primarily driven by the break in loading mechanism after an arrest of pressure changes. Creep scenarios: in order to evaluate the effects of time-dependent compaction, three different scenarios of E/E1 ratios of final time-dependent (creep) strain and elastic strain are considered in the Kelvin-Chain model (Table 2). The model parameters have been chosen such that the total final strain in all models is the same. The models with E/E1 = 0.5 and 1 time-dependent compaction closely match laboratory experiments (e.g. Hol et al.2015; Pruiksma et al.2015), whereas the E/E1 = 5 case mimics the compaction model of Mossop, where time-dependent compaction is more dominant than elastic deformation. The adopted models predict relatively large creep in the first years after a sudden halt in pressure change compared to rate type and linear compaction models (e.g. Pruiksma et al.2015), for the same amount of final creep. In this sense, the Kelvin-Chain model can be considered most pessimistic in terms of reservoir strains and stress changes associated with these strains. Pseudo-elastic approach: under uniaxial loading conditions, the total strain for the Kelvin Chain is:   $${\varepsilon _z}(t)\ = \ {\varepsilon ^c}\!\!_z(t) + {\varepsilon ^e}\!\!_z(t)$$ (9)with direct elastic strain εez(t) corresponding to the instantaneous reacting spring in the Kelvin Chain:   $${\varepsilon ^e}\!\!_z(t) = {C_m}\,\,\Delta P(t),\quad{C_m} = \frac{{(1 - v - 2{v^2})/(1 - v)}}{E})$$ (10) The resulting stress–strain relationship can be written as:   $${\varepsilon _z}(t) = \widetilde {{C_m}}{\rm{\ }}\left( {\rm{t}} \right){\rm{\ \Delta }}P(t)$$ (11)with   \begin{equation*}\widetilde {{C_m}} = {\varepsilon _{{\rm{ratio}}}}(t)\ {C_m}\end{equation*}  \begin{equation*}{\varepsilon _{{\rm{ratio}}}}(t) = \frac{{{\varepsilon _z}(t)}}{{{\varepsilon ^e}\!\!_z(t)}}\end{equation*} At infinite time, εratio(t) corresponds to E/E1. The scaling of the stress–strain tensor relationship is scalar in the creep function, invariant for v. Therefore, in the stress–strain tensor relationship in a finite-element formulation, εratio(t) can be used to obtain a time-dependent pseudo-compaction coefficient or corresponding pseudo-stiffness for the reservoir zone (Young's modulus referred to as Ep(t)). This allows us to relate relative stress responses as a function of time-dependent compaction from a pseudo-elastic analysis adopting the pseudo-stiffness for reservoir and E for surrounding rock. A benefit of this pseudo-elastic approach is it allows numerically to be coupled in finite-element models with Mohr–Coulomb failure, whereas for the Kelvin-Chain model this is not supported in the workflow we use. 3.2 Model geometry and parametrization A plane-strain elastic finite-element model of a depleting gas reservoir is constructed and used to analyse the sensitivity of total stress scaling and ΔCFF to different Kelvin-Chain parameters. For the finite-element modeling, DIANA FEA (http://dianafea.com) is used, a commercially available finite-element code that is widely used for geomechanical analysis in oil and gas reservoirs, including the analysis of Mohr–Coulomb failure and slip. More details of the finite-element modeling procedure are given in Orlic & Wassing (2013). The numerical model incorporates two reservoir compartments, separated by a fault with a throw of half the reservoir thickness (Fig. 5 and Table 1), which can be considered representative for major fault zones in the central area of the Groningen field, where most seismic events have occurred (Fig. 2; Bourne et al.2014). The depth and thickness range, and underlying parametrization is in accordance with earlier 2-D models for the Groningen gas field (e.g. Van den Bogert 2015). The plane strain finite-element mesh models a 2-D section of 6 km width and 6 km depth. The model consists of roughly 135 000 nodes which comprise approximately 44 000 quadratic quadrilateral plane-strain elements representing the bulk rock and ∼250 quadratic quadrilateral interface elements representing the fault. The mesh is finest around the fault with a resolution of 1 m at reservoir depth and gets gradually coarser away from the fault and the reservoir (Fig. 3). The model resolution has been chosen sufficiently high to capture high bending stresses close to the fault and at the interface of reservoir and adjacent rock. Figure 3. View largeDownload slide Geometry and finite-element mesh for the model setup. Figure 3. View largeDownload slide Geometry and finite-element mesh for the model setup. Table 1. Parameters of the geomechanical model, based on earlier studies for the Groningen field (cf. Mulders 2003; Van Wees et al. 2014; Van den Bogert 2015). Parameter  Value  Unit  E  18  GPa  V  0.2  –  Top reservoir  2900  m  Thickness reservoir  150  m  Fluid density  1150  kg m−3  Rock density  2260  kg m−3  Friction angle  30  deg  Fault dip  70  deg  Depletion  25  MPa  K0eff  0.45  –  Parameter  Value  Unit  E  18  GPa  V  0.2  –  Top reservoir  2900  m  Thickness reservoir  150  m  Fluid density  1150  kg m−3  Rock density  2260  kg m−3  Friction angle  30  deg  Fault dip  70  deg  Depletion  25  MPa  K0eff  0.45  –  View Large The plane-strain elements evaluate stresses and pressures on eight integration points per element, while the interface elements contain six integration points per element. The model is laterally constrained (zero displacement in horizontal direction) on the sides and vertically constrained (zero displacement in vertical direction) at the bottom, and subject to gravity loading. This effectively leads to a stress response on the sides of the model which corresponds to an infinite lateral extension of the reservoir. For each loading step, the systems of equations are solved using an iterative solution method which applies incomplete LU-decomposition preconditioning. For the Kelvin-Chain creep model, the elements representing the reservoir are governed by a Kelvin-Chain constitutive model, while the surrounding rock is modeled to be fully linearly elastic. In order to prevent the development of shear stresses and shear strains in the model due to the viscoelastic response of the reservoir to the initialization load, the model is initialized with ρrock = ρfluid. This results in a model with increasing pore pressure with depth, but without effective stresses. The reservoir is then depleted, resulting in non-zero effective stresses. The stress changes of this model are superimposed on the results of a fully elastic model which represents the initial stress state. For the initial stress state, the reservoir is assumed overpressured relative to hydrostatic conditions (Bourne et al.2014). The initial ratio of effective horizontal-to-vertical stress is set to K0eff = 0.45 (cf. Van Wees et al.2014). In the pseudo-elastic approach, we adopt the pseudo-stiffness Ep(t) for the reservoir, resulting in contrasting stiffness of the reservoir and the surrounding rock. The models are initialized in two steps to prevent development of shear stresses and shear strains in the pre-depletion state. First, the model is initialized with a homogeneous stiffness all throughout the model and a prescribed K0eff. The resulting stress state is then transferred to a model which contains the stiffness contrast. The final depletion in the reservoir is 25 MPa, approximately in agreement with pressure depletion at present day in the field. Depletion is implemented to occur simultaneously, at the same rate, in both reservoir compartments in yearly load steps. The pressure depletion of 25 MPa is built-up linearly in a time period of 60 yr, followed by 40 yr of no pressure change. The decrease in pore pressure is the same for all values of the horizontal coordinate, in agreement with high diffusion rates observed in the Groningen gas field (van Thienen-Visser & Breunese 2015), marked by a delay time of a couple of months to equilibrate pressures from well cluster in the central area The linear pressure depletion rate has been constrained by the observed linear pressure reduction in the last 30 yr prior to reduction in production (Bourne et al.2014). The time period of pressure depletion has been extended to 60 yr relative to approximately 50 yr of actual production in the field, to compensate for higher pressure rate changes in the first 20 yr of production (Bourne et al.2014). The fault zone is assumed to be subject to the reservoir fluid pressure when reservoir rock is present on both sides of the fault. For the portion of the fault where reservoir rock and surrounding rock are juxtaposed, two scenarios for pressure in the fault are considered: actual reservoir pressure and ambient (initial reservoir) pressure. Obviously, a single fault geometry cannot capture the complex response expected from many mapped faults in the area with variable throw, dip and strike (Fig. 2; Sanz et al.2015). Uncertainty in the location of seismic events hampers an allocation of seismic events to particular faults. Fault strike, dip, throw, natural stress variability and variability in mechanical properties, can all have a significant control on in situ stress and spatial and temporal variability in seismicity (e.g. Van Wees et al.2014; Sanz et al.2015). However, these aspects are not expected to have a prime control on the relative change in seismic response upon production stop, as they do not affect strongly the ΔCFF as a function of pressure change. 3.3 Coulomb failure function (ΔCFF) evolution The evolution of ΔCFF is constructed from the resulting stress changes of the finite-element model for the three time-dependent compaction scenarios from Table 2, with creep to elastic deformation ratios (E/E1) of 0.5, 1 and 5 at infinite time. Fig. 4 shows an example of the adopted pressure change, and resulting effective normal stress and shear stress changes for E/E1 = 1 for both the creep (Kelvin Chain) and pseudo-elastic finite-element model. Fig. 4 also shows the ΔCFF for all three scenarios. During production, the ΔCFF increases almost linearly with production rate. The total ΔCFF increase differs, depending on the relative contribution of time-dependent compaction and is in close agreement with the stress response in the pseudo-elastic approach (eq. 8). Figure 4. View largeDownload slide Finite-element model solutions for Kelvin-Chain creep (solid lines) and pseudo-elastic models (dashed lines). Panels (a) and (b) for E/E1 = 1 (a) effective normal stress on the fault assuming reservoir pore pressure, (b) shear stress on the fault. (c)–(e) Change in Coulomb Failure Function for E/E1 = 0.5, 1 and 5, respectively. The sign of shear stress is positive if shear sense is normal faulting. The time in years relative to production stop is indicated by t. The depth extent of the reservoir zones is indicated on both sides of the fault in grey. Figure 4. View largeDownload slide Finite-element model solutions for Kelvin-Chain creep (solid lines) and pseudo-elastic models (dashed lines). Panels (a) and (b) for E/E1 = 1 (a) effective normal stress on the fault assuming reservoir pore pressure, (b) shear stress on the fault. (c)–(e) Change in Coulomb Failure Function for E/E1 = 0.5, 1 and 5, respectively. The sign of shear stress is positive if shear sense is normal faulting. The time in years relative to production stop is indicated by t. The depth extent of the reservoir zones is indicated on both sides of the fault in grey. Table 2. Viscoelastic parameters in the Kelvin-Chain model (Fig. 3). Scenario  E (GPa)  E1 (GPa)  λ1 = η1/E1 (yr)  50 per cent creep  13.5  27  7.3  100 per cent creep  18  18  7.3  500 per cent creep  54  10.8  7.3  Scenario  E (GPa)  E1 (GPa)  λ1 = η1/E1 (yr)  50 per cent creep  13.5  27  7.3  100 per cent creep  18  18  7.3  500 per cent creep  54  10.8  7.3  Notes: Viscosity η1 is chosen such that it agrees with a relaxation time λ1 = 7.3 yr (Mossop 2012). View Large In the 40 yr after production, ΔCFF increases only moderately by a few per cent relative to the ΔCFF build-up during production, and the curves are hard to discriminate in Fig. 4. The moderate amount of total ΔCFF change after production stop, and the abrupt change in ΔCFF evolution, are easily understood from the fundamentally different mechanisms contributing to stress changes. Prior to production stop, both the pressure change (eq. 3) and the progressive growth of εratio(t) is contributing to ΔCFF, whereas after production stop the only contribution to ΔCFF is the moderate increase of εratio(t) as a function of progressive compaction (cf. eq. 11). In Fig. 4, all the creep models show that ΔCFF has a distinct peak relative to the pseudo-elastic approach at the structural interfaces of the reservoir and surrounding rock. These are related to the creep causing higher and localized deviatoric stresses at the structural interfaces, during depletion. It highlights the limitations of the pseudo-elastic approach to reproduce the exact value of stress changes at structural interfaces. Fig. 5 shows the relative change in ΔCFF, after production stop. Here, only the creep model results are presented, as it was recognized that the pseudo-elastic model is not well capable of reproducing the relative ΔCFF changes of the Kelvin-Chain model, in the first years immediately after production stop. The change in ΔCFF after production stop is lower than the change in ΔCFF just prior to production stop for all cases considered—with the most significant difference experienced for the strongest time-dependent compaction scenarios. Figure 5. View largeDownload slide (a)–(c) Change in Coulomb failure function on the fault around productions for E/E1 = 0.5, 1 and 5, respectively, assuming reservoir pore pressure on the fault. The time in years relative to production stop is indicated by t. Figure 5. View largeDownload slide (a)–(c) Change in Coulomb failure function on the fault around productions for E/E1 = 0.5, 1 and 5, respectively, assuming reservoir pore pressure on the fault. The time in years relative to production stop is indicated by t. Figs 4 and 5 consistently show that ΔCFF is monotonously increasing before and after production stop for the portion of the faults, which is subject to failure. This indicates that the stress path remains similar. Furthermore as the total ΔCFF change after production stop is so small compared with the production period, it is expected that potential unsolicited effects of reversibility of creep, which is inherent of the Kelvin-Chain model, are not significant. Summarizing, for all creep scenarios considered, the predicted increase of ΔCFF after a production stop is few per cent relative to the ΔCFF which had been built-up during production, marked by a strong break in stress build-up rates, at the production cut. The scenarios also show a very quick decay in the build-up of ΔCFF in the first few years after the end of production. Furthermore, the models show a close resemblance of stress responses of the Kelvin-Chain model and the pseudo-elastic approach, except at structural interfaces and in the years immediately after production stop. The break in stressing rate is in agreement with the observed break in seismic event rates (Fig. 2; Nepveu et al.2016). 4 SEISMIC MOMENT The previous sections highlight the generic relationship between ΔCFF with variations in time-dependent compaction. The bearings of modeled ΔCFF changes on seismicity are evaluated in this section, in terms of the relative effect on seismic moment increase. As mentioned before, the finite-element model approach used here does not support a combined Kelvin-Chain and Mohr–Coulomb material model, for calculation of fault slip and moment (cf. eq. 1). To overcome this limitation, the relative growth of seismic moment is determined by a novel semi-analytical approach for the construction of the progressive increase of slip area and displacement of the fault, based on the elastic stress solution only. 4.1 Semi-analytical model to resolve moment from the elastic stress solution For a circular fault area, a closed relationship exists between area S, radius r, average slip u, stress drop Δσ and seismic moment M0 (Eshelby 1957; Keilis-Borok 1959; Udias et al.2014):   $$M0 = \frac{16}{7\pi} \Delta \sigma rS$$ (12) This equation is commonly used by seismologists to estimate Δσof an earthquake after determining M0 and source dimensions of the rupture surface (e.g. Udias et al.2014), and has been generalized for elliptic surfaces of fault slip (Madariaga 1979). In this case, the slip and moment are scaled by a correction factor which takes into account both the ellipticity and the slip direction on the fault. For plane-strain dip-slip faulting conditions, the correction factor cancels the term $$\frac{16}{7\pi }$$. Considering a slipping fault section with length l in strike and dip direction, and adopting S = l2, and $$r = \frac{l}{{\sqrt {\rm{\pi }}}}$$ for the radius of the equivalent circular fracture and substituting in eq. (9) gives the fault moment density [N] per m strike:   $$M{0_m} = \ \Delta \sigma \frac{{\ {l^2}}}{{\sqrt \pi }}$$ (13) In the approach adopted here, moment is determined from the determination of slip length and average stress drop from the elastic stress solution. Initially, l = l0 is estimated to correspond to the portion of the fault where CFF > 0, and CFF corresponding to Δσ. In reality, the stress drop would cause the fault length to grow progressively beyond the fault tips, with a stress redistribution resulting in CFF ≤ 0 on the slipping portion of the faults. We designed an algorithm for the estimation of the increased slip length based on the assumption that the elastic energy is redistributed by shear. To this end, the CFF is discretized in patches along the dip of the fault at a resolution of Δx = 1 m, and l and Δσ are solved by the following steps: Determine the range of patches [istart.iend] where CFFi > 0. For this range calculate the integral of excess Coulomb stress $$Y_{\rm{sum\ }} = \mathop \sum \nolimits_{{i_{{\rm{start}}}}}^{{i_{{\rm{end}}}}} {\rm{CF}}{{\rm{F}}_i}{\rm{\ \Delta \ }}x$$. Set Yistart = Yiend =  0.5 Ysum Beyond the tips of the fault distribute the excess stress Yistart and Yiend following a function of the distance di of the centre of each patch to the nearest tip of the fault (corresponding to istart or iend). Increment CFFi with ai =  Adi − 1.3. A is chosen such that the integral of the function equates the fraction of Y to be distributed beyond the tip under consideration. In the range of patches [istart.iend] set CFFi = 0. For incremented CFFi in the other patches, include patch i in the range of slipping patches if it results in failure (CFFi > 0) and recalculate the integral of excess Coulomb stress $$Y{\rm{\ }} = \mathop \sum \nolimits_{{i_{{\rm{start}}}}}^{{i_{{\rm{end}}}}} {\rm{CF}}{{\rm{F}}_i}{\rm{\ \Delta \ }}x$$. Divide Y over Yistart and Yiend according the portion originated from each of the tips of the fault. If Y > 0, repeat steps 2–4. Set $$l{\rm{\ }} = {\rm{\ }}( {{i_{{\rm{end}}}} - {i_{{\rm{start}}}}} ){\rm{\Delta }}x,\Delta \sigma {\rm{\ }} = \frac{{Y{\rm{sum}}}}{l}{\rm{\ }}$$. In the algorithm, a key assumption is that the area integral of Coulomb stress change remains unchanged upon slip (cf. Okada 1992; Wassing et al.2014). The characteristics of the power function for the Coulomb stress redistribution in step 2 have been determined from a numerical analysis of stress response for l in the range of 10–100 s of metres, for a steeply dipping fault at 3000 m depth, infinite strike length and adopting the analytical solution of slip for a constant stress drop (e.g. Madariaga 1979). The stress response has been calculated in a regular grid with 10 m spacing from an analytical solution of stress for rectangular sources of slip following the grid resolution (cf. Madariaga 1979; Okada 1992). The average input slip and output stress response of each patch of the grid is determined using two-point Gauss–Legendre Integration, which gave a numerically convergent result relative to a finer grid resolution for stress distribution. The stress redistribution is not sensitive to elastic parameters, nor to the exact orientation of the slip plane in the stress field, as can be derived from the stress equations (Okada 1992). Fig. 6 depicts the magnitude of the stress increase decaying away from the fault tip. It demonstrates that the redistribution of stress affects a considerable area, and can be well represented by a power function of distance from the tip for the decay. Figure 6. View largeDownload slide Predicted Coulomb stress change for analytical stress response (cf. Okada 1992, friction angle cf. Table 1), as a function of distance from the fault tip, for an elliptical displacement along a normal fault oriented favourably in the stress field and with a maximum of 0.01 m displacement for l = 30, 110 and 210 m. The dashed line represents the best fit x−1.3 line for l = 110 m (cut-off at x = 300 m). Figure 6. View largeDownload slide Predicted Coulomb stress change for analytical stress response (cf. Okada 1992, friction angle cf. Table 1), as a function of distance from the fault tip, for an elliptical displacement along a normal fault oriented favourably in the stress field and with a maximum of 0.01 m displacement for l = 30, 110 and 210 m. The dashed line represents the best fit x−1.3 line for l = 110 m (cut-off at x = 300 m). In order to test the robustness of the method for estimating seismic moment, various benchmarks comparing finite-element models resolving slip and slip length along faults incorporating Mohr–Coulomb failure (cf. Orlic & Wassing 2013) have been run, applying eq. (1) and comparing these to the results of the semi-analytical approach. For the model setup, the initial geometry, properties, boundary and initial conditions used for the previous models are applied. Static elastic properties through time have been adopted, with offset of the fault in the range (0.5–1)-D, and both ambient and reservoir pressure in the faults have been considered. Two model results for reservoir pressure are shown for 0.5-D and 1-D displacement in Fig. 7. This figure shows an excellent correspondence of fault moment density obtained by the finite-element method and the analytical method. In all model runs, the moment predicted from the semi-analytical model runs show less than 5 per cent deviation from the moment in the finite-element solutions, whereas the uncorrected l = l0 would significantly underestimate the moment by 25–70 per cent. From this, it is concluded that the semi-analytical approach provides a good estimator of moment growth from an elastic stress solution. Figure 7. View largeDownload slide Comparison of moment predicted by finite-element analyses incorporating Mohr–Coulomb failure and the semi-analytical model, for 0.5-D and 1-D offsets prior to production stop, adopting uniform elastic properties (cf. Table 1, reservoir pressure in fault). Figure 7. View largeDownload slide Comparison of moment predicted by finite-element analyses incorporating Mohr–Coulomb failure and the semi-analytical model, for 0.5-D and 1-D offsets prior to production stop, adopting uniform elastic properties (cf. Table 1, reservoir pressure in fault). The strength of the semi-analytical approach is that it provides a direct relationship between relative changes in ΔCFF and moment evolution, which is not biased to potential changes in response due to variability of in situ stress, specific orientation of faults and contrasts in elastic and frictional properties. Furthermore, it avoids potentially ambiguous results of finite-element slip with respect to the shear modulus to be used in eq. (1) to construct moment, if the model is marked by contrasting material behaviour on both sides of the fault. 4.2 Post-production time-dependent compaction and seismic moment increase In Fig. 8, the evolution of the Kelvin-Chain compaction and fault seismic moment density determined from the Coulomb failure function by the algorithm in Section 4.1 is plotted, normalized to its value at the end of production. Furthermore, the moment density has been constructed in accordance to the heuristic correlation of compaction strain and the logarithm of seismic moment increase (cf. Bourne et al.2014, median in their fig. 18a). These have been constructed for a reservoir segment of 2 km wide in map view, representative for average fault spacing in the central area. The Kelvin-Chain compaction in Fig. 8 for the E/E1 = 0.5 and 1 scenarios shows a discontinuity after production stop. This is supported by a clear discontinuity in GPS subsidence data in the central area of the Groningen field occurring nine weeks after production stop (Pijpers & Van der Laan 2016). Markedly, the ratio of seismic moment increase after and before production stop in the geomechanical models is 2–10 times lower and more than 10 times lower than those predicted from the approach using compaction as proxy (Bourne et al.2014), for the first few years and decades after shut-in, respectively (Fig. 8). This holds for all three time-dependent compaction scenarios that are considered here, as well as for both the reservoir and ambient pressure conditions on the fault. For a model that does not incorporate creep, the normalized stress changes and increase in seismic moment would be very similar to the E/E1 = 0.5 scenario and would arrest when production stops. The geomechanical predictions of all creep scenarios are in accordance with the strong reduction of induced events after 80 per cent cut in production rates in the central area in the Groningen field (cf. Fig. 2, Nepveu et al.2016). For reservoir pressure condition on the fault, the predicted rate of fault moment density increase for the E/E1 = 0.5 and 1 scenarios is about 0.5 GN yr−1, 1.5 yr prior to production stop. The central area is marked by about 5 × 1013 N m moment increase per year in the same period based on the seismic events. In the geomechanical model, 100 km of fault length would suffice to generate the observed moment rate. For the ambient pressure scenario less than 5 km of fault length would be sufficient. Cumulative fault length in the central area is in the order of 50 km (Fig. 2). The evolution of seismic moment prior to production stop is showing less accelerated growth as observed in the field (Bourne et al.2014). This may be explained by the fact that the 2-D model fails to capture important effects contributing to increased growth, including spatial variability in fault strike, dip, reservoir offset and thickness and initial stress. These can jointly contribute to a progressively increasing fault length involved in slip and would result in a more convex pattern of seismic moment growth. Figure 8. View largeDownload slide Left vertical axis: compaction of the reservoir for different creep models. On the right vertical axis: seismic moment density normalized to the seismic moment at t = 0 using the Bourne et al. (2014) heuristic correlation of compaction and seismic moment (α in the legend refers to the moment partitioning coefficient at t = 0, see Bourne et al. (2014) and the text for details) and using the analytical method described in this paper (moment density in the legend refers to the moment density at t = 0). Figure 8. View largeDownload slide Left vertical axis: compaction of the reservoir for different creep models. On the right vertical axis: seismic moment density normalized to the seismic moment at t = 0 using the Bourne et al. (2014) heuristic correlation of compaction and seismic moment (α in the legend refers to the moment partitioning coefficient at t = 0, see Bourne et al. (2014) and the text for details) and using the analytical method described in this paper (moment density in the legend refers to the moment density at t = 0). 5 DISCUSSION AND CONCLUSIONS The Kelvin Chain creep and the geometry of the model do not allow an in-depth assessment of the sensitivities of the results to variability in geometrical aspects of the reservoir and alternative creep behaviour. However, the prediction of a significant reduction of seismicity after cessation of production is expected to hold as well for alternative geometries and constitutive laws, as the first-order effects are determined by the arrest in pressure change (eq. 3). In the model developed, it is assumed that pressure does not change after production stop. Because of continued production in the other production clusters of the Groningen field (Fig. 2), in two years it is expected that reservoir pressures will be further reduced in the central area, due to hydraulic diffusion (van Thienen-Visser & Breunese 2015), which is likely to increase the likelihood of seismicity, compared to the present state. The prediction of the model, building from the relative increase in expected moment, needs to be interpreted with care in view of earthquake processes. A critically stressed fault at the end of the production period, which is incrementally loaded by an additional few per cent of CFF predicted by this model after stop of production, could be sufficient to trigger a relatively large event (i.e. M > 3). Therefore, it is important to note that the model does not preclude the possibility of a large event. Concluding, our geomechanical stress-based models incorporating reservoir creep, predict that seismic moment evolution after a production stop is subject to a very moderate increase, in the order of 2–10 times and more than 10 times less than the values predicted by the approach correlating seismic moment to progressive compaction strain (Bourne et al.2014), for the first years and decades following cessation of production, respectively. This difference is easily understood by the fact that the strain-based model of Bourne et al. (2014) ignores stress rate changes, which occur at a production stop. Our stress-based model is in agreement with the low seismicity encountered in the central area of the Groningen field, where production has been reduced by 80 per cent, on 2014 January 17. The model findings support scope for an immediate mitigation of induced seismicity through adjusting production scenarios, as the seismicity rates are expected to be correlated primarily with rates of pressure change. The findings of this study may well be generalized to other depleting reservoirs in similar structural settings, as the strong jump of stress loading behaviour on faults as a function of arrest of pressure changes (eq. 3) is universal, and prolonged stress effects of time-dependent compaction are typically one order of magnitude lower than those of pressure change. A number of issues were beyond the scope of this study but need to be addressed to augment to the robustness of our findings. These include: more realistic constitutive laws for reservoir creep; adoption of the full 3-D complexity of reservoir structure and of pressure depletion; pressure diffusion and creep in the over and underburden of the reservoir and effects of the rupture process in terms of a rate and state friction model. ACKNOWLEDGEMENTS The data used in this manuscript are publically available: seismicity data from the Royal Netherlands Meteorological Institute KNMI (www.knmi.nl); monthly production data from production wells in the Groningen field, available from the Geological Survey of the Netherlands TNO (www.nlog.nl) and numerous data and reports on the Groningen gas field (namplatform.nl). The authors would like to thank Jaap Breunese, Dries Hegen, Bogdan Orlic and Brecht Wassing for constructive comments on earlier versions of the manuscript. The manuscript benefited considerably from constructive comments from Alexander Droujinine and an anonymous reviewer REFERENCES Aki K., Richards P., 2002. Quantitative Seismology , University Science Books. Altmann J.B., Müller T.M., Müller B.I.R., Tingay M.R.P., Heidbach O., 2010. Poroelastic contribution to the reservoir stress path, Int. J. Rock Mech. Min. Sci. , 47( 7), 1104– 1113. https://doi.org/10.1016/j.ijrmms.2010.08.001 Google Scholar CrossRef Search ADS   Baranova V., Mustaqeem A., Bell S., 1999. A model for induced seismicity caused by hydrocarbon production in the Western Canada Sedimentary Basin, Can. J. Earth Sci. , 36( 1), 47– 64. https://doi.org/10.1139/e98-080 Google Scholar CrossRef Search ADS   Bazant Z.P., Dougill J., Huet C., Tsubaki T., Wittmann F.H., 1988. Material models for structural creep analysis, in Mathematical Modelling of Creep and Shrinkage of Concrete, ed. Bazant Z.P.J., Wiley & Sons, Chichester, England, pp. 275– 310. Bourne S., 2014. A seismological model for earthquakes induced by fluid extraction from a subsurface reservoir, J. geophys. Res. , 119, 8991– 9015. Google Scholar CrossRef Search ADS   Chang C., Mallman E., Zoback M., 2014. Time-dependent subsidence associated with drainage-induced compaction in Gulf of Mexico shales bounding a severely depleted gas reservoir, AAPG Bull. , 98( 6), 1145– 1159. https://doi.org/10.1306/11111313009 Google Scholar CrossRef Search ADS   Dahm T., Krüger F., Stammler K., Klinge K., Kind R., Wylegalla K., Grasso J.-R., 2007. The 2004 Mw 4.4 Rotenburg, Northern Germany, earthquake and its possible relationship with gas recovery, Bull. seism. Soc. Am. , 97( 3), 691– 704. https://doi.org/10.1785/0120050149 Google Scholar CrossRef Search ADS   De Borst R., Van den Boogaard A.H., 1994. Finite-element modelling of deformation and cracking in early-age concrete, J. Engin. Mech. , 120( 12), 2519– 2534. Google Scholar CrossRef Search ADS   De Waal J., 1986. On the rate type compaction behaviour of sandstone reservoir rock, PhD thesis , Technical University Delft, Delft, Netherlands. Dieterich J., 1994. A constitutive law for rate of earthquake production and its application to earthquake clustering, J. geophys. Res. , 99( B2), 2601– 2618. https://doi.org/10.1029/93JB02581 Google Scholar CrossRef Search ADS   Eshelby J.D., 1957. The determination of the elastic field of an ellipsoidal inclusion, and related problems, Proc. R. Soc. A: Math. Phys. Eng. Sci. , 241( 1226), 376– 396. https://doi.org/10.1098/rspa.1957.0133 Google Scholar CrossRef Search ADS   Fokker P.A., Van Thienen-Visser K., 2016. Inversion of double-difference measurements from optical leveling for the Groningen gas field, Int. J. Appl. Earth Observ. Geoinform. , 49, 1– 9. https://doi.org/10.1016/j.jag.2016.01.004 Google Scholar CrossRef Search ADS   Hettema M.H.H., Schutjens P.M.T.M., Verboom B.J.M., Gussinklo H.J. et al.  , 2000. Production-induced compaction of a sandstone reservoir: the strong influence of stress path, SPE Reserv. Eval. Eng. , 3( 04), 342– 347. https://doi.org/10.2118/65410-PA Google Scholar CrossRef Search ADS   Hol S., 2015. Long-term compaction behavior of Permian sandstones - an investigation into the mechanisms of subsidence in the Dutch Wadden Sea, in 49th US Rock Mechanics/Geomechanics Symposium , American Rock Mechanics Association, San Francisco, CA. Keilis-Borok V., 1959. On estimation of the displacement in an earthquake source and of source dimensions, Ann. Geofys. , 12, 205– 214. Kostrov V., 1974. Seismic moment and energy of earthquakes, and seismic flow of rock, Phys. Solid Earth , 1, 13– 21. Madariaga R., 1979. On the relation between seismic moment and stress drop in the presence of stress and strength heterogeneity, J. geophys. Res. , 84( B5), 2243– 2250. https://doi.org/10.1029/JB084iB05p02243 Google Scholar CrossRef Search ADS   Marketos G., Govers R., Spiers C.J., 2015. Ground motions induced by a producing hydrocarbon reservoir that is overlain by a viscoelastic rocksalt layer: a numerical model, Geophys. J. Int. , 203( 1), 198– 212. https://doi.org/10.1093/gji/ggv294 Google Scholar CrossRef Search ADS   McGarr A., 2014. Maximum magnitude earthquakes induced by fluid injection, J. geophys. Res. , 119, 1008– 1019. Google Scholar CrossRef Search ADS   Mossop A., 2012. An explanation for anomalous time dependent subsidence, in 46th US Rock Mechanics/Geomechanics Symposium , American Rock Mechanics Association, Chicago, IL. Mulders F., 2003. Modelling of stress development and fault slip in and around a producing gas reservoir, PhD thesis , Technical University of Delft, Delft, The Netherlands. Nepveu M., Van Thienen-Visser K., Sijacic D., 2016. Statistics of seismic events at the Groningen field, Bull. Earthq. Eng. , 14( 12), 3343– 3362. Google Scholar CrossRef Search ADS   Okada Y., 1992. Internal deformation due to shear and tensile faults in a half-space, Bull. seism. Soc. Am. , 82, 1018– 1040. Orlic B., Wassing B.B.T., 2013. A study of stress change and fault slip in producing gas reservoirs overlain by elastic and viscoelastic caprocks, Rock Mech. Rock Eng. , 46, 421– 435. https://doi.org/10.1007/s00603-012-0347-6 Google Scholar CrossRef Search ADS   Pijpers F., Van der laan F., 2016. Trend changes in ground subsidence in Groningen. Available at: https://www.sodm.nl/binaries/staatstoezicht-op-de-mijnen/documenten/publicaties/2016/06/21/23—cbs—trend-changes-in-ground-subsidence-in-groningen/23-cbs-trend-changes-in-ground-subsidence-in-groningen-update-may-2016-052016.pdf, last accessed 2 October 2017. Statistics Netherlands | 2016 Scientific paper. Pruiksma J., 2015. Isotach formulation of the rate type compaction model for sandstone, Int. J. Rock Mech. Min. Sci. , 78, 127– 132. Roest J., Kuilman W., 1994. Geomechanical analysis of small earthquakes at the Eleveld gas reservoir, Rock Mechanics in Petroleum Engineering, 29–31 August, Delft, Netherlands, doi:10.2118/28097-MS. Ruina A., 1983. Slip instability and state variable friction laws, J. geophys. Res. , 88, 10 359– 10 370. https://doi.org/10.1029/JB088iB12p10359 Google Scholar CrossRef Search ADS   Rutqvist J., Rinaldi A.P., Cappa F., Moridis G.J. et al., 2013. Modeling of fault reactivation and induced seismicity during hydraulic fracturing of shale-gas reservoirs, J. Petrol. Sci. Eng. , 107, 31– 44. https://doi.org/10.1016/j.petrol.2013.04.023 Google Scholar CrossRef Search ADS   Sanz P. et al., 2015. Geomechanical analysis to evaluate production-induced fault reactivation at Groningen gas field, in SPE Annual Technical Conference and Exhibition, 28–30 September, Houston, Texas, USA., doi:10.2118/174942-MS. Segall P., 1989. Earthquakes triggered by fluid extraction, Geology , 17, 942– 946. https://doi.org/10.1130/0091-7613(1989)017%3c0942:ETBFE%3e2.3.CO;2 Google Scholar CrossRef Search ADS   Soltanzadeh H., Hawkes C.D., 2008. Semi-analytical models for stress change and fault reactivation induced by reservoir production and injection, J. Petrol. Sci. Eng. , 60, 71– 85. https://doi.org/10.1016/j.petrol.2007.05.006 Google Scholar CrossRef Search ADS   Sone H., Zoback M., 2014. Time-dependent deformation of shale gas reservoir rocks and its long-term effect on the in situ state of stress, Int. J. Rock Mech. Min. Sci. , 69, 120– 132. Udias A., Madariaga R., Buforn E., 2014. Source Mechanisms of Earthquakes: Theory and Practice , Cambridge Univ. Press. Google Scholar CrossRef Search ADS   Van Eijs R.M.H.E., Mulders F., Nepveu M., Kenter C.J., Scheffers B.C. et al., 2006. Correlation between hydrocarbon reservoir properties and induced seismicity in the Netherlands, Eng. Geol. , 84, 99– 111. https://doi.org/10.1016/j.enggeo.2006.01.002 Google Scholar CrossRef Search ADS   Van den Bogert P., 2015. Impact of various modelling options on the onset of fault slip and the fault slip response using 2-dimensional Finite Element modelling. Report No SR.15.11455, Nederlandse Aardolie Maatschappij, Assen. Available at: https://nam-feitenencijfers.data-app.nl/embed/component/?id=onderzoeksrapporten, last accessed 2 October 2017. Van Thienen-Visser K., Breunese J.N., 2015. Induced seismicity of the Groningen gas field: history and recent developments, Leading Edge , 34, 664– 671. https://doi.org/10.1190/tle34060664.1 Google Scholar CrossRef Search ADS   Van Wees J.D., Mijnlieff H., Lutgert J., Breunese J., Bos C., Rosenkranz P., Neele F. et al., 2008. A Bayesian belief network approach for assessing the impact of exploration prospect interdependency: a application to predict gas discoveries in the Netherlands, AAPG Bull. , 92, 1315– 1336. https://doi.org/10.1306/06040808067 Google Scholar CrossRef Search ADS   Van Wees J.D., Buijze L., Van Thienen-Visser K., Nepveu M., Wassing B.B.T., Orlic B., Fokker P.A., 2014. Geomechanics response and induced seismicity during gas field depletion in the Netherlands, Geothermics , 52, 206– 219. https://doi.org/10.1016/j.geothermics.2014.05.004 Google Scholar CrossRef Search ADS   Wassing B.B.T., Van Wees J.D., Fokker P.A., 2014. Coupled continuum modeling of fracture reactivation and induced seismicity during enhanced geothermal operations, Geothermics , 52, 153– 164. https://doi.org/10.1016/j.geothermics.2014.05.001 Google Scholar CrossRef Search ADS   Zoback M., 2010. Reservoir Geomechanics , Cambridge Univ. Press. © The Author(s) 2017. Published by Oxford University Press on behalf of The Royal Astronomical Society. http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png Geophysical Journal International Oxford University Press

# Reservoir creep and induced seismicity: inferences from geomechanical modeling of gas depletion in the Groningen field

, Volume 212 (3) – Mar 1, 2018
11 pages

/lp/ou_press/reservoir-creep-and-induced-seismicity-inferences-from-geomechanical-0gVLk0FVla
Publisher
Oxford University Press
ISSN
0956-540X
eISSN
1365-246X
D.O.I.
10.1093/gji/ggx452
Publisher site
See Article on Publisher Site

### Abstract

Summary The Groningen gas field in the Netherlands experienced an immediate reduction in seismic events in the year following a massive cut in production. This reduction is inconsistent with existing models of seismicity predictions adopting compaction strains as proxy, since reservoir creep would then result in a more gradual reduction of seismic events after a production stop. We argue that the discontinuity in seismic response relates to a physical discontinuity in stress loading rate on faults upon the arrest of pressure change. The stresses originate from a combination of the direct poroelastic effect through the pressure changes and the delayed effect of ongoing compaction after cessation of reservoir production. Both mechanisms need to be taken into account. To this end, we employed finite-element models in a workflow that couples Kelvin-Chain reservoir creep with a semi-analytical approach for the solution of slip and seismic moment from the predicted stress change. For ratios of final creep and elastic compaction up to 5, the model predicts that the cumulative seismic moment evolution after a production stop is subject to a very moderate increase, 2–10 times less than the values predicted by the alternative approaches using reservoir compaction strain as proxy. This is in agreement with the low seismicity in the central area of the Groningen field immediately after reduction in production. The geomechanical model findings support scope for mitigating induced seismicity through adjusting rates of pressure change by cutting down production. Geomechanics, Induced seismicity 1 INTRODUCTION Time-dependent compaction of hydrocarbon reservoirs is a well-known phenomenon. Laboratory experiments indicate that time-dependent compaction can be of the same order of magnitude as direct (elastic) compaction for clastic and shale reservoirs (De Waal 1986; Hettema et al.2000; Sone & Zoback 2014; Pruiksma et al.2015). Results from subsidence inversion demonstrate that reservoir compaction rates are delayed relative to pressure drop, in line with laboratory data (Hettema et al.2000; Mossop 2012; Fokker & Van Thienen-Visser 2016). Compaction of clastic reservoir rocks is predicted to continue for decades after a potential shut-in of production or end of field life, and creep may constitute a significant portion (up to 10 per cent) of the compaction prior to cessation of production (Mossop 2012; Pruiksma et al.2015). In addition, time-dependent strain of overburden and underburden rock can contribute to time-dependent compaction and subsidence (Orlic & Wassing 2013; Chang et al.2014; Marketos et al.2015). Bourne et al. (2014) introduced a predictive model for seismicity due to fluid extraction from hydrocarbon reservoirs, using compaction strains as a proxy for potential stress change and seismic moment, in analogy to the strain–stress seismic moment proxy introduced by Kostrov (1974) and used for induced seismicity by McGarr (2014). For time-dependent compaction following a production stop, this model implies that seismicity rates remain significant and reduces slowly. This outcome was unexpected to us, as it does not agree with the most recent observations in the Groningen field in the Netherlands: a significant cut in production rates in the central area that was applied in 2014 was immediately followed by a statistically significant reduction in seismic event rates in the next year (Nepveu et al. 2016). In this paper, we demonstrate that the unexpected outcome of the proxy model can be attributed to a break which arises in the compaction strain–stress relationship, and which is related to the cut in production. The fundamental difference between our study and that of Bourne et al. (2014) is that we consider the stresses acting on the activated faults as the driving parameters of induced seismicity, rather than using the reservoir compaction as a proxy (Fig. 1). These stresses originate from a combination of the direct poroelastic effect through the pressure changes and from the delayed effect of ongoing compaction after cessation of reservoir production. Both mechanisms need to be taken into account, and this can explain the observed break. To demonstrate this, a simplified geomechanical model for the central area in the Groningen field is employed. Section 2 introduces the relevant data of the Groningen gas field and their significance in unraveling the relative control of pressure change and time-dependent compaction on seismicity. In Section 3, the geomechanical modeling approach is introduced. It is based on a viscoelastic (Kelvin Chain) finite-element solution for a two-dimensional reservoir section. The model evaluates the stress response of the reservoir, including time-dependent reservoir compaction, and highlights the reduction of stress build-up upon termination of production. In Section 4, the results in terms of seismic moment increase are evaluated using a novel analytical approach to predict seismic moment from the viscoelastic stress solution. The analytical technique is required as the finite-element code we use does not allow for combining Kelvin-Chain creep and Mohr–Coulomb failure. It is demonstrated that the model predicts a significant reduction of seismic moment increase after production stop, qualitatively in agreement with the observed reduction of event rates. We argue that the findings of this study are expected to hold for depleting reservoirs in general, assuming similar reservoir conditions. Figure 1. View largeDownload slide Schematic diagram of the geomechanical workflow in this paper, marked by blue lines. The model predicts a physical discontinuity in stress loading rate on faults upon the arrest of pressure change after cessation of reservoir production, with some delayed effect on stress change due to creep. The orange line represents the alternative model using compaction strain as proxy to predict seismic moment (Bourne et al.2014). The numbers refer to sections in the paper. Figure 1. View largeDownload slide Schematic diagram of the geomechanical workflow in this paper, marked by blue lines. The model predicts a physical discontinuity in stress loading rate on faults upon the arrest of pressure change after cessation of reservoir production, with some delayed effect on stress change due to creep. The orange line represents the alternative model using compaction strain as proxy to predict seismic moment (Bourne et al.2014). The numbers refer to sections in the paper. 2 GRONINGEN FIELD DATA The Groningen gas field is a very large gas field located onshore in the Netherlands. Gas production from the Groningen field started in 1963. Induced seismicity of the field first was recorded in 1991 (ML 2.4). During the subsequent 10 yr, induced seismicity stayed at a rate of about five events (ML ≥ 1.5) per year, increasing afterwards (e.g. van Thienen-Visser & Breunese 2015). In 2012, the largest event (ML 3.6) occurred, which caused the most damage to date. As a consequence, many studies have been carried out to look at the relationship between induced seismicity and gas depletion. These include predictive approaches based on compaction strain (Bourne et al.2014), and 2-D and 3-D finite-element-based mechanical studies (e.g. Sanz et al.2015; Van den Bogert 2015). These studies focused on the period prior to an 80 per cent production reduction, which was implemented on the 2014 January 17 to the production in the clusters in the central area of the Groningen field (Fig. 2). This area, constituting over 25 per cent of total field production prior to reduction, had shown the majority of induced events over the lifetime of the field (van Thienen-Visser & Breunese 2015). The seismic data of the central area of the field are shown in Fig. 2. The seismic event rate is marked by a statistically significant reduction of seismic events according to a recent study of Nepveu et al. (2016). This strongly suggests that the rate of seismic moment increase is considerably reduced after the production cut. Figure 2. View largeDownload slide Outline of the Groningen field (red line denotes gas water contact), coloured with pressure change predicted from NAM’s reservoir model 1.5 yr prior (left) and after (right) reduction in production. Induced events of the central area within the green polygon is shown by red dots, for the remainder of the field by grey dots (source www.knmi.nl), production clusters by triangles. The hydraulic diffusivity of the reservoir is such that pressure gets equilibrated in between production clusters in the central area within months to half a year, whereas pressure diffusion from the edges of the field to the centre would take over a few years (van Thienen-Visser & Breunese 2015). Figure 2. View largeDownload slide Outline of the Groningen field (red line denotes gas water contact), coloured with pressure change predicted from NAM’s reservoir model 1.5 yr prior (left) and after (right) reduction in production. Induced events of the central area within the green polygon is shown by red dots, for the remainder of the field by grey dots (source www.knmi.nl), production clusters by triangles. The hydraulic diffusivity of the reservoir is such that pressure gets equilibrated in between production clusters in the central area within months to half a year, whereas pressure diffusion from the edges of the field to the centre would take over a few years (van Thienen-Visser & Breunese 2015). The reduction in production in the central area of the Groningen field, the effects of which are closely monitored via the recorded seismicity, provides a unique data set documenting the break in trend of seismicity which can be attributed to relative change in rates of time-dependent compaction and pressure change in the reservoir. Past studies and associated data sets on gas depletion in other fields do not provide such data for the following reasons: A mid-field-life production cut at the scale applied in the central area of the Groningen field, is exceptional. Induced seismicity is not always monitored by a dedicated local network prior to the occurrence of induced events leading to stop of production, such as in the case of the Rotenburg Field in Germany (Dahm et al.2007). Therefore, seismic catalogues generally contain an insufficient number of events for statistical analysis. Gas fields are typically marked by a gradual reduction in production and the associated pressure change decays prior to production stop (e.g. Van Wees et al.2008). This appears to be consistent with a gradual reduction of seismic event rates occurring in the tail-end of production, as observed in the Strachan field in Canada (Baranova et al.1999). In fields where production has stopped, seismicity appears to return to background rates, as observed in the Strachan Field and the Roswinkel and Bergermeer fields in the Netherlands, with 25 and 4 seismic events (M > 1.5) during production and no events for 9 and 8 yr after production stop, respectively, in accordance with absence of natural seismicity (Van Wees et al. 2014). In the aforementioned fields, there was considerable time to allow a gradual cessation of time-dependent compaction prior to production stop, because of the decay in pressure rate changes. This renders these cases not suitable for the analysis pursued here. A strong aquifer support can result in equilibration of reservoir pressure on short timescales and the abundance of natural seismicity can render it difficult to analyse post-production evolution in view of induced seismicity. In summary, the Groningen field is marked by a sharp and prolonged drop in production by 80 per cent and pressure change rates in the central area (Fig. 2), of which the effects in terms of a sharp reduction in seismic event rates are well documented in the seismic catalogue (Nepveu et al. 2016). The field is located in a tectonically stable area, and characterized by limited aquifer pressure support. Therefore, it qualifies as an excellent case study to disentangle the pressure- and time-dependent compaction control on induced seismicity. 3 GEOMECHANICAL MODEL Semi-analytical and finite-element-based geomechanical models demonstrate that a progressive increase in induced events relates to progressive stress perturbations due to depletion, and that the relative magnitude of stress perturbations is strongly dependent on initial stresses, pressure change, reservoir geometry (Segall 1989; Roest & Kuilman 1994; Mulders 2003), reservoir and overburden creep (Orlic & Wassing 2013). In finite-element geomechanical models predicting displacement on faults, the expected evolution of cumulative seismic moment evolution (M0) can be determined via (e.g. Aki & Richards 2002):   $$M0 = \int G\ u\ dS\$$ (1)where S is surface area of faults, u is displacement and G is the shear modulus. The onset of slip is determined by the Coulomb failure function (CFF > 0) (e.g. Zoback 2010):   $${\rm CFF}\, = {{\sigma}_s}\, - \, \mu {{\sigma}_n}\!\!^{'}$$ (2)where σs is shear stress, σn' effective normal stress on the fault and μ the friction coefficient. Slip causes the fault to return to the yield condition (CFF = 0). This can include slip weakening effects and/or rate and state friction effects during rupture events (Ruina 1983; Dieterich 1994; Rutqvist et al.2013; Wassing et al.2014). Part of the slip may occur aseismically. In this simplified approach, the dynamic effects of slip and slip weakening are discarded, and it is assumed that all incremental slip is released seismically. In elastic models, the change in the CFF (ΔCFF) scales linearly to pressure change ΔP through changes in the total stress tensor on the fault (Mulders 2003; Soltanzadeh & Hawkes 2008):   $$\Delta \ {\sigma _{ij}} = \ \alpha ({\gamma _{ij}} - {\delta _{ij}})\ \Delta P$$ (3)where σij are the total stress tensor components, and γij normalized stress path parameters, describing the stress path in terms of change in stress tensor as a function of pressure change. α corresponds to Biot's coefficient. The stress path parameters are dependent on poroelastic parameters, fault geometry and offset, and can be determined semi-analytically from linear poroelasticity or numerically in finite-element models (Mulders 2003; Soltanzadeh & Hawkes 2008). They vary over the fault surface, in case of fault offset, but are static through time under spatially uniform depletion. Prior to the production cut, the Groningen field production strategy targeted at maintaining a uniform pressure change over the field. The central area of the field is marked by high permeability capable to equilibrate pressures in between well clusters in months, rather than years (van Thienen-Visser & Breunese 2015). For this reason, we ignore effects of spatial changes in pressure depletion which can cause stress perturbations (e.g. Altmann et al.2010). In the Groningen field, the pressure change has been decreasing almost linearly during the last 30 yr, prior to the cut in production (e.g. Bourne et al.2014; Sanz et al.2015). The associated stress paths are well capable of predicting an accelerating growth of seismic moment as a function of pressure change (Van Wees et al.2014). Mulders (2003) and Van Eijs et al. (2006) demonstrate that the ΔCFF is amplified with increasing contrast in the elastic properties of reservoir and surrounding rock. 3.1 Adopting time-dependent compaction Insufficient consensus of the exact creep mechanism in the Groningen field is currently available to justify a particular choice for the time-dependent compaction constitutive law (Mossop 2012; Hol et al.2015; Pruiksma et al.2015). For this study, a Kelvin-Chain model has been adopted. This model is based on a specific form of the generalized creep formulation (cf. Bazant et al.1988; De Borst & Van den Boogaard 1994):   $${\boldsymbol{\varepsilon}}(t) = \mathop \int \nolimits_0^t J\left( {t,\tau } \right){\boldsymbol{C\dot{\sigma }}}\left( \tau \right)d\tau$$ (4)ε(t) is strain tensor at time t, $${\boldsymbol{\dot{\sigma }}}$$ is stress rate tensor, J(t, τ) is creep function; τ  is time of load application and the components of C are defined by (De Borst & Van den Boogaard 1994):   $${C_{ijkl}} = \frac{1}{2}\left( {1 + \nu } \right)( {{\delta _{ik}}{\delta _{jl}} + {\delta _{il}}{\delta _{jk}}} ) - \ \nu \ {\delta _{ij}}{\delta _{kl}}$$ (5) In the Kelvin-Chain model, a single creep element in series with the elastic response is adopted (Fig. 4). The creep function, J(t, τ) corresponds to:   $$J\left( {t,\tau } \right) = {E^{ - 1\ }}\ + {E_{1\ }}^{ - 1}\left[ {1 - {e^{ - \left( {\frac{{t - \tau }}{{{\lambda _1}}}} \right)}}} \right]$$ (6) The spring with stiffness E allows for an elastic stress response, which should be in accordance to the static Young's modulus. Lab measurements document E typically in the range of 8–25 GPa for rock samples which are representative for the Groningen gas field (Hol et al.2015; Pruiksma et al.2015). The Kelvin-Chain component combining the spring with constant E1 and the dashpot with Newtonian viscosity η1 represents the time-dependent compaction. For uniaxial compaction, the creep component εcz(t) is constructed by:   $$\varepsilon _z^c(t) = \int_{0}^{t}{{\left[ {1 - {e^{ - \left( {\frac{{t - \tau }}{{{\lambda _1}}}} \right)}}} \right]}}{C_{{\rm{tinf}}}}\dot{\Delta }P\left( \tau \right)d\tau$$ (7)where $${C_{{\rm{tinf}}}} = \ \frac{{( {1 - v - 2{v^2}} )/( {1 - v} )}}{{{E_{1\ }}}}$$. It follows that the uniaxial Kelvin-Chain model corresponds to the decay model for time-dependent compaction introduced by Mossop (2012):   $$\ {\varepsilon ^c}\!\!_z\ (t) = f'\ *g\ = \ {{\cal F}^{ - 1}}\left[ {{\cal F}\left[ {f'\ } \right]\ {\cal F}\left[ g \right]} \right]$$ (8)with   \begin{equation*}f'(t) = \frac{d}{{dt}}\ \left[ {1 - {e^{ - \left( {\frac{{t - \tau }}{{{\lambda _1}}}} \right)}}} \right]\ = \frac{1}{{{\lambda _1}}}{e^{ - t/{\lambda _1}}}\end{equation*}  \begin{equation*}g\ (t) = {C_{{\rm{tinf}}}}{\rm{\ \Delta }}P(t)\end{equation*} where $${\cal F}$$ denotes the Fourier Transform in the time domain and $${{\cal F}^{ - 1}}$$ its inverse, ΔP represents pressure change and $${\lambda _1} = \frac{{{\eta _1}}}{{{E_1}}}$$ is the relaxation time of time-dependent compaction. The Kelvin-Chain model cannot account for irreversible deformation, and stress paths in alternative constitutive models will deviate from this paper's model results. However, the relative change in stress response upon cessation of depletion is not likely to change significantly, as this is primarily driven by the break in loading mechanism after an arrest of pressure changes. Creep scenarios: in order to evaluate the effects of time-dependent compaction, three different scenarios of E/E1 ratios of final time-dependent (creep) strain and elastic strain are considered in the Kelvin-Chain model (Table 2). The model parameters have been chosen such that the total final strain in all models is the same. The models with E/E1 = 0.5 and 1 time-dependent compaction closely match laboratory experiments (e.g. Hol et al.2015; Pruiksma et al.2015), whereas the E/E1 = 5 case mimics the compaction model of Mossop, where time-dependent compaction is more dominant than elastic deformation. The adopted models predict relatively large creep in the first years after a sudden halt in pressure change compared to rate type and linear compaction models (e.g. Pruiksma et al.2015), for the same amount of final creep. In this sense, the Kelvin-Chain model can be considered most pessimistic in terms of reservoir strains and stress changes associated with these strains. Pseudo-elastic approach: under uniaxial loading conditions, the total strain for the Kelvin Chain is:   $${\varepsilon _z}(t)\ = \ {\varepsilon ^c}\!\!_z(t) + {\varepsilon ^e}\!\!_z(t)$$ (9)with direct elastic strain εez(t) corresponding to the instantaneous reacting spring in the Kelvin Chain:   $${\varepsilon ^e}\!\!_z(t) = {C_m}\,\,\Delta P(t),\quad{C_m} = \frac{{(1 - v - 2{v^2})/(1 - v)}}{E})$$ (10) The resulting stress–strain relationship can be written as:   $${\varepsilon _z}(t) = \widetilde {{C_m}}{\rm{\ }}\left( {\rm{t}} \right){\rm{\ \Delta }}P(t)$$ (11)with   \begin{equation*}\widetilde {{C_m}} = {\varepsilon _{{\rm{ratio}}}}(t)\ {C_m}\end{equation*}  \begin{equation*}{\varepsilon _{{\rm{ratio}}}}(t) = \frac{{{\varepsilon _z}(t)}}{{{\varepsilon ^e}\!\!_z(t)}}\end{equation*} At infinite time, εratio(t) corresponds to E/E1. The scaling of the stress–strain tensor relationship is scalar in the creep function, invariant for v. Therefore, in the stress–strain tensor relationship in a finite-element formulation, εratio(t) can be used to obtain a time-dependent pseudo-compaction coefficient or corresponding pseudo-stiffness for the reservoir zone (Young's modulus referred to as Ep(t)). This allows us to relate relative stress responses as a function of time-dependent compaction from a pseudo-elastic analysis adopting the pseudo-stiffness for reservoir and E for surrounding rock. A benefit of this pseudo-elastic approach is it allows numerically to be coupled in finite-element models with Mohr–Coulomb failure, whereas for the Kelvin-Chain model this is not supported in the workflow we use. 3.2 Model geometry and parametrization A plane-strain elastic finite-element model of a depleting gas reservoir is constructed and used to analyse the sensitivity of total stress scaling and ΔCFF to different Kelvin-Chain parameters. For the finite-element modeling, DIANA FEA (http://dianafea.com) is used, a commercially available finite-element code that is widely used for geomechanical analysis in oil and gas reservoirs, including the analysis of Mohr–Coulomb failure and slip. More details of the finite-element modeling procedure are given in Orlic & Wassing (2013). The numerical model incorporates two reservoir compartments, separated by a fault with a throw of half the reservoir thickness (Fig. 5 and Table 1), which can be considered representative for major fault zones in the central area of the Groningen field, where most seismic events have occurred (Fig. 2; Bourne et al.2014). The depth and thickness range, and underlying parametrization is in accordance with earlier 2-D models for the Groningen gas field (e.g. Van den Bogert 2015). The plane strain finite-element mesh models a 2-D section of 6 km width and 6 km depth. The model consists of roughly 135 000 nodes which comprise approximately 44 000 quadratic quadrilateral plane-strain elements representing the bulk rock and ∼250 quadratic quadrilateral interface elements representing the fault. The mesh is finest around the fault with a resolution of 1 m at reservoir depth and gets gradually coarser away from the fault and the reservoir (Fig. 3). The model resolution has been chosen sufficiently high to capture high bending stresses close to the fault and at the interface of reservoir and adjacent rock. Figure 3. View largeDownload slide Geometry and finite-element mesh for the model setup. Figure 3. View largeDownload slide Geometry and finite-element mesh for the model setup. Table 1. Parameters of the geomechanical model, based on earlier studies for the Groningen field (cf. Mulders 2003; Van Wees et al. 2014; Van den Bogert 2015). Parameter  Value  Unit  E  18  GPa  V  0.2  –  Top reservoir  2900  m  Thickness reservoir  150  m  Fluid density  1150  kg m−3  Rock density  2260  kg m−3  Friction angle  30  deg  Fault dip  70  deg  Depletion  25  MPa  K0eff  0.45  –  Parameter  Value  Unit  E  18  GPa  V  0.2  –  Top reservoir  2900  m  Thickness reservoir  150  m  Fluid density  1150  kg m−3  Rock density  2260  kg m−3  Friction angle  30  deg  Fault dip  70  deg  Depletion  25  MPa  K0eff  0.45  –  View Large The plane-strain elements evaluate stresses and pressures on eight integration points per element, while the interface elements contain six integration points per element. The model is laterally constrained (zero displacement in horizontal direction) on the sides and vertically constrained (zero displacement in vertical direction) at the bottom, and subject to gravity loading. This effectively leads to a stress response on the sides of the model which corresponds to an infinite lateral extension of the reservoir. For each loading step, the systems of equations are solved using an iterative solution method which applies incomplete LU-decomposition preconditioning. For the Kelvin-Chain creep model, the elements representing the reservoir are governed by a Kelvin-Chain constitutive model, while the surrounding rock is modeled to be fully linearly elastic. In order to prevent the development of shear stresses and shear strains in the model due to the viscoelastic response of the reservoir to the initialization load, the model is initialized with ρrock = ρfluid. This results in a model with increasing pore pressure with depth, but without effective stresses. The reservoir is then depleted, resulting in non-zero effective stresses. The stress changes of this model are superimposed on the results of a fully elastic model which represents the initial stress state. For the initial stress state, the reservoir is assumed overpressured relative to hydrostatic conditions (Bourne et al.2014). The initial ratio of effective horizontal-to-vertical stress is set to K0eff = 0.45 (cf. Van Wees et al.2014). In the pseudo-elastic approach, we adopt the pseudo-stiffness Ep(t) for the reservoir, resulting in contrasting stiffness of the reservoir and the surrounding rock. The models are initialized in two steps to prevent development of shear stresses and shear strains in the pre-depletion state. First, the model is initialized with a homogeneous stiffness all throughout the model and a prescribed K0eff. The resulting stress state is then transferred to a model which contains the stiffness contrast. The final depletion in the reservoir is 25 MPa, approximately in agreement with pressure depletion at present day in the field. Depletion is implemented to occur simultaneously, at the same rate, in both reservoir compartments in yearly load steps. The pressure depletion of 25 MPa is built-up linearly in a time period of 60 yr, followed by 40 yr of no pressure change. The decrease in pore pressure is the same for all values of the horizontal coordinate, in agreement with high diffusion rates observed in the Groningen gas field (van Thienen-Visser & Breunese 2015), marked by a delay time of a couple of months to equilibrate pressures from well cluster in the central area The linear pressure depletion rate has been constrained by the observed linear pressure reduction in the last 30 yr prior to reduction in production (Bourne et al.2014). The time period of pressure depletion has been extended to 60 yr relative to approximately 50 yr of actual production in the field, to compensate for higher pressure rate changes in the first 20 yr of production (Bourne et al.2014). The fault zone is assumed to be subject to the reservoir fluid pressure when reservoir rock is present on both sides of the fault. For the portion of the fault where reservoir rock and surrounding rock are juxtaposed, two scenarios for pressure in the fault are considered: actual reservoir pressure and ambient (initial reservoir) pressure. Obviously, a single fault geometry cannot capture the complex response expected from many mapped faults in the area with variable throw, dip and strike (Fig. 2; Sanz et al.2015). Uncertainty in the location of seismic events hampers an allocation of seismic events to particular faults. Fault strike, dip, throw, natural stress variability and variability in mechanical properties, can all have a significant control on in situ stress and spatial and temporal variability in seismicity (e.g. Van Wees et al.2014; Sanz et al.2015). However, these aspects are not expected to have a prime control on the relative change in seismic response upon production stop, as they do not affect strongly the ΔCFF as a function of pressure change. 3.3 Coulomb failure function (ΔCFF) evolution The evolution of ΔCFF is constructed from the resulting stress changes of the finite-element model for the three time-dependent compaction scenarios from Table 2, with creep to elastic deformation ratios (E/E1) of 0.5, 1 and 5 at infinite time. Fig. 4 shows an example of the adopted pressure change, and resulting effective normal stress and shear stress changes for E/E1 = 1 for both the creep (Kelvin Chain) and pseudo-elastic finite-element model. Fig. 4 also shows the ΔCFF for all three scenarios. During production, the ΔCFF increases almost linearly with production rate. The total ΔCFF increase differs, depending on the relative contribution of time-dependent compaction and is in close agreement with the stress response in the pseudo-elastic approach (eq. 8). Figure 4. View largeDownload slide Finite-element model solutions for Kelvin-Chain creep (solid lines) and pseudo-elastic models (dashed lines). Panels (a) and (b) for E/E1 = 1 (a) effective normal stress on the fault assuming reservoir pore pressure, (b) shear stress on the fault. (c)–(e) Change in Coulomb Failure Function for E/E1 = 0.5, 1 and 5, respectively. The sign of shear stress is positive if shear sense is normal faulting. The time in years relative to production stop is indicated by t. The depth extent of the reservoir zones is indicated on both sides of the fault in grey. Figure 4. View largeDownload slide Finite-element model solutions for Kelvin-Chain creep (solid lines) and pseudo-elastic models (dashed lines). Panels (a) and (b) for E/E1 = 1 (a) effective normal stress on the fault assuming reservoir pore pressure, (b) shear stress on the fault. (c)–(e) Change in Coulomb Failure Function for E/E1 = 0.5, 1 and 5, respectively. The sign of shear stress is positive if shear sense is normal faulting. The time in years relative to production stop is indicated by t. The depth extent of the reservoir zones is indicated on both sides of the fault in grey. Table 2. Viscoelastic parameters in the Kelvin-Chain model (Fig. 3). Scenario  E (GPa)  E1 (GPa)  λ1 = η1/E1 (yr)  50 per cent creep  13.5  27  7.3  100 per cent creep  18  18  7.3  500 per cent creep  54  10.8  7.3  Scenario  E (GPa)  E1 (GPa)  λ1 = η1/E1 (yr)  50 per cent creep  13.5  27  7.3  100 per cent creep  18  18  7.3  500 per cent creep  54  10.8  7.3  Notes: Viscosity η1 is chosen such that it agrees with a relaxation time λ1 = 7.3 yr (Mossop 2012). View Large In the 40 yr after production, ΔCFF increases only moderately by a few per cent relative to the ΔCFF build-up during production, and the curves are hard to discriminate in Fig. 4. The moderate amount of total ΔCFF change after production stop, and the abrupt change in ΔCFF evolution, are easily understood from the fundamentally different mechanisms contributing to stress changes. Prior to production stop, both the pressure change (eq. 3) and the progressive growth of εratio(t) is contributing to ΔCFF, whereas after production stop the only contribution to ΔCFF is the moderate increase of εratio(t) as a function of progressive compaction (cf. eq. 11). In Fig. 4, all the creep models show that ΔCFF has a distinct peak relative to the pseudo-elastic approach at the structural interfaces of the reservoir and surrounding rock. These are related to the creep causing higher and localized deviatoric stresses at the structural interfaces, during depletion. It highlights the limitations of the pseudo-elastic approach to reproduce the exact value of stress changes at structural interfaces. Fig. 5 shows the relative change in ΔCFF, after production stop. Here, only the creep model results are presented, as it was recognized that the pseudo-elastic model is not well capable of reproducing the relative ΔCFF changes of the Kelvin-Chain model, in the first years immediately after production stop. The change in ΔCFF after production stop is lower than the change in ΔCFF just prior to production stop for all cases considered—with the most significant difference experienced for the strongest time-dependent compaction scenarios. Figure 5. View largeDownload slide (a)–(c) Change in Coulomb failure function on the fault around productions for E/E1 = 0.5, 1 and 5, respectively, assuming reservoir pore pressure on the fault. The time in years relative to production stop is indicated by t. Figure 5. View largeDownload slide (a)–(c) Change in Coulomb failure function on the fault around productions for E/E1 = 0.5, 1 and 5, respectively, assuming reservoir pore pressure on the fault. The time in years relative to production stop is indicated by t. Figs 4 and 5 consistently show that ΔCFF is monotonously increasing before and after production stop for the portion of the faults, which is subject to failure. This indicates that the stress path remains similar. Furthermore as the total ΔCFF change after production stop is so small compared with the production period, it is expected that potential unsolicited effects of reversibility of creep, which is inherent of the Kelvin-Chain model, are not significant. Summarizing, for all creep scenarios considered, the predicted increase of ΔCFF after a production stop is few per cent relative to the ΔCFF which had been built-up during production, marked by a strong break in stress build-up rates, at the production cut. The scenarios also show a very quick decay in the build-up of ΔCFF in the first few years after the end of production. Furthermore, the models show a close resemblance of stress responses of the Kelvin-Chain model and the pseudo-elastic approach, except at structural interfaces and in the years immediately after production stop. The break in stressing rate is in agreement with the observed break in seismic event rates (Fig. 2; Nepveu et al.2016). 4 SEISMIC MOMENT The previous sections highlight the generic relationship between ΔCFF with variations in time-dependent compaction. The bearings of modeled ΔCFF changes on seismicity are evaluated in this section, in terms of the relative effect on seismic moment increase. As mentioned before, the finite-element model approach used here does not support a combined Kelvin-Chain and Mohr–Coulomb material model, for calculation of fault slip and moment (cf. eq. 1). To overcome this limitation, the relative growth of seismic moment is determined by a novel semi-analytical approach for the construction of the progressive increase of slip area and displacement of the fault, based on the elastic stress solution only. 4.1 Semi-analytical model to resolve moment from the elastic stress solution For a circular fault area, a closed relationship exists between area S, radius r, average slip u, stress drop Δσ and seismic moment M0 (Eshelby 1957; Keilis-Borok 1959; Udias et al.2014):   $$M0 = \frac{16}{7\pi} \Delta \sigma rS$$ (12) This equation is commonly used by seismologists to estimate Δσof an earthquake after determining M0 and source dimensions of the rupture surface (e.g. Udias et al.2014), and has been generalized for elliptic surfaces of fault slip (Madariaga 1979). In this case, the slip and moment are scaled by a correction factor which takes into account both the ellipticity and the slip direction on the fault. For plane-strain dip-slip faulting conditions, the correction factor cancels the term $$\frac{16}{7\pi }$$. Considering a slipping fault section with length l in strike and dip direction, and adopting S = l2, and $$r = \frac{l}{{\sqrt {\rm{\pi }}}}$$ for the radius of the equivalent circular fracture and substituting in eq. (9) gives the fault moment density [N] per m strike:   $$M{0_m} = \ \Delta \sigma \frac{{\ {l^2}}}{{\sqrt \pi }}$$ (13) In the approach adopted here, moment is determined from the determination of slip length and average stress drop from the elastic stress solution. Initially, l = l0 is estimated to correspond to the portion of the fault where CFF > 0, and CFF corresponding to Δσ. In reality, the stress drop would cause the fault length to grow progressively beyond the fault tips, with a stress redistribution resulting in CFF ≤ 0 on the slipping portion of the faults. We designed an algorithm for the estimation of the increased slip length based on the assumption that the elastic energy is redistributed by shear. To this end, the CFF is discretized in patches along the dip of the fault at a resolution of Δx = 1 m, and l and Δσ are solved by the following steps: Determine the range of patches [istart.iend] where CFFi > 0. For this range calculate the integral of excess Coulomb stress $$Y_{\rm{sum\ }} = \mathop \sum \nolimits_{{i_{{\rm{start}}}}}^{{i_{{\rm{end}}}}} {\rm{CF}}{{\rm{F}}_i}{\rm{\ \Delta \ }}x$$. Set Yistart = Yiend =  0.5 Ysum Beyond the tips of the fault distribute the excess stress Yistart and Yiend following a function of the distance di of the centre of each patch to the nearest tip of the fault (corresponding to istart or iend). Increment CFFi with ai =  Adi − 1.3. A is chosen such that the integral of the function equates the fraction of Y to be distributed beyond the tip under consideration. In the range of patches [istart.iend] set CFFi = 0. For incremented CFFi in the other patches, include patch i in the range of slipping patches if it results in failure (CFFi > 0) and recalculate the integral of excess Coulomb stress $$Y{\rm{\ }} = \mathop \sum \nolimits_{{i_{{\rm{start}}}}}^{{i_{{\rm{end}}}}} {\rm{CF}}{{\rm{F}}_i}{\rm{\ \Delta \ }}x$$. Divide Y over Yistart and Yiend according the portion originated from each of the tips of the fault. If Y > 0, repeat steps 2–4. Set $$l{\rm{\ }} = {\rm{\ }}( {{i_{{\rm{end}}}} - {i_{{\rm{start}}}}} ){\rm{\Delta }}x,\Delta \sigma {\rm{\ }} = \frac{{Y{\rm{sum}}}}{l}{\rm{\ }}$$. In the algorithm, a key assumption is that the area integral of Coulomb stress change remains unchanged upon slip (cf. Okada 1992; Wassing et al.2014). The characteristics of the power function for the Coulomb stress redistribution in step 2 have been determined from a numerical analysis of stress response for l in the range of 10–100 s of metres, for a steeply dipping fault at 3000 m depth, infinite strike length and adopting the analytical solution of slip for a constant stress drop (e.g. Madariaga 1979). The stress response has been calculated in a regular grid with 10 m spacing from an analytical solution of stress for rectangular sources of slip following the grid resolution (cf. Madariaga 1979; Okada 1992). The average input slip and output stress response of each patch of the grid is determined using two-point Gauss–Legendre Integration, which gave a numerically convergent result relative to a finer grid resolution for stress distribution. The stress redistribution is not sensitive to elastic parameters, nor to the exact orientation of the slip plane in the stress field, as can be derived from the stress equations (Okada 1992). Fig. 6 depicts the magnitude of the stress increase decaying away from the fault tip. It demonstrates that the redistribution of stress affects a considerable area, and can be well represented by a power function of distance from the tip for the decay. Figure 6. View largeDownload slide Predicted Coulomb stress change for analytical stress response (cf. Okada 1992, friction angle cf. Table 1), as a function of distance from the fault tip, for an elliptical displacement along a normal fault oriented favourably in the stress field and with a maximum of 0.01 m displacement for l = 30, 110 and 210 m. The dashed line represents the best fit x−1.3 line for l = 110 m (cut-off at x = 300 m). Figure 6. View largeDownload slide Predicted Coulomb stress change for analytical stress response (cf. Okada 1992, friction angle cf. Table 1), as a function of distance from the fault tip, for an elliptical displacement along a normal fault oriented favourably in the stress field and with a maximum of 0.01 m displacement for l = 30, 110 and 210 m. The dashed line represents the best fit x−1.3 line for l = 110 m (cut-off at x = 300 m). In order to test the robustness of the method for estimating seismic moment, various benchmarks comparing finite-element models resolving slip and slip length along faults incorporating Mohr–Coulomb failure (cf. Orlic & Wassing 2013) have been run, applying eq. (1) and comparing these to the results of the semi-analytical approach. For the model setup, the initial geometry, properties, boundary and initial conditions used for the previous models are applied. Static elastic properties through time have been adopted, with offset of the fault in the range (0.5–1)-D, and both ambient and reservoir pressure in the faults have been considered. Two model results for reservoir pressure are shown for 0.5-D and 1-D displacement in Fig. 7. This figure shows an excellent correspondence of fault moment density obtained by the finite-element method and the analytical method. In all model runs, the moment predicted from the semi-analytical model runs show less than 5 per cent deviation from the moment in the finite-element solutions, whereas the uncorrected l = l0 would significantly underestimate the moment by 25–70 per cent. From this, it is concluded that the semi-analytical approach provides a good estimator of moment growth from an elastic stress solution. Figure 7. View largeDownload slide Comparison of moment predicted by finite-element analyses incorporating Mohr–Coulomb failure and the semi-analytical model, for 0.5-D and 1-D offsets prior to production stop, adopting uniform elastic properties (cf. Table 1, reservoir pressure in fault). Figure 7. View largeDownload slide Comparison of moment predicted by finite-element analyses incorporating Mohr–Coulomb failure and the semi-analytical model, for 0.5-D and 1-D offsets prior to production stop, adopting uniform elastic properties (cf. Table 1, reservoir pressure in fault). The strength of the semi-analytical approach is that it provides a direct relationship between relative changes in ΔCFF and moment evolution, which is not biased to potential changes in response due to variability of in situ stress, specific orientation of faults and contrasts in elastic and frictional properties. Furthermore, it avoids potentially ambiguous results of finite-element slip with respect to the shear modulus to be used in eq. (1) to construct moment, if the model is marked by contrasting material behaviour on both sides of the fault. 4.2 Post-production time-dependent compaction and seismic moment increase In Fig. 8, the evolution of the Kelvin-Chain compaction and fault seismic moment density determined from the Coulomb failure function by the algorithm in Section 4.1 is plotted, normalized to its value at the end of production. Furthermore, the moment density has been constructed in accordance to the heuristic correlation of compaction strain and the logarithm of seismic moment increase (cf. Bourne et al.2014, median in their fig. 18a). These have been constructed for a reservoir segment of 2 km wide in map view, representative for average fault spacing in the central area. The Kelvin-Chain compaction in Fig. 8 for the E/E1 = 0.5 and 1 scenarios shows a discontinuity after production stop. This is supported by a clear discontinuity in GPS subsidence data in the central area of the Groningen field occurring nine weeks after production stop (Pijpers & Van der Laan 2016). Markedly, the ratio of seismic moment increase after and before production stop in the geomechanical models is 2–10 times lower and more than 10 times lower than those predicted from the approach using compaction as proxy (Bourne et al.2014), for the first few years and decades after shut-in, respectively (Fig. 8). This holds for all three time-dependent compaction scenarios that are considered here, as well as for both the reservoir and ambient pressure conditions on the fault. For a model that does not incorporate creep, the normalized stress changes and increase in seismic moment would be very similar to the E/E1 = 0.5 scenario and would arrest when production stops. The geomechanical predictions of all creep scenarios are in accordance with the strong reduction of induced events after 80 per cent cut in production rates in the central area in the Groningen field (cf. Fig. 2, Nepveu et al.2016). For reservoir pressure condition on the fault, the predicted rate of fault moment density increase for the E/E1 = 0.5 and 1 scenarios is about 0.5 GN yr−1, 1.5 yr prior to production stop. The central area is marked by about 5 × 1013 N m moment increase per year in the same period based on the seismic events. In the geomechanical model, 100 km of fault length would suffice to generate the observed moment rate. For the ambient pressure scenario less than 5 km of fault length would be sufficient. Cumulative fault length in the central area is in the order of 50 km (Fig. 2). The evolution of seismic moment prior to production stop is showing less accelerated growth as observed in the field (Bourne et al.2014). This may be explained by the fact that the 2-D model fails to capture important effects contributing to increased growth, including spatial variability in fault strike, dip, reservoir offset and thickness and initial stress. These can jointly contribute to a progressively increasing fault length involved in slip and would result in a more convex pattern of seismic moment growth. Figure 8. View largeDownload slide Left vertical axis: compaction of the reservoir for different creep models. On the right vertical axis: seismic moment density normalized to the seismic moment at t = 0 using the Bourne et al. (2014) heuristic correlation of compaction and seismic moment (α in the legend refers to the moment partitioning coefficient at t = 0, see Bourne et al. (2014) and the text for details) and using the analytical method described in this paper (moment density in the legend refers to the moment density at t = 0). Figure 8. View largeDownload slide Left vertical axis: compaction of the reservoir for different creep models. On the right vertical axis: seismic moment density normalized to the seismic moment at t = 0 using the Bourne et al. (2014) heuristic correlation of compaction and seismic moment (α in the legend refers to the moment partitioning coefficient at t = 0, see Bourne et al. (2014) and the text for details) and using the analytical method described in this paper (moment density in the legend refers to the moment density at t = 0). 5 DISCUSSION AND CONCLUSIONS The Kelvin Chain creep and the geometry of the model do not allow an in-depth assessment of the sensitivities of the results to variability in geometrical aspects of the reservoir and alternative creep behaviour. However, the prediction of a significant reduction of seismicity after cessation of production is expected to hold as well for alternative geometries and constitutive laws, as the first-order effects are determined by the arrest in pressure change (eq. 3). In the model developed, it is assumed that pressure does not change after production stop. Because of continued production in the other production clusters of the Groningen field (Fig. 2), in two years it is expected that reservoir pressures will be further reduced in the central area, due to hydraulic diffusion (van Thienen-Visser & Breunese 2015), which is likely to increase the likelihood of seismicity, compared to the present state. The prediction of the model, building from the relative increase in expected moment, needs to be interpreted with care in view of earthquake processes. A critically stressed fault at the end of the production period, which is incrementally loaded by an additional few per cent of CFF predicted by this model after stop of production, could be sufficient to trigger a relatively large event (i.e. M > 3). Therefore, it is important to note that the model does not preclude the possibility of a large event. Concluding, our geomechanical stress-based models incorporating reservoir creep, predict that seismic moment evolution after a production stop is subject to a very moderate increase, in the order of 2–10 times and more than 10 times less than the values predicted by the approach correlating seismic moment to progressive compaction strain (Bourne et al.2014), for the first years and decades following cessation of production, respectively. This difference is easily understood by the fact that the strain-based model of Bourne et al. (2014) ignores stress rate changes, which occur at a production stop. Our stress-based model is in agreement with the low seismicity encountered in the central area of the Groningen field, where production has been reduced by 80 per cent, on 2014 January 17. The model findings support scope for an immediate mitigation of induced seismicity through adjusting production scenarios, as the seismicity rates are expected to be correlated primarily with rates of pressure change. The findings of this study may well be generalized to other depleting reservoirs in similar structural settings, as the strong jump of stress loading behaviour on faults as a function of arrest of pressure changes (eq. 3) is universal, and prolonged stress effects of time-dependent compaction are typically one order of magnitude lower than those of pressure change. A number of issues were beyond the scope of this study but need to be addressed to augment to the robustness of our findings. These include: more realistic constitutive laws for reservoir creep; adoption of the full 3-D complexity of reservoir structure and of pressure depletion; pressure diffusion and creep in the over and underburden of the reservoir and effects of the rupture process in terms of a rate and state friction model. ACKNOWLEDGEMENTS The data used in this manuscript are publically available: seismicity data from the Royal Netherlands Meteorological Institute KNMI (www.knmi.nl); monthly production data from production wells in the Groningen field, available from the Geological Survey of the Netherlands TNO (www.nlog.nl) and numerous data and reports on the Groningen gas field (namplatform.nl). The authors would like to thank Jaap Breunese, Dries Hegen, Bogdan Orlic and Brecht Wassing for constructive comments on earlier versions of the manuscript. The manuscript benefited considerably from constructive comments from Alexander Droujinine and an anonymous reviewer REFERENCES Aki K., Richards P., 2002. Quantitative Seismology , University Science Books. Altmann J.B., Müller T.M., Müller B.I.R., Tingay M.R.P., Heidbach O., 2010. Poroelastic contribution to the reservoir stress path, Int. J. Rock Mech. Min. Sci. , 47( 7), 1104– 1113. https://doi.org/10.1016/j.ijrmms.2010.08.001 Google Scholar CrossRef Search ADS   Baranova V., Mustaqeem A., Bell S., 1999. A model for induced seismicity caused by hydrocarbon production in the Western Canada Sedimentary Basin, Can. J. Earth Sci. , 36( 1), 47– 64. https://doi.org/10.1139/e98-080 Google Scholar CrossRef Search ADS   Bazant Z.P., Dougill J., Huet C., Tsubaki T., Wittmann F.H., 1988. Material models for structural creep analysis, in Mathematical Modelling of Creep and Shrinkage of Concrete, ed. Bazant Z.P.J., Wiley & Sons, Chichester, England, pp. 275– 310. Bourne S., 2014. A seismological model for earthquakes induced by fluid extraction from a subsurface reservoir, J. geophys. Res. , 119, 8991– 9015. Google Scholar CrossRef Search ADS   Chang C., Mallman E., Zoback M., 2014. Time-dependent subsidence associated with drainage-induced compaction in Gulf of Mexico shales bounding a severely depleted gas reservoir, AAPG Bull. , 98( 6), 1145– 1159. https://doi.org/10.1306/11111313009 Google Scholar CrossRef Search ADS   Dahm T., Krüger F., Stammler K., Klinge K., Kind R., Wylegalla K., Grasso J.-R., 2007. The 2004 Mw 4.4 Rotenburg, Northern Germany, earthquake and its possible relationship with gas recovery, Bull. seism. Soc. Am. , 97( 3), 691– 704. https://doi.org/10.1785/0120050149 Google Scholar CrossRef Search ADS   De Borst R., Van den Boogaard A.H., 1994. Finite-element modelling of deformation and cracking in early-age concrete, J. Engin. Mech. , 120( 12), 2519– 2534. Google Scholar CrossRef Search ADS   De Waal J., 1986. On the rate type compaction behaviour of sandstone reservoir rock, PhD thesis , Technical University Delft, Delft, Netherlands. Dieterich J., 1994. A constitutive law for rate of earthquake production and its application to earthquake clustering, J. geophys. Res. , 99( B2), 2601– 2618. https://doi.org/10.1029/93JB02581 Google Scholar CrossRef Search ADS   Eshelby J.D., 1957. The determination of the elastic field of an ellipsoidal inclusion, and related problems, Proc. R. Soc. A: Math. Phys. Eng. Sci. , 241( 1226), 376– 396. https://doi.org/10.1098/rspa.1957.0133 Google Scholar CrossRef Search ADS   Fokker P.A., Van Thienen-Visser K., 2016. Inversion of double-difference measurements from optical leveling for the Groningen gas field, Int. J. Appl. Earth Observ. Geoinform. , 49, 1– 9. https://doi.org/10.1016/j.jag.2016.01.004 Google Scholar CrossRef Search ADS   Hettema M.H.H., Schutjens P.M.T.M., Verboom B.J.M., Gussinklo H.J. et al.  , 2000. Production-induced compaction of a sandstone reservoir: the strong influence of stress path, SPE Reserv. Eval. Eng. , 3( 04), 342– 347. https://doi.org/10.2118/65410-PA Google Scholar CrossRef Search ADS   Hol S., 2015. Long-term compaction behavior of Permian sandstones - an investigation into the mechanisms of subsidence in the Dutch Wadden Sea, in 49th US Rock Mechanics/Geomechanics Symposium , American Rock Mechanics Association, San Francisco, CA. Keilis-Borok V., 1959. On estimation of the displacement in an earthquake source and of source dimensions, Ann. Geofys. , 12, 205– 214. Kostrov V., 1974. Seismic moment and energy of earthquakes, and seismic flow of rock, Phys. Solid Earth , 1, 13– 21. Madariaga R., 1979. On the relation between seismic moment and stress drop in the presence of stress and strength heterogeneity, J. geophys. Res. , 84( B5), 2243– 2250. https://doi.org/10.1029/JB084iB05p02243 Google Scholar CrossRef Search ADS   Marketos G., Govers R., Spiers C.J., 2015. Ground motions induced by a producing hydrocarbon reservoir that is overlain by a viscoelastic rocksalt layer: a numerical model, Geophys. J. Int. , 203( 1), 198– 212. https://doi.org/10.1093/gji/ggv294 Google Scholar CrossRef Search ADS   McGarr A., 2014. Maximum magnitude earthquakes induced by fluid injection, J. geophys. Res. , 119, 1008– 1019. Google Scholar CrossRef Search ADS   Mossop A., 2012. An explanation for anomalous time dependent subsidence, in 46th US Rock Mechanics/Geomechanics Symposium , American Rock Mechanics Association, Chicago, IL. Mulders F., 2003. Modelling of stress development and fault slip in and around a producing gas reservoir, PhD thesis , Technical University of Delft, Delft, The Netherlands. Nepveu M., Van Thienen-Visser K., Sijacic D., 2016. Statistics of seismic events at the Groningen field, Bull. Earthq. Eng. , 14( 12), 3343– 3362. Google Scholar CrossRef Search ADS   Okada Y., 1992. Internal deformation due to shear and tensile faults in a half-space, Bull. seism. Soc. Am. , 82, 1018– 1040. Orlic B., Wassing B.B.T., 2013. A study of stress change and fault slip in producing gas reservoirs overlain by elastic and viscoelastic caprocks, Rock Mech. Rock Eng. , 46, 421– 435. https://doi.org/10.1007/s00603-012-0347-6 Google Scholar CrossRef Search ADS   Pijpers F., Van der laan F., 2016. Trend changes in ground subsidence in Groningen. Available at: https://www.sodm.nl/binaries/staatstoezicht-op-de-mijnen/documenten/publicaties/2016/06/21/23—cbs—trend-changes-in-ground-subsidence-in-groningen/23-cbs-trend-changes-in-ground-subsidence-in-groningen-update-may-2016-052016.pdf, last accessed 2 October 2017. Statistics Netherlands | 2016 Scientific paper. Pruiksma J., 2015. Isotach formulation of the rate type compaction model for sandstone, Int. J. Rock Mech. Min. Sci. , 78, 127– 132. Roest J., Kuilman W., 1994. Geomechanical analysis of small earthquakes at the Eleveld gas reservoir, Rock Mechanics in Petroleum Engineering, 29–31 August, Delft, Netherlands, doi:10.2118/28097-MS. Ruina A., 1983. Slip instability and state variable friction laws, J. geophys. Res. , 88, 10 359– 10 370. https://doi.org/10.1029/JB088iB12p10359 Google Scholar CrossRef Search ADS   Rutqvist J., Rinaldi A.P., Cappa F., Moridis G.J. et al., 2013. Modeling of fault reactivation and induced seismicity during hydraulic fracturing of shale-gas reservoirs, J. Petrol. Sci. Eng. , 107, 31– 44. https://doi.org/10.1016/j.petrol.2013.04.023 Google Scholar CrossRef Search ADS   Sanz P. et al., 2015. Geomechanical analysis to evaluate production-induced fault reactivation at Groningen gas field, in SPE Annual Technical Conference and Exhibition, 28–30 September, Houston, Texas, USA., doi:10.2118/174942-MS. Segall P., 1989. Earthquakes triggered by fluid extraction, Geology , 17, 942– 946. https://doi.org/10.1130/0091-7613(1989)017%3c0942:ETBFE%3e2.3.CO;2 Google Scholar CrossRef Search ADS   Soltanzadeh H., Hawkes C.D., 2008. Semi-analytical models for stress change and fault reactivation induced by reservoir production and injection, J. Petrol. Sci. Eng. , 60, 71– 85. https://doi.org/10.1016/j.petrol.2007.05.006 Google Scholar CrossRef Search ADS   Sone H., Zoback M., 2014. Time-dependent deformation of shale gas reservoir rocks and its long-term effect on the in situ state of stress, Int. J. Rock Mech. Min. Sci. , 69, 120– 132. Udias A., Madariaga R., Buforn E., 2014. Source Mechanisms of Earthquakes: Theory and Practice , Cambridge Univ. Press. Google Scholar CrossRef Search ADS   Van Eijs R.M.H.E., Mulders F., Nepveu M., Kenter C.J., Scheffers B.C. et al., 2006. Correlation between hydrocarbon reservoir properties and induced seismicity in the Netherlands, Eng. Geol. , 84, 99– 111. https://doi.org/10.1016/j.enggeo.2006.01.002 Google Scholar CrossRef Search ADS   Van den Bogert P., 2015. Impact of various modelling options on the onset of fault slip and the fault slip response using 2-dimensional Finite Element modelling. Report No SR.15.11455, Nederlandse Aardolie Maatschappij, Assen. Available at: https://nam-feitenencijfers.data-app.nl/embed/component/?id=onderzoeksrapporten, last accessed 2 October 2017. Van Thienen-Visser K., Breunese J.N., 2015. Induced seismicity of the Groningen gas field: history and recent developments, Leading Edge , 34, 664– 671. https://doi.org/10.1190/tle34060664.1 Google Scholar CrossRef Search ADS   Van Wees J.D., Mijnlieff H., Lutgert J., Breunese J., Bos C., Rosenkranz P., Neele F. et al., 2008. A Bayesian belief network approach for assessing the impact of exploration prospect interdependency: a application to predict gas discoveries in the Netherlands, AAPG Bull. , 92, 1315– 1336. https://doi.org/10.1306/06040808067 Google Scholar CrossRef Search ADS   Van Wees J.D., Buijze L., Van Thienen-Visser K., Nepveu M., Wassing B.B.T., Orlic B., Fokker P.A., 2014. Geomechanics response and induced seismicity during gas field depletion in the Netherlands, Geothermics , 52, 206– 219. https://doi.org/10.1016/j.geothermics.2014.05.004 Google Scholar CrossRef Search ADS   Wassing B.B.T., Van Wees J.D., Fokker P.A., 2014. Coupled continuum modeling of fracture reactivation and induced seismicity during enhanced geothermal operations, Geothermics , 52, 153– 164. https://doi.org/10.1016/j.geothermics.2014.05.001 Google Scholar CrossRef Search ADS   Zoback M., 2010. Reservoir Geomechanics , Cambridge Univ. Press. © The Author(s) 2017. Published by Oxford University Press on behalf of The Royal Astronomical Society.

### Journal

Geophysical Journal InternationalOxford University Press

Published: Mar 1, 2018

## You’re reading a free preview. Subscribe to read the entire article.

### DeepDyve is your personal research library

It’s your single place to instantly
that matters to you.

over 18 million articles from more than
15,000 peer-reviewed journals.

All for just $49/month ### Explore the DeepDyve Library ### Search Query the DeepDyve database, plus search all of PubMed and Google Scholar seamlessly ### Organize Save any article or search result from DeepDyve, PubMed, and Google Scholar... all in one place. ### Access Get unlimited, online access to over 18 million full-text articles from more than 15,000 scientific journals. ### Your journals are on DeepDyve Read from thousands of the leading scholarly journals from SpringerNature, Elsevier, Wiley-Blackwell, Oxford University Press and more. All the latest content is available, no embargo periods. DeepDyve ### Freelancer DeepDyve ### Pro Price FREE$49/month
\$360/year

Save searches from
PubMed

Create lists to

Export lists, citations