Induced seismicity in the Groningen gas field – arrest of ruptures by fault plane irregularities

Abstract

From dynamic rupture simulations, we reveal under which conditions a rupture in the Groningen gas field stops along fault dip or along fault strike after it starts on a fault in the reservoir. The simulations focus on the capabilities of fault plane irregularities to arrest ruptures. Such irregularities can be recognised in sandstone outcrops. Fault planes in the Groningen field, extracted from the 3D seismic data by seismic attribute extraction methods, show similar irregularities. A detailed surface of a major fault plane in the field indicates that steps and jogs of tenths of metres are possible. Although these irregularities are close to seismic resolution and could be partially artificial, we investigated their effect on rupture arrest.

For typical current stresses in the Groningen field, jogs and steps of this length scale are found to be remarkably effective to stop ruptures in the reservoir. Also, a significant increase in the fault dip along fault strike can stop these ruptures but a kink in the fault under a constant fault dip not.

Including non-planar fault features and pressure diffusion in the Carboniferous, the simulations in this paper follow trends of previous simulations in the literature using 2D planar faults. In particular, the horizontal stress in this formation and the strength of the Carboniferous fault zone are important for rupture propagation. If the fault would have been reactivated in the Neogene or Quaternary and poorly healed in clay-rich parts, rupture propagation into the Carboniferous can only be prevented by jogs of sufficient size and lateral continuity under the present estimate of the horizontal field stress.

Introduction

The Groningen gas field in the Netherlands is the largest onshore gas field in Europe. It suffers from induced seismicity since the early 1990s. The seismicity follows from continuous gas production since 1963 leading to reservoir compaction and more stress on faults in the reservoir. Because of the resulting damage and public unrest in a densely populated area, a large number of scientific studies aimed to understand the different aspects related to induced seismicity in the Groningen field. Making use of a dense shallow borehole geophone and ground accelerometer network in the Groningen area and detailed field and rock data, many reports and articles about the Groningen seismicity have been published.

Despite all this work, it is less clear how the ruptures in the field stop or are arrested. This is important in relation to the possibility of large earthquakes in the field, see NAM (2016) and NAM (2022). Potentially, such earthquakes could be caused or triggered by a substantial rupture in the reservoir. So far, rupture arrest in the field has been addressed for idealised flat or planar two-dimensional (2D) fault planes by Wentinck (2018a), Buijze et al. (2019) and Buijze (2020).

Before starting additional simulations, we applied recent analytical models for rupture arrest on three-dimensional (3D) planar faults from Galis et al. (2019) for typical but simplified stress conditions in the Groningen reservoir and Carboniferous underburden. According to these models, rupture penetration into the Carboniferous cannot be ruled out when the cohesion strength along the fault would be low because of fault reactivation.

Fortunately, real fault planes are not flat but have kinks and transitions in fault dip along fault strike, jogs and steps, and it is known that such geometrical irregularities are capable to arrest ruptures. For this reason, we included these features in a limited number of generic dynamic rupture simulations presented in this paper using realistic field and rock data. The simulations also include pressure diffusion into the lower part of the Zechstein and Upper Carboniferous following from gas production, non-uniform rock properties and Ohnaka’s constitutive model for fault slip. Some relevant field data for these simulations, including a detailed fault surface of one major NW–SE fault in the central part of the field, are presented first.

Field data

The 30 × 40 km gas reservoir at about 3 km depth comprises aeolian and fluvial-deltaic sandstones of the Upper Rotliegend Group (Permian) and intercalated shale layers of the Ten Boer Claystone and Ameland Claystone, see de Jager and Visser (2017). In the southern part of the field, the reservoir also includes fluvial sandstones of the Upper Carboniferous. Thick Zechstein evaporites act as seal for the reservoir. The reservoir thickness gradually increases towards the North from ∼0.1 to 0.3 km. The subsurface geometry of the Groningen field has been determined by NAM with high lateral resolution down to the horizon between the reservoir and the Carboniferous underburden. Prominent are major, primarily NW–SE extensional fault zones in the centre and SW part of the field, a few of them with considerable vertical throw along several subsided and uplifted blocks, elongated in NW–SE direction. Most of these faults align with the regional NW–SE striking trough between two early Carboniferous (Dinantian) highs in the SW and NE parts of the field suggesting that these highs influence the fault structure within the reservoir. Many of the NW–SE faults extend below the reservoir to the Lower Carboniferous strata but not into the ductile Zechstein salt above the reservoir and this also holds for several E–W faults, see Kortekaas et al. (2021). Forming a large negative flower structure, several NW–SE faults gradually bend with a lower dip angle at greater depth.

Several large E–W striking faults with relatively little throw at top reservoir are part of a system of E–W striking fault systems in this part Europe. When crossing the NW–SE trough, the E–W striking faults seem to disappear, at least at reservoir level. Together with several pop-up blocks, the discontinuities in the E–W striking faults may be an expression of compressional step-overs along the NW–SE faults under the N–S tectonic compression.

The faults in the field experienced a long complex tectonic history, see e.g. Doornenbal and Stevenson (2010) and de Keijzer (2008). Relevant for this work and affecting the current fault strength and field stress in the Groningen field is the Cenozoic tectonics, parallel and in front of the Alpine arc including the ∼0.1 km uplift of the Rhenish Massif ∼300 km south of the Groningen field during the Quaternary. The analysis of river terraces indicates an uplift wave, elongating parallel to the Alpine chain and migrating from the southern margin of the massif in the early Pleistocene to the northern Ardennes and Eifel today, see Demoulin and Hallot (2009) and Demoulin and Bourbon (2022). According to these authors, the uplift mechanism seems mainly related to lithospheric folding from the Alpine compression rather than to a local mantle plume under the Eifel. In the Groningen field, reverse fault activation by the Alpine compression was probably minor, leading to some pop-up blocks between NW–SE and E–W striking faults, see de Keijzer (2008). Some inversion during this period from N–S compressive forces cannot be excluded as isochore maps of the formations in the field show local thinning of the upper Paleogene sediments, see Logeman (2017), herein Fig. 29.

In the past 5 years, many other small faults with minor throw in the field have been identified from the seismic data using the Ant-tracking seismic attribute from Petrel^TM, see Kortekaas and Jaarsma (2017). Detailed fault surfaces at reservoir level with extension into the deeper strata below the reservoir have been extracted from the reprocessed and depth-imaged 3D seismic volume of the Groningen field.

Figure 1 shows the central part of the field with faults from the NAM fault database and earthquake epicentres until May 2022 and highlights the variation of the fault dip over the field. More detailed and statistical data about fault strike orientation, fault dip and fault throw can be found in, e.g. de Keijzer (2008), Kortekaas and Jaarsma (2017) and Wentinck and Kortekaas (2022). Figure 2 shows a detailed surface of one of the major NW–SE faults in the field as derived from Ant-tracking. Figure 3 shows cross-sections of this fault at various depths. The figure highlights a few kinks along fault strike and possible expressions of jogs along fault dip noting that these irregularities are close to seismic resolution and could be partly artificial, see Wentinck and Kortekaas (2022). Even so, such features can be seen in outcrops of comparable sandstone formations, e.g. in the national parks in the western part of the US and are known to occur in transtensional faulting, see e.g. Fossen et al. (2013). Further, jogs may originate just from normal faulting when planes of weakness parallel to the bedding plane fail, first forming discontinuous wing cracks which later connect to form a fault plane, as shown on a small scale in tri-axial anisotropic rock failure tests, see Acosta and Violay (2020). So far, we have not correlated the appearance of possible jogs with important changes in the lithology.

Figure 1. Central part of the Groningen field with most induced seismicity. The part shown covers an area of 20 $\times$ 20 km. The faults at the top of the Rotliegend reservoir are shown by the coloured lines and originate from the NAM fault database which includes also data about fault dip angle and fault throw. The colouring elucidate variations in the fault dip angle along fault strike in the range 60–90 $^\circ$ related to the last 3D simulation for a fault with variable fault dip along fault strike, see Fig. 12. The blue and red dots show the epicentres of earthquakes 1.5 ≤ M_L < 3 and M_L ≥ 3 until May 2022, respectively. The smaller blue ones show the M_L < 2.0 earthquakes. The epicentres are from the KNMI public earthquake database unless they were improved by Spetzler and Dost (2017), Willacy et al. (2018), Dost et al. (2020) and Smith et al. (2020), see for details Wentinck and Kortekaas (2022), herein §2.4. The dashed black line indicates the location of the fault plane shown in detail in Figs 2 and 3. The thick red arrow shows the direction of the regional maximum horizontal stress of about 340° with respect to the North. The X- and Y-coordinates are from the Dutch Rijksdriehoeksstelsel coordinate system.

Figure 2. Detailed surface of a major NW–SE fault in the central part of the Groningen field of about 7 km length and over a depth of 2.7–3.3 km. The surface is constructed using the Ant-tracking seismic attribute from Petrel $^{{\rm TM}}$ , explained in Wentinck and Kortekaas (2022), herein App. E. The green surface gives an impression where the fault crosses about the reservoir between a depth of 2.9 and 3.1 km. The fault is indicated in Fig. 1 by the dashed black line. The fault shows possible expressions of sharp kinks along fault strike and of jogs along fault dip, further elucidated by cross-sections at a few depths in Fig. 3.

Figure 3. Cross-sections derived from Ant-tracking at about 2.9, 2.95, 3.0, 3.05 and 3.1 km depth of the NW–SE fault, shown in Fig. 2. The blue and red lines are 20 m below and 20 m above these depths. The solid and dashed lines show the mean and minimum/maximum Y-coordinates of the fault plane cross-sections, respectively. The blue arrows point to significant separations up to $\sim$ 100 m between the solid blue and red lines of tenths of metres that could be expressions of compressional jogs. Although not directly visible, steps may be associated with two opposing kinks relatively close to each other. For reference: the reservoir is between a depth of 2.9 and 3.1 km. The base of the Ten Boer claystone, the Upper and Lower Slochteren are at about 2.94, 3.02 and 3.15 km depth, respectively, see also Table 1 and Figs. 5 and 6.

Table 1. Input parameters for the analytic and dynamic rupture models. The depths, field stress and pressure gradients and friction coefficients are considered to be typical for the reservoir and Carboniferous underburden in the central part of the Groningen gas field.

¹ This value is slightly lower than Sanz et al. (2015) and van den Bogert (2018b) use for their models which is 16 MPa/km leading to a $\sigma_{\text{h}}$ = 48 MPa at 3 km depth.

Most of the M_L ≥ 1.5 earthquake hypocentres are in the reservoir and coincide with faults from the NAM fault database, and a substantial number of them are on the major NW–SE faults. Most rake angles indicate fault slip predominantly in vertical direction. For the main period of induced seismicity 1995–2021, we observe no trends over time in fault dip or fault throw of fault segments nearest to the earthquakes epicentres although the induced stress from gas production increases with time and systematically varies with fault dip and fault throw, see Wentinck and Kortekaas (2022), herein App. F. Estimates for the stress drop during rupture Δτ and rupture plane dimension L have been derived by Ameri et al. (2020) and Dost et al. (2016). The stress drop Δτ = τ ₀ − τ _r [Pa] is the difference between the stress on the fault τ ₀ before rupture and the high-velocity residual slip resistance τ _r during rupture. They supposed that the earthquakes originate from a circular rupture plane with radius R [m] and that the rupture propagates in the fault plane with the shear velocity. In addition, we reprocessed data from ground motions of a number of earthquakes in a similar way and the derived stress drops agree quite well with those of Ameri et al. (2020), herein Fig. 5b, see Wentinck and Kortekaas (2022), herein Fig. 14. For M _L > 2 earthquakes, Δτ = 1–5 MPa. Because of strong dependency of Δτ on the supposed rupture velocity, Δτ values can easily vary with a factor 2–3. To facilitate a comparison of the various results, the relation between the mean corner frequency and rupture plane dimension in these references was based on Brune’s source model but note that for non-symmetric ruptures, the constants used herein can differ, see e.g. Kaneko and Shearer (2014).

From variations in the ground motion spectra over the source–receiver azimuth, Ameri et al. (2020) and Oates et al. (2020) conclude that several ruptures in the field have propagated predominantly in one direction along fault strike. For the Zeerijp M_L 3.4 earthquake in 2018, the same was concluded from variations in the lower corner frequency and apparent source time function durations over the source–receiver azimuth by Wentinck (2018b). An indication for unidirectional rupture propagation along fault strike implies the existence of distinct rupture barriers along fault strike that stop the rupture in the opposing fault strike direction. From the relative variation of the duration of these functions over the source–receiver azimuth, the rupture velocity along fault strike was roughly estimated ∼1 km/s.

Rupture arrest

Before starting with dynamic rupture simulations, we applied recent analytical rupture arrest and runaway models of Galis et al. (2019) to determine which parameters are most relevant for rupture arrest and runaway in the Groningen field, see Wentinck and Kortekaas (2022). These models are for uniform and planar fault planes with 2D rectangular asperities or nucleation zones of various shape. The models result from a balance between the mechanical energy released by the rupture and the fracture energy, required to generate fractured surface near the fracture tip at the circumference of the rupture plane. While the rupture develops, the mechanical energy available and transported to the fracture tip changes with the rupture velocity and with the size of the fracture formed. The energy balance includes the important breakdown length scale L _f [m] over which the stress breaks down around the fracture tip, i.e. L _f = 2μG _c /τ _b ² where G _c [J/m²] is the fracture energy per unit area required, τ _b [Pa] the breakdown stress drop and μ [Pa] is the shear modulus of the rock surrounding the fault zone. τ _b = τ _p − τ _r is the difference between the peak or maximal slip resistance of the fault plane τ _p [Pa], and the high-velocity residual slip resistance τ _r [Pa]. G _c can be expressed as G _c = 1/2D _c τ _b . This expression is exact for linear strain-weakening slip behavior and holds approximately for Ohnaka’s constitutive model for slip, see Ohnaka (2013), herein §2.2.5.

The analytical models of Galis et al. (2019) align with those found in textbooks, see e.g. Aki and Richards (2009) and Udias et al. (2014) which are based on work of Kostrov (1966) and Husseini et al. (1975) among others. In the case that the breakdown stress τ _b equals the stress drop Δτ, the fracture energy per unit length required G _c [J/m²] must increase, if stepwise, with δG _c > G _c (L − L _f )/L _f to arrest a 1D rupture of length L, see e.g. Husseini et al. (1975) or Aki and Richards (2009), herein Ch. 11. If L ≫ L _f , δG _c /G _c > L/L _f and only a substantial fracture energy barrier can stop such ruptures on a flat fault plane. Alternatively, a rupture can stop by a reduction of the load on the fault which leads to a reduction of the stress drop. But it can be shown that this reduction must be also substantial, see Aki and Richards (2009), herein Ch. 11. For typical values for the Groningen reservoir for 2 < M_L< 3 earthquakes, D _c = 0.01 m, μ = 3 GPa and τ _b = 5 MPa, we have L _f ∼ 6 m. According to these expressions, a 2D rupture over 100 m or more in the reservoir would require a substantial increase in the fracture energy when penetrating into the Carboniferous before it stops.

The other important model parameter for rupture arrest in the analytical models is the so-called strength parameter S = (τ _p − τ ₀)/(τ ₀ − τ _r ) [ − ] where τ ₀ is the tangential stress or load on the fault plane. The related parameter S + 1 = (τ _p − τ _r )/(τ ₀ − τ _r ) compares the load or stress drop Δτ = τ ₀ − τ _r [Pa] relative to the breakdown stress drop τ _b = τ _p − τ _r . For a critically loaded fault, τ ₀ → τ _p and S→ 0.

The input data for the analytical models as applied to the Groningen field conditions are extensively explained by Wentinck and Kortekaas (2022), herein Ch. 3. The data have been derived from several recent reports and articles and recently processed seismic data, see below. We concluded that M_L > 2 earthquakes in the past decades on faults with dip angles <70^∘ and with a relatively low cohesion strength in the fault zone had potential to penetrate into the Carboniferous. A low cohesion strength in the Carboniferous NW–SE faults could be a result of reactivation of clay-rich parts of these faults during the N–S tectonic compression in the Quaternary if these faults would have functioned as compressional-stepovers. Hence, we performed the 3D simulations below primarily for this unfavourable condition.

Dynamic rupture simulations

2D dynamic rupture models for induced seismicity in the Groningen field have been presented by Buijze et al. (2020), Wentinck (2018a), van den Bogert (2018a), van den Bogert (2018b), Buijze et al. (2019) and Buijze (2020) for a broad range of input parameters. These models for ruptures on planar faults show that ruptures mostly stop in the reservoir but can also propagate into the Carboniferous underburden. Rupture propagation into the Carboniferous is promoted by a large tangential stress on the fault plane, a large stress drop, a low strength or small fracture energy and little or no fault throw.

Considering the difficulty to rule out poorly healed faults in the Carboniferous with a low cohesion strength under unfavourable stress conditions in the field, the purpose of the present dynamic rupture simulations is to investigate the effectiveness of geometrical irregularities to arrest ruptures. In particular, we present here 2D and 3D generic simulations for faults with jogs along fault dip and steps, kinks and fault dip transitions along fault strike. The simulations include non-uniform reservoir rock and a macroscopic constitutive model for slip resistance from Ohnaka (2013), herein §4.3. This model includes the cohesion strength in the fault zone and a smooth transition from elastic to non-elastic deformation. In our case for predominantly normal faulting and following Scholz (2002), herein §3.5 and Udias et al. (2014), herein §11.5, jogs are deviations along fault dip and steps are deviations along fault strike. Jogs induce additional compressive or extensional forces on the fault during slip, and these irregularities are recognised to be effective barriers to impede or stop ruptures.

Set-up

The 2D and 3D simulations have been done with the finite element method solver COMSOL®. Details about the set-up of this solver are given in Wentinck and Kortekaas (2022), herein App. C. For convenience, we summarise here key solver features. The mechanical behaviour of the rock outside the fault zone is solved by the ‘Solid Mechanics’ module under default settings. ‘Boundary Similarity Component Couplings’ between the hanging and foot walls of the fault are used to calculate the relative displacements or slip between these planes. The couplings are used to calculate the normal and tangential forces on both planes. The poro-elastic pressure response is included using the ‘Thermal Expansion’ submodule, replacing the expression for thermal expansion by an equivalent poro-elastic one. The constitutive slip resistance model and rock properties are implemented as explicit analytical functions using ’Variables’ and ’Analytical Functions’.

From seismics, we expect jogs and steps with sizes of several tenths of metres in the Groningen field as indicated by Fig. 3. These fault plane irregularities and also kinks along fault strike have been modelled using 1D and 2D parametric sigmoid-type functions in ‘Parametric Surfaces’. For example, a cross-section of a jog along fault dip is defined by adding the second term in x = x(s) = x ₀ + (s − z ₀)cotan(δ _main) − w _jog/(1+ exp((s − z _jog)/δz _jog*)). Herein, x and z are the fault plane coordinates in the xz-plane and s is the parameter. δ _main is the dip angle of the main part of the fault away from the jog, w _jog [m] is the width of the jog and z _jog [m] the depth of the jog. The parameters δz _jog* [m] and w _jog determine the minimum dip or slope angle of the jog δ _jog* [degr.]. Such functions lead to a sufficient smooth fault plane to avoid numerical problems. Herewith, the mesh size along the fault plane follows not only from the usual rupture propagation requirements in relation to the rupture length scale L _f but also from capturing the curvature of the fault plane irregularity properly. For a few cases, the mesh was refined and coarsened to be confident on a minor impact on arrest conditions.

The simulations consist of two time-dependent ‘Studies’. In the first study using a variable time step, the fault is quasi-statically loaded by the field stress and the pressure drop in the reservoir and Carboniferous until fault slip starts. Nucleation of a rupture is detected automatically by the solver, monitoring a local steep increase of fault slip. In this work with a focus on rupture arrest and runaway, the nucleation of the rupture around a predefined stress condition is realised by including a relatively weak patch in the fault plane in the Upper Slochteren formation in the reservoir with an artificial lower strength of cohesion. After nucleation, the second study starts with a maximum time step of 1 ms over a period of 0.2–0.5 s.

All faults simulated have a typical moderate fault throw of 40 m and a fault dip of 70^∘ unless mentioned otherwise. The forces on the upper and lower fault planes are determined by the relative displacement or slip in the fault zone between these two planes. The subsurface formations attached to the upper and lower fault planes are shifted upwards and downwards half the value of the fault throw, respectively. The input parameters to model the field stress on these faults and rock properties have been derived from several recent reports and articles and from recently processed seismic data and have been extensively discussed in Wentinck and Kortekaas (2022), herein Ch. 2. The field stress gradients used originate from van Eijs (2015). The vertical field stress has been derived from the gravitational load of the sediments. The horizontal field stress in the field results from the gravitational load, the tectonic Alpine compression and ridge push from the Mid-Atlantic Ridge. It increases with depth except at the top of the reservoir where it drops with several MPa as the Zechstein salt above the reservoir has a much higher Poisson ratio. The maximum horizontal stress has been inferred from differential strain relaxation measurements in cores, sonic logs in boreholes, borehole break-outs and a seismic noise interferometry measurement. The horizontal field stress anisotropy is modest and the combined field and induced stress favours normal faulting and broadly agrees with the regional mean horizontal stress.

Conservatively, we suppose that the minimal horizontal field stress gradient also holds for the Carboniferous shale. However, a more favourable higher value is possible in this clay-rich formation because of viscous differential stress relaxation over geological time. According to minifrac tests in the BLF-104/105 wells near Blija north-west of the Groningen field, the minimum horizontal field stress gradient at ∼ 2.7 km depth is ∼ 10% higher in the Ten Boer claystone than in the underlying Slochteren sandstone, see NAM (2022), presentation L. Buijze. Another example is the Wolfcamp shale sequence in the Permian Midland basin in West Texas. At ∼ 2.4 km depth, the increase in the minimum horizontal field stress gradient is ∼ 10–15% in parts with ∼ 35 wt% clay when compared to more sandy strata, see Kohli and Zoback (2021).

The reservoir pressure before production was ∼ 35 MPa, ∼2 MPa above a hydrostatic pressure of 33 MPa at 3 km depth. Gas production started 1959 and the reservoir pressure declined after 1970 to about 8 MPa in the central part of the field. The pressure difference over most faults in this part of the field was less than 1 MPa. Because of excellent reservoir permeability, the pressure is practically uniform over depth in the reservoir, increasing hydrostatically below the gas–water contact. Over the production time until 2020, the pressure drop probably has diffused tenths to hundreds metres into the Carboniferous underburden.

Rock data

The Slochteren sandstone from the Stedum and Zeerijp (ZRP) wells has been rigorously analysed under in situ conditions in tri-axial cells by Hol et al. (2018), Pijnenburg (2019) and Hunfeld et al. (2020). For a considerable number of reservoir and Carboniferous cores, the quasi-static uniaxial compression modulus C _m [Pa ^{− 1}], Poisson ratio ν [ − ] and onset of inelastic deformation have been determined. These sandstone properties quite depend on the porosity ϕ [ − ], see Fig. 4. Furthermore, C _m has been derived from strain in the reservoir using a distributed optical fibre cable over 18 and 41 months, see Kole et al. (2020). The data agree well with those derived from the core data.

Figure 4. Left: Uniaxial compaction coefficient C_m [1/Pa] of reservoir sandstone versus porosity $\phi$ . The experimental data shown originate from core data from the ZRP-3a well from Shell Global Solutions International (SGS-I) and from Pijnenburg (2019) (UU) (black dots) and from strain measurements by the distributed strain sensing (DSS) optical fibre cable in the ZRP-3a well (green dots), see Kole et al. (2020). The red line follows from the correlation for the apparent uniaxial compaction coefficient from van Eijs and van der Wal (2017), i.e. $C'_{m}(\phi) = 10^{-4}(267.3\phi^3 - 68.72\phi^2 + 9.85\phi + 0.21)$ with $C'_{m}$ values in MPa $^{-1}$ . In the simulations, we use $C_m = C'_m/\alpha$ where the Biot constant $\alpha$ [−] is calculated from the elastic constants as explained in Wentinck and Kortekaas (2022), herein App. A. Right: Unconfined compressive strength UCS [Pa] of reservoir sandstone versus porosity as derived from a scratch test on cores from the ZRP-3a well (dots). The strength steeply drops for $\phi \gt $ 17%. The blue line is an empirical fit function used in simulations, i.e. $\text{UCS}(\phi) = (4 + 20 \exp(-\text{max}(0,(\phi-6))/8)^{1/0.15}/20))$ . It originates from a combination of two reasonable correlations or empirical fits between the sandstone permeability $k$ [mD] and porosity $\phi$ and the sandstone permeability and the unconfined compressive strength, i.e., $\phi(k) = 6 + 8 k^{0.15} \quad \text{and} \quad \text{UCS}(k) = 4 + 20 \exp(-k/20)$ . Since the fit between permeability and porosity is not well constrained, another reasonable fit UCS( $\phi$ ) was explored and resulted in similar arrest conditions, see Wentinck and Kortekaas (2022), herein App. A and Fig. 10.

Reservoir compaction involves both poro-elastic and time independent, permanent strain. The latter follows from consolidation and shear of clay films between the sandstone grains, while grain failure occurs at higher stresses. This elasto-plastic behaviour is well described by a Modified Cam-Clay (MCC) model. But in the models in this paper, we use for the reservoir sandstone and the Carboniferous shale elastic properties, which are directly derived from the compaction data. According to Buijze (2020), herein §3.4, these so-called ‘apparent elastic’ properties mimic the MCC behaviour quite well in dynamic rupture simulations.

From scratch tests on sandstone cores of the ZRP-3a well, the unconfined compressive strength UCS [Pa] is in the range 5–25 MPa. Figure 4 shows the unconfined compressive strength UCS versus porosity. The strength steeply drops for ϕ> 17%. The corresponding strength of cohesion S ₀ of intact Slochteren sandstone is in the range 1–10 MPa, using S ₀ ∼ UCS/3. The sandstone properties considerably vary with the porosity, and the porosity considerably varies with depth. In the central part of the field, the average porosity in the reservoir is 18–22%, and there is lateral continuity of porosity streaks over hundreds of metres.

From scratch tests on Carboniferous cores from the ZRP-3 well, the unconfined compressive stress UCS = 50–100 MPa, see van der Linden et al. (2020). The lower values are for relatively weak clay-rich zones with more than 63 wt% clay. The higher ones are for quartz-rich zones with less than 58 wt% clay. The corresponding cohesion strength S ₀ varies in the range 20–30 MPa. The Carboniferous cores have a porosity of 1–3%, and the uniaxial compaction coefficient C _m is in the range 2–8 10 ^{− 5} MPa ^{− 1}. The data are consistent with values derived from in situ strain measurements in the Carboniferous, see Kole et al. (2020), herein Fig. 4. The mechanical properties of the rather heterogeneous Ten Boer claystone in the upper part of the reservoir are in between those of the Carboniferous shales and Slochteren sandstones, but more similar to the Carboniferous shales. For intact Ten Boer claystone, UCS = 50–100 MPa.

The peak or maximal slip resistance of the fault plane follows from τ _p = S ₀ + μ _s ${\sigma_{n}^{\prime}}$ [Pa] where ${\sigma_{n}^{\prime}}$ = σ _n − p [Pa] is the effective normal stress with respect to rock failure. μ _s [ − ] is the quasi-static or static friction coefficient and S ₀ [Pa] is the strength of cohesion in the fault zone. In general, the strength of cohesion in the fault zone can differ considerably from the one in the intact rock around it. At least it depends much on the healing of the fault after fault reactivation. This healing is a slow process and depends on many factors, such as pressure, temperature and rock and pore fluid chemistry. For the sandstone reservoir, we assume that the fault zone is healed in line with healing experiments of Hunfeld et al. (2020) and that, without further knowledge, the cohesion strength is the same as in the corresponding intact rock around it. Expecting that a rupture propagates through the weakest material in the fault zone and the actual slip zone is thin compared to the fault zone, we assume that primarily the lower value of an unknown fault zone strength and nearby intact rock strength matters. Hence, in the reservoir, S ₀ follows from the aforementioned core scratch tests and from the resistive stress just before failure in tri-axial cells, see Fig. 4. For faults with fault throw, we choose to calculate S ₀ = S ₀(ϕ) from the mean porosity of the rock on both sides of the fault plane.

Extrapolating the healing experiments of Hunfeld et al. (2020) to the clay-rich Carboniferous, it cannot be ruled out that a fault in the Carboniferous, activated a few millions year ago, was not healed. It is for this reason that we also investigated rupture propagation and arrest in this formation for a few cases that the cohesion strength and herewith the peak resistance are much lower than could be expected from intact rock.

The macroscopic constitutive model, which relates the mechanical or tangential slip resistance of the fault zone τ [Pa] to the relative displacement or slip D [m] over the fault zone, is from Ohnaka (2013), herein §4.3, Eq. 4.20. This model is a bit more complicated than the frequently used linear strain-weakening model because it allows for a natural transition between linear and non-linear resistance and includes fault strengthening before fault failure. The input parameters are the breakdown slip distance D _c [m] over which the slip resistance decreases to the residual value and the peak and high-velocity residual resistances τ _p and τ _r [Pa]. The latter are calculated using the quasi-static and high-velocity residual friction coefficients μ _s and μ _r and the aforementioned strength of cohesion S ₀.

As said, τ _p = S ₀ + μ _s ${\sigma_{n}^{\prime}}$ . The high-velocity residual slip resistance follows from τ _r = μ _r ${\sigma_{n}^{\prime}}$ . According to measurements on representative fault gouge powders, μ _r is in the range 0.2–0.3 for the Rotliegend reservoir sandstone and Carboniferous shale, see Hunfeld et al. (2021). The other input parameters are the ratio τ _i /τ _p between the slip resistance at the onset of inelastic deformation τ _i and the peak resistance τ _p and the ratio D _i /D _c between the slip distance at the onset of inelastic deformation D _i and the breakdown distance D _c . τ _i /τ _p is in the range 0.7–0.9 corresponding to similar ratios for quasi-static deformation of intact sandstone and Carboniferous cores under tri-axial loading, see e.g. Pijnenburg (2019), herein Fig. 2.3 and 3.3 and van der Linden et al. (2020), herein Fig. 4. D _i /D _c is in the range 0.1–0.4. D _{p, inel}/D _c ∼ 0.2 where D _{p, inel} = D _p − D _i hardly changes with the effective normal stress and peak resistance, in line with other rock failure data, see Ohnaka (2013), herein §4.4.

For pre-existing faults, D _c depends on the breaking of interlocking asperities in the irregular fault plane, ploughing of hard rock fragments in softer rock, fault zone dilatation and on cohesion forces between the grains in the fault zone. As a consequence, D _c scales with these self-similar irregularities in the fault zone. According to theory and observations of natural earthquakes, these irregularities, and herewith D _c , scale with the seismic moment M₀ as D _c ∝ M₀ ^1/3 over 6–7 orders of magnitude, see Ohnaka (2013), herein Fig. 5.21. For 1.5 ≤ M_L ≤ 3.5 earthquakes of interest, M₀ varies in the range 1–1000 TJ and typical values for D _c increase from ∼0.3 to ∼3 cm. Following common practice in dynamic rupture modelling, we suppose D _c is constant although it depends on the magnitude of slip or the magnitude of the earthquake. Since this work is about rupture arrest of representative 2 < M_L < 3 earthquakes and not about rupture nucleation, this is to our opinion not a serious simplification. We use D _c = 1 cm. Figure 7 shows τ _res for typical effective normal stresses and two extreme values of the rock porosity.

The formations on both sides of the fault have porosity profiles over depth and in lateral direction representative for the central part of the Groningen field. A detailed well log of the ZRP-3a well in the Slochteren sandstone has been used to generate porosity variations along depth, see Fig. 5. The rock properties vary with the porosity as shown in Fig. 4 and lead to a heterogeneous reservoir and Carboniferous with realistic lateral porosity variations. The cohesion strength along the fault follows from S ₀ = S ₀(ϕ). In the Lower Slochteren formation in the reservoir where ϕ< 17%, S ₀∼10 MPa. Figure 6 shows the geometry of a fault with a single step along fault strike and the porosity variations used in the 3D simulations.

Figure 5. Left: porosity variations in the 210 m thick Slochteren sandstone in the reservoir over a true vertical depth (TVD) of 2.94–3.15 km (blue line). The data are from a well log in the ZRP-3a well with sample intervals of 0.15 m. This porosity profile has been used to calculate representative ones along fault dip, such as shown on the right, as explained in Wentinck and Kortekaas (2022). The two red dots indicate the depth of the top and the base of the reservoir. Right: porosity $\phi$ (black line), corresponding shear modulus $G$ (blue line) and cohesion strength S ₀ (red line) versus depth along the foot wall of the fault used in one of the 2D simulations. Because the fault has a throw of 40 m, the porosity profile is shifted 20 m upwards relative to the porosity profile in the left figure. The maximal value of S ₀ is taken 12 MPa.

Figure 6. Foot wall of normal fault with one step along fault strike of 10 m width in the reservoir and upper part of the carboniferous down to a depth of 3.3 km. The formations indicated are the ten Boer claystone (ROCLT) at 2.88–2.94 km depth, the Upper Slochteren sandstone (ROSLU) at 2.94–3.02 km depth and the lower Slochteren sandstone (ROSLL) at 3.02–3.15 km depth. The fault dip is 70 $^\circ$ , fault strike azimuth is 120 $^\circ$ and fault throw is 40 m. The blue colour contouring indicates the porosity.

A 1D analytical solution has been implemented in the solver to calculate the pressure drop diffusion into the Zechstein and Carboniferous. Assuming that the Carboniferous shale is fully saturated with brine and a porosity and permeability in vertical direction ϕ∼ 0.05 and k ∼ 100 nD, the pressure diffusion coefficient used is D _p = k/(ϕμ _f κ _f ) ∼ 8 10 ^{− 6} m²/s where κ _f [Pa ^{− 1}] and μ _f [Pa.s] are the compressibility and viscosity of the pore fluid, respectively. For 7 Molar or 38 wt% salinity brine at an in situ temperature of 100 ^∘C, μ _f ∼ 0.5 mPa.s and κ _f ∼ 5 10 ^{− 10} Pa ^{− 1}. A pressure drop at the bottom of the reservoir penetrates into the Carboniferous in a period Δt over a distance $\sqrt{D_p \Delta t}$ . Over 50 years, ∼0.1 km.

A permeability of 100 nD for this type of shale is not unusual, but it could also be a factor 10 higher or lower. The value for D _p does not conflict with the pressure drop estimated in the Carboniferous near Zeerijp by Kole et al. (2020), see Wentinck and Kortekaas (2022), herein App. D. We note that de Zeeuw and Geurtsen (2018), herein Ch. 4, use a considerably higher vertical permeability in the Carboniferous regarding strong heterogeneities in this formation. According to Kole et al. (2020), herein Fig. 7, also the Zechstein shows a noticeable compression albeit minimal in the dense anhydrites of the Basal Zechstein. Since the ratio between the vertical and horizontal field stresses is less than in the reservoir and fault healing in the Zechstein is fast, rupture runaway into the Zechstein is not likely. Pressure diffusion in this almost impermeable formation is modelled similarly as for the Carboniferous but with more uncertainty about the permeability and porosity. Assuming that pores in the Zechstein are filled with a similar high-salinity brine as in the Carboniferous and taking ϕ = 3% and k = 1 nD, we use D _p ∼ 1.3 10 ^{− 7} m²/s. Because of the presence of permeable sandy layers in the Ten Boer claystone, small possible pressure differences between the Slochteren sandstone and the Ten Boer claystone are ignored (Fig. 7).

Figure 7. Left: slip resistance $\tau _{{\rm res}}$ . Right: effective friction coefficient $\mu _{{\rm res}} = \tau _{{\rm res}}/\sigma '_n$ . The peaks of the curves equal the peak resistance $\tau _p = S_0 + \mu _s\sigma '_n$ . The difference between the values of the two red peaks for the 5 and 25% porosity curves correspond to the difference in S ₀. The values at large slip are equal to the high-velocity residual slip resistance $\tau _r = \mu _r\sigma '_n$ . Both figures have been calculated for the fault zone at 3 km depth in the reservoir rock using Ohnaka’s constitutive model, see Ohnaka (2013), herein §4.3. The parameters are given in Table 1. The effective normal stress $\sigma'_n$ is 25, 30 and 35 MPa. For a fault with a dip angle of 70°, this corresponds to a reservoir pressure of 18.8, 13.8 and 8.8 MPa. The fault zone consists of healed rock of 5 and 25% porosity and a corresponding cohesion strength S ₀. The latter is shown as a function of porosity by the solid line in Fig. 4, using $S_0 \sim \text{UCS}/3$ . The transition from elastic to inelastic slip resistance is supposed to occur at 80% of the peak resistance $\tau_p$ and 10% of the breakdown slip distance D_c. Using $\eta = 3/D_c$ , the parameters $\alpha$ and $\beta$ in Ohnaka’s constitutive model for slip are calculated from these values.

Results of 2D simulations

Figure 8 shows the tangential slip resistance and effective normal stress on this fault for the aforementioned reservoir pressures just before rupture and the corresponding strength parameter S. Before rupture the tangential slip resistance is equal but opposite to the tangential or shear stress on the fault plane. Despite a considerable reservoir pressure drop of ∼10 MPa after 1991, the reduction in S is modest relative to the one before 1991. The tangential stress on the faults in the reservoir considerably varies with the porosity. It considerably differs from the stress calculated for uniform reservoirs with no pressure depletion in the Carboniferous, as has been done in earlier simulations. In particular, the present tangential stress profile lacks the two stress peaks at the base of the reservoir.

Figure 8. Results of 2D simulations on a planar fault with no jog with a $70^\circ$ dip angle and 40 m throw and for three different reservoir pressures. These are 35, 20 and 10 MPa and correspond in time with 1958 (black lines) before production, 1992 just before the onset of earthquakes in the field (blue lines) and recent years (red lines), respectively. The resulting stress conditions are typical for all 2D and 3D simulations. The two red dots in plot a indicate the depth of the top and the base of the reservoir. From left to right: (A) tangential slip resistance $\tau_{\text{res}}$ of foot wall (pointing upwards when positive), (B) effective normal stress in relation to failure $\sigma'_n$ , (C) the stress drop $\Delta \tau$ (solid lines) and breakdown stress drop $\tau_b$ (dashed lines) and (D) the strength parameter $S$ along the foot wall of the fault. The corresponding porosity profile is shown in Fig. 7 and the cohesion strength in the fault zone $S_0 = S_0(\phi)$ follows from the solid line in Fig. 4, using $S_0 \sim \text{UCS}/3$ . The tangential slip resistance of the foot wall is upwards when positive. The initial stress variations in the formations above the reservoir originate from variations in the Poisson ratio, see also Table 1. In the lower and upper Zechstein above the reservoir, the tangential stress on the fault of 1–2 MPa due to the high Poisson ratio in these formations. This leads to a non-physical negative stress drop $\Delta \tau$ in plot C and is not shown. Despite a considerable reservoir pressure drop of 10 MPa after 1991, the reduction in $S$ (and increase in fault loading) is modest when compared to the reduction in $S$ before 1991. For this fault with a considerable cohesion strength $S_0(\phi)$ for $\phi \lt $ 17% in the Lower Slochteren, the breakdown stress drop $\tau_b \sim$ 7–12 MPa. This potential stress drop is higher than the observed stress drop $\Delta \tau$ for $\text{M}_{\text{L}} \gt $ 2 earthquakes, see Wentinck and Kortekaas (2022), herein Fig. 19. For a constant S ₀ = 2 MPa (not shown), $\tau_b$ would be 5–7 MPa. In relation to the observed stress drop (equal to the difference between the tangential slip resistance just before slip $\tau_{\text{res}}$ and the high-velocity residual slip resistance), the calculated stress drop $\Delta \tau = \tau_{\text{res}} - \tau_r$ (solid lines) is more relevant. The calculated one is the range 3–5 MPa in the period between 1991 and recent years.

The simulated ruptures preferably nucleate where the induced stress is high, as in the Upper Slochteren sandstone or in sandy layers in the Ten Boer claystone. In the Zechstein, stress conditions are unfavourable for nucleation. According to Fig. 8, the tangential stress on the fault in the Zechstein is relatively low, and the strength parameter S is high as a result of a relatively high Poisson ratio and presumable high cohesion strength of this rock. According to laboratory measurements, reactivated anhydrite and salty fault gouges easily heal and recover in strength.

The maximum slip over the fault zone following from a nucleation at 3.0 km depth is ∼0.2 m when the rupture reaches the bottom of the reservoir at 3.15 km depth. The mean and maximum rupture propagation velocities between 3.0 and 3.1 km depth are 1.1 and 1.6 km/s, respectively, considerably lower than the shear velocity V _s = 2–2.5 km/s in the reservoir. Comparable low values are obtained by Buijze et al. (2019), herein Fig. 9. Also on a much larger scale, substantial variations in V _r over the fault plane have been found for natural earthquakes, see Wen et al. (2012) and Zhang et al. (2019).

According to Fig. 8, the stress drop Δτ is in the range 3–5 MPa between 1991 and recent years. Δτ = τ _res − τ _r is equal to the difference between the tangential slip resistance before slip and the high-velocity residual slip resistance. As concluded before, this can be achieved if μ _r is in the range 0.2–0.3, see Wentinck (2018a). For ruptures penetrating into the Carboniferous, the contribution of the cumulative slip in the Carboniferous to the earthquake magnitude quite depends on the stress drop in this formation.

Rupture arrest on compressional jogs in the lower part of the reservoir and of various dimensions shows the following. The ruptures are either arrested by a jog, hardly affected by it or only impeded by it for a short period. In the latter case, the rupture still propagates towards the Zechstein on the upper side of the nucleation patch during this period. Rupture arrest not only depends on the geometrical properties of the jog but also on the slip, the stress drop and loading of the fault. Rupture arrest by a compressional jog is favoured by a higher fracture energy G _c and a lower load on the fault or, equivalently, a higher strength parameter S around the jog. The main reason for this is the higher normal stress on the jog from an increasing contribution of the vertical field stress.

Rupture arrest on jogs of various dimensions and located 50–200 m below the rupture nucleation patch at 3.0 km depth has been classified in Fig. 9. To generate the largest rupture plane in the reservoir possible, the deeper jogs are up to the base of the reservoir. The width and slope angle with respect to the horizontal plane of the jogs were varied with w _jog = 5–40 m and δ _jog* = 2–50^∘. The data include also simulations done for lower cohesion strength along the fault plane in the range 2–5 MPa.

Figure 9. Tenths of rupture runaway and arrest conditions for downwards propagating ruptures for a fault with a jog along fault dip and with representative porosity profiles along dip. The conditions are plotted against the two dimensionless parameters $\Lambda$ and $\Gamma$ related to geometry and energy, see text for a detailed description of these parameters. The blue dots show arrest, the red ones runaway without time delay and the pink ones runaway with a time delay $\gt$ 100 ms during which the rupture still grows towards the Zechstein above the nucleation zone. The jogs are of various shapes with $w_{{\rm jog}}$ = 5–40 m and $\delta _{{\rm jog}}^*$ = 2–50 $^\circ$ and located 50–200 m below rupture nucleation at 3.0 km depth. The cohesion strength along the fault follows from $S_0 = S_0(\phi )$ or is maximised by $S_{0,{\rm max}}$ = 2–5 MPa. Loading and fault strength variations lead to strong variations in the strength parameter $S$ . The ratio $G_{{\rm jog}} = \Delta G_{c,{\rm jog}}/\bar G_{c,{\rm jog}}$ is in the range 0.3–0.9. The black large dot is a reference marker leading to rupture arrest and for the following typical parameter values: $w_{{\rm jog}}$ = 20 m, $\delta ^*_{{\rm jog}}$ = 10 $^\circ$ , $L^{*}$ = 200 m, $S_{{\rm jog}}$ = 1.0 and $R_{{\rm jog}}$ = 0.8. Rupture runaway occurs for low values of $\Lambda$ and/or $\Gamma$ . To help the eye, the grey hyperbole more or less follows a diffuse boundary between rupture arrest and runaway conditions.

The conditions for rupture arrest and runaway are plotted in this figure against two dimensionless parameters Λ and Γ [ − ], aiming for a simple boundary between these conditions with minimal overlap. Rupture runaway occurs for low values of Λ and/or Γ. Following dominant terms or factors in analytical expressions for rupture arrest from Galis et al. (2019), Λ is based on geometrical parameters and Γ on the fracture energy and the loading on the fault and jog. The relevant geometrical parameters for the jog and the rupture reaching the jog are the width w _jog, minimal dip of the jog δ _jog* and the size of the rupture plane L* [m] when reaching or passing the jog. L* is a good measure of the earthquake magnitude developed at these stages. We choose Λ = w _jog/(T _jog ² L*) where T _jog = tan (δ _jog*) is the tangent of the minimum dip angle of the jog.

Γ [ − ] is based on the relative increase of the fracture energy G _jog and the loading on the fault at and around the jog. G _jog = ΔG _{c, jog}/Ḡ _{c, jog} where ΔG _{c, jog} = max(G _{c, jog}) − Ḡ _{c, jog} and Ḡ _{c, jog} and ̅S _jog are the mean values of G _c and S in 10 m depth intervals 20 m above and 20 m below the jog. G _jog increases when the jog becomes more horizontal, i.e. when the jog dip angle δ _jog* decreases due to the increasing contribution of the vertical field stress to the normal stress on the jog.

Following expressions in the analytical models for rupture runaway, we have included the effect of fault loading by adding the factor (̅S _jog+ 1)² to Γ. ̅S _jog + 1 relates to the breakdown stress drop and the stress drop in the fault zone around the jog Δτ as ̅S _jog + 1 = τ _b /(τ _res − τ _r ) = τ _b /Δτ. For a critically loaded jog, τ _res → τ _p , ̅S _jog + 1 → 1 and ̅S _jog → 0. Herewith, Γ = (̅S _jog+ 1)² G _jog

Considering all variations, Fig. 9 shows a reasonably simple but diffuse boundary between all rupture arrest and runaway conditions examined. Using equal powers for the factors jog and L*, the breakdown length scale L _f cancels in the expression for the parameter Λ. Using different powers for these factors made the boundary not sharper. The load on the fault, expressed by the strength parameter S, has a strong effect on the capability of a jog to arrest a rupture. A jog which could arrest a rupture at moderate loading may fail to do so when the loading further increases (Fig. 10).

Figure 10. Top: snapshots of the slip magnitude $D$ during ruptures in a fault plane with a single vertical step. $D$ is expressed by colour and by a deformation out of the fault plane and proportional to $D$ . The reservoir pressure is 10 MPa. The nucleation of the earthquake at 3.0 km depth and 50 m away from the step is promoted by slightly increasing the porosity in a small patch in the reservoir. In the three cases shown, the breakdown slip distance and the stress condition or strength parameter $S$ have been varied; the latter by varying the high-velocity residual friction coefficient $\mu_r$ . For $\mu_r$ = 0.3 and D_c = 0.01 m, the step stops the rupture, also when the rupture growths further or the nucleation occurs farther away from the step. Reducing the breakdown slip distance from D_c = 0.01 m to 0.008 m, the step impedes the rupture but after a short time the rupture passes the step. The same happens but after a shorter time for $\mu_r$ = 0.25 and D_c = 0.01 m. Bottom: strength parameter $S$ and fracture energy G_c in the fault plane for $\mu_r$ = 0.3. In the Upper Slochteren formation in the reservoir, $S$ values are relatively low (dark blue colour) because of the stress build-up by gas production and the low cohesion strength S ₀ of the sandstone, as shown in Figs 4 and 5. Oscillations with low values of $S$ around the weak jog at a depth of 3.3 km are a numerical artifact. As the focus in this simulation was on rupture arrest at the step in the reservoir, the mesh was considerably coarsened at the depth of this jog and deeper to reduce computation time. For $\mu_r$ = 0.25 and for the same load, $S$ values are on average $\sim$ 20% lower than for $\mu_r$ = 0.3. In the step impeding or arresting the rupture along fault strike, $S \gt $ 5. Contrary to a compressive jog, the increase in the fracture energy is modest here, i.e. $\lt $ 15%.

Results of 3D simulations

From a multitude of possible 3D fault plane irregularities and a limited number of simulations, rupture propagation and arrest on three representative generic fault geometries have been selected to show in this paper the impact of steps, kinks and fault dip transitions along fault strike and of jogs along fault dip on rupture arrest. The orientation of the fault planes is representative for the major NW–SE faults in the central region. The fault dip is 70^∘, and the fault strike azimuth is 120^∘. To show the impact of stress drop or strength parameter S on rupture propagation and arrest, the residual high-velocity friction coefficient μ _r in the reservoir and Carboniferous has been varied in the range 0.2–0.3.

Without steps along fault strike, the rupture would continue to propagate along fault strike in the high-porosity Upper Slochteren sandstone unless the fault dip significantly changes. It is not stopped by representative porosity variations in the reservoir or by kinks along fault strike, at least up to changes in fault strike azimuth of 45^∘; an example is shown in Fig. 11. As for the 2D simulations, downwards propagation of the rupture can be impeded or arrested by compressive jogs but a rupture can circumvent the jog if it is of limited size along fault strike as elucidated this figure. Relative to compressive jogs, steps along fault strike appear to be less effective to arrest ruptures. The effective normal stress due to the field and induced stresses in the plane of the step is less, especially if the step plane is nearly vertical.

Figure 11. Top: three snapshots of slip in the fault plane with two 30 $^\circ$ opposing kinks along fault strike, 150 m apart and with a local jog along fault dip at 3.1 km depth. The kinks, generated by 2D parametric sigmoid functions, are practically instantaneous. The strength of cohesion S ₀ around the jog is in this simulation only $\sim$ 3 MPa. The snapshots are taken 0.14, 0.16 and 0.24 s after nucleation at 3.0 km depth and 50 m left from one of the kinks. In the last snapshot, the rupture has generated a earthquake of magnitude M $\sim$ 3.2. Before reaching the end of the jog, the rupture does not pass the jog. It circumvents it with minor slip on the jog hereafter. The two kinks along fault strike form no barriers arresting the rupture. Not shown, this also holds for 45 $^\circ$ kinks. Bottom: Strength parameter $S$ (left), stress drop $\Delta \tau$ (centre) just before nucleation and rake angle $\lambda$ (right) during rupture. The nucleation zone is indicated by the dark blue spots in the left and centre figures. On the plane of the jog $S \gt $ 5. In the Upper Slochteren in plane B, $S$ is lower than in the other planes because of a local high porosity streak and herewith a lower cohesion strength S ₀. The shadow in the figures to elucidate the fault geometry might give a wrong suggestion that this also holds for whole plane B. $\lambda$ is in the range − 100 to − 107° in most of plane A $_1$ , − 99 to − 102° in most of planes A $_2$ and C and − 96 to − 101° in most of plane B. Where the kinks approach the Zechstein, deviations of $\lambda$ from normal faulting are larger.

The simplest 3D geometry shown here is a fault with a single step along fault strike of 10 m width as shown in Fig. 6. The nucleation zone at 3 km depth is 100 m away from the step. Figure 10 shows a few snapshots of the slip on the fault after nucleation for a few of these simulations and the strength parameter S and fracture energy G _c over the fault plane just before rupture. The strength parameter S substantially increases at the step primarily because of reduced induced stress along fault dip. Relative to the compressive jog, the effective normal stress is low and the increase of the fracture energy G _c at the step is minor. Again, the stress drop, affecting the value of S, is important for rupture arrest. For μ _r = 0.3, the rupture is arrested at the step but for μ _r = 0.25 or lower, the rupture is only impeded by the step and continues to propagate along fault strike.

As for 2D simulations, the rupture velocity V _r is considerably lower than the shear velocity V _s . It varies in direction and in the fault plane depending on fault strength and stress on the fault. For μ _r = 0.3, along fault strike in the high-porosity Upper Slochteren, V _r ∼ 1.1 km/s in reasonable agreement with a value derived for the Zeerijp M_L 3.4 earthquake in 2018. Along fault dip in the Lower Slochteren, V _r ∼ 0.75 km/s and towards the Zechstein, V _r ∼ 0.9 km/s. As a result, the rupture plane is non-circular. The velocity along fault dip is considerably lower than for the 2D simulations. One reason is that dynamic stress transfer to the curved rim just in front of the rupture is less than for a straight rim of a 2D rupture. For μ _r = 0.2, V _r increases with ∼20%.

Although downwards rupture propagation can be impeded or arrested by jogs, a rupture can circumvent the jog if the jog is of limited size along fault strike and the strength of cohesion is low, see Fig. 11. The jog starts in the centre of the fault plane shown and develops in positive X-direction. The snapshots illustrate that the rupture is not arrested by the kinks of 30^∘ along fault strike and circumvents the jog when the fault dip remains constant. This also holds for simulations with 45^∘ kinks along fault strike (not shown). The figure also shows the strength parameter and stress drop just before nucleation and the rake angle λ during rupture. λ varies over the fault plane, mostly in the range − 96 to − 107^∘. Where kinks approach the Zechstein, local deviations of the slip direction from normal faulting increase. So far, this does not explain that a few observed rake angles strongly deviate from normal faulting. Such rake angles can only be expected for faults with large fault dip angles and a higher horizontal field stress anisotropy, i.e. σ _H/σ _h> 1.1.

According to Fig. 1, the dip angle of the faults in the reservoir considerably varies along fault strike in the central part of the Groningen field. For this reason, we show as an example, a rupture on a part of the fault plane with a fault dip transition from a dip of 70^∘ to a dip of 80^∘ along fault strike. Faults with dip angle variations in this range are not uncommon in the field. The transition is over a length of about 100 m in the reservoir. Figure 12 shows the variation of the strength parameter S over the fault and several snapshots of the fault slip after nucleation. The rupture does not propagate into the left part of the fault plane with a dip angle of 80^∘ because of the considerably lower load on this part of the fault as expressed by the much higher value of the strength parameter S.

Figure 12. Fault plane with a gradual transition of the fault dip from 70 $^\circ$ on the right side to 80 $^\circ$ on the left side of the transition over a length of about 100 m in the reservoir along fault strike. Left: The strength parameter $S$ just before nucleation with higher values on the left side where the loading on the fault is less. The nucleation starts in weak rock in the transition of the fault dip in the reservoir indicated by the dark blue patch. Centre and Right: two snapshots of slip in the fault plane which clearly show that the rupture does not penetrate into the part with the higher fault dip angle on the left side.

Discussion

For natural earthquakes geometrical irregularities such as jogs and steps have long be recognised as local structural controls to stop or impede ruptures, leading to the radiation of high-frequency energy and strong ground motions. A recent statistical analysis on 29 strike-slip ruptures shows that the rupture magnitude correlates well with the type and width of the jog (or step) but not with their number, see Li and Zhou (2018). Jogs were introduced in 2D dynamic rupture models for large earthquakes in the 1990s, see e.g. Harris et al. (1991). In recent years, jogs and kinks have been included in dynamic rupture models by, e.g. Yikilmaz et al. (2015) and Hisakawa et al. (2020).

So far, fault plane irregularities have not been addressed in models for induced seismicity by gas production and rupture arrest in the field has only been addressed for idealised 2D planar faults by Wentinck (2018a), Buijze et al. (2019) and Buijze (2020). In the Groningen field, the mean dip angle of the faults in the reservoir varies along fault strike, as shown in Fig. 1 and sandstone outcrops show kinks, jogs and steps along fault planes. Such features are usually related to transtensional faulting. Although possible jogs and steps indicated by the seismic attribute Ant-tracking in the NW–SE fault are close to or beyond seismic resolution and the observed roughness may result from seismic noise, fault roughness on the Groningen faults in the form of jogs and steps is not unlikely.

To investigate the impact of these geometrical irregularities on rupture arrest in the Groningen gas field, we have performed a small series of dynamic rupture simulations. This also enabled us to look to effects of horizontal field stress anisotropy and pressure diffusion in the Carboniferous on rupture propagation. Because of the data richness for the field, we believe that the input parameters for the dynamic rupture models could be reasonably well constrained. Although there are no indications that penetration of ruptures into the Carbonifeorus has happened so far, the simulations include a low cohesion strength in the Carboniferous NW–SE faults. Low values could be the result of N-S tectonic compression in the Quaternary if these faults would have functioned as compressional-stepovers, and these faults were reactivated. The clay-rich parts of these faults could have been poorly healed. Conservatively, we have disregarded a possible viscous relaxation of the field stress in the Carboniferous over geological time. If the minimal horizontal field stress in the Carboniferous would be ∼5% or 2–3 MPa higher than estimated from well data from the reservoir as a result of this relaxation, rupture runaway into the Carboniferous would be considerably less likely. A minifrac test in the Carboniferous would resolve this.

Jogs, steps and fault dip variations are plausible mechanisms to explain relatively few M_L> 3 earthquakes in the field despite the expected lateral continuity of induced stress along the faults and the observed uni-directional propagation of some ruptures along fault strike under the constrained input parameters. The pressure drop diffusion into the Carboniferous is important to include in the simulations. It stabilises the fault because of a higher effective normal stress and on the other hand it reduces the induced stress by gas production at the base of the reservoir, important for faults with considerable throw.

Since a substantial rupture penetration into the Carboniferous may lead to much stronger earthquakes in the Groningen field than observed so far, the questions stand out whether the horizontal field stress in the Carboniferous is somewhat higher than assumed and whether there are sufficient jogs and steps or fault dip changes in the faults in the field to prevent this penetration. From sandstone outcrops and Ant-tracking of two major NW–SE faults in the field, we believe that several jogs per km fault with widths of a few tenths of metres are not unlikely. Future studies on jog and step densities, and methods to derive them either from seismic attributes and/or from sandstone outcrops and on possible correlations between their appearance and substantial changes in the lithology, such as at the base of the reservoir, are recommended.

Conclusions

Fault planes in the Groningen field are irregular and show fault dip variations along fault strike. Presumably, the fault planes contain kinks, jogs and steps. Such features can be frequently seen in outcrops of sandstone formations, e.g. in the western part of the US. Ant-tracking of two representative major NW–SE faults in the central region of the Groningen field suggests the same albeit that the irregularities found are close to seismic resolution. This work is a first attempt to show the impact of fault plane irregularities on rupture arrest in the reservoir from dynamic rupture simulations. It is found that jogs and steps of tenths of metres in the fault plane and fault dip transitions can stop ruptures in the depleted gas reservoir. This holds for rupture propagation along fault dip and for rupture propagation along fault strike. In both cases, the strength parameter S increases around these irregularities. The simulations indicate that a considerable stress drop can only be realised if the high-velocity residual friction coefficient is in the range 0.2–0.3, as observed in the laboratory.

If parts of the major NW–SE faults were reactivated in the last millions of years, it cannot be excluded that the clay-rich Ten Boer and Carboniferous have a relatively low cohesion strength in the fault zone in these formations because of poor healing. If so, analytical models for rupture arrest indicate that M_L > 2 earthquakes on poorly healed planar faults and with fault dips equal or less than 70–75^∘ had potential to propagate into the Carboniferous in the last decades. But, so far there are no indications that ruptures have propagated substantially into the Carboniferous.

This suggests that apart from a possible reduction of the induced stress at the base of the reservoir, fault zones in the Carboniferous have sufficient cohesion strength and/or the horizontal stress in the Carboniferous is higher than assumed and/or, as highlighted by the simulations in this work, jogs and steps and fault dip transitions have contributed to stop the ruptures before entering the Carboniferous.