This thesis investigates the mechanism of an induced earthquake associated with subsurface reservoir depletion, focusing on numerical simulations that incorporate the real-world reservoir geometry of the Groningen gas field in the Netherlands. It begins with a review of the poroelastic theory and its relevance to stress changes induced by reservoir depletion. The study then examines how fault offset and three-dimensional structural complexities, particularly fault intersections and horst formations, influence stress localisation and fault reactivation. While earlier studies typically consider a single fault in simplified reservoir settings, this work demonstrates that accurate modelling of the full 3-D fault system is critical for capturing realistic rupture behaviour of induced earthquakes.

To quantify these effects, we perform 3-D geomechanical simulations incorporating a faulted reservoir model based on the Groningen field, including two intersecting faults and the resulting horst structure. The study specifically focuses on the 2018 M_L Zeerijp earthquake, using numerical simulations to calculate the stress evolution over the reservoir’s production history and the fault slip during the induced earthquake. Synthetic seismic data are generated and benchmarked against field observations, including event magnitude, depletion level at reactivation, waveforms and the inverted focal mechanism.

The results demonstrate that the fault intersection angle influences not only the depletion level required for reactivation, but also the location of the slip initiation and the resulting rupture pattern. In the subsequent simulation of the 2018 Zeerijp earthquake, we observed a rupture pattern consistent with that seen in the sensitivity study for a similar intersection angle. The model also reproduces similar depletion levels, local magnitude, waveform characteristics, and focal mechanisms. These results demonstrate that current poroelastic models, when combined with realistic geological and structural representations, are capable of capturing key features of induced seismicity.

We also investigate the relationship between the inferred hypocentre location and the frequency content of the input waveforms used in inversion. This analysis is based on the simulated rupture of the 2018 Zeerijp earthquake, using both synthetic and field-observed waveforms. We observe that the estimated hypocentre shifts from the centre of the slip patch to the initial slip area when higher frequency components are included. This shift is attributed to the fact that faster slip during the rupture generates higher frequency seismic waves, a behaviour previously observed in large tectonic earthquakes. Our results show that this effect is also detectable in moderate-magnitude induced events, suggesting the potential of frequency-dependent waveform analysis to resolve rupture histories and source dynamics of reservoir-depletion-induced earthquakes. ... Read more

Local slopes carry useful information about the directionality of the predominant events in a seismic dataset and therefore can be used to steer the reconstruction process of sparsely sampled data. However, in the presence of spatial aliasing (for example, in the crossline direction of streamer data), conventional algorithms fail to provide a reliable estimate of such slopes and only low-frequency, smooth versions of the slope field can be produced. We show that provided the availability of multi-component data, and more precisely the pressure wavefield and its first-order gradient, such slopes are naturally embedded in the data and can be easily obtained by smoothed division of those wavefields. We further show that the estimated slopes can be used as regularization in a multi-channel sparse interpolation problem, providing additional guidance to the reconstruction process compared just using the pressure data and its gradient at the available traces. Numerical examples on 2D and 3D datasets confirm the effectiveness of the proposed two-stage process for multi-channel seismic data reconstruction. ... Read more

To investigate the physical processes behind induced seismicities due to, for example, production of hydrocarbons from a reservoir, most of the earlier studies performed geomechanical simulations on a simple reservoir geometry. The effect of fluid depletion is, in general, simulated for such a simple geometry. Neglecting the contribution of realistic 3-D reservoir geometries can lead to a wrong estimation of the incremental stress field. A reliable estimate of the induced stress field is key to producing meaningful simulation results. We perform geomechanical simulations on a simple fault model as well as a more realistic model based on the known geological structures at the earthquake source-region in Zeerijp region, the Netherlands. Our results demonstrate that the angle of the fault intersection affects the incremental stress field, including the effective normal stress, the shear stress, and hence, the Coulomb stress and the SCU value. Our results also show a shift in the rupture pattern and the location of the maximum slip on the fault plane. We conclude that, to properly evaluate the effects of production activities and to simulate precisely the in-situ stress field and the induced seismicity, the incorporation of a realistic reservoir structure in modelling is essential. ... Read more

Intersecting faults are often ignored in the geomechanical simulation of induced seismicity. To investigate the effects of fault intersection and the resulting reservoir geometry on induced seismicity, caused, for instance, by gas extraction, we have developed 3D geomechanical models considering two intersecting normal faults and the surrounding horst structure. We simulate the stress field and the dynamic fault reactivation in a uniformly depleted reservoir. We observe that a smaller intersection angle increases the incremental Coulomb stress at the lower reservoir juxtaposition, thus changing the temporal rupture pattern of the seismic event. In our dynamic simulation, the rupture propagates from the main fault to the secondary fault. We conclude that the fault intersection has important effects on the induced seismicity and should be taken into account when evaluating the seismicity risk in a specific region. ... Read more

Geomechanical modelling is generally used to simulate the nucleation of induce d earthquakes in, for instance the Groningen gas field. We apply quasi static simulation to investigate the stress changes from gas production. When a fault reaches a critical state, dynamic simulation provides information on the dynamic rupture during ea rthqu ake nucleation and the resulting wavefield . With the use of geomechanical modelling, it is possible to investigate the effects of the model parameters, e.g., depletion pattern and friction parameters. I n the modelling, the dynamic rupture at a finite fault is simulated both in space and time. The generated seismic wavefield from such a finite source is considered to be more realistic than the resulting wavefield from a point source. T he latter is often assumed in previous studies on the inversion of in duced earthquake data in the Groningen area. To link the wavefield generated by a geomechanically simulated finite source to the field seismic data for an earlier earthquake, we apply the same full moment tensor inversion to the waveform of a finite and of a point source . The inverted moment tensor from the field seismic observation provides a constraint to our geomechanical simulation. This allows us to perform a more realistic simulation of an induced earthquake. ... Read more

Estimating earthquake parameters, including their uncertainty, requires probabilistic sampling or inversion using Bayesian algorithms. One such Bayesian algorithm known to be highly efficient is the Hamiltonian Monte Carlo (HMC) algorithm, and modifying the algorithm with an additional linearization step can further increase this efficiency. However, the modified HMC relies heavily on accurate prior information to effectively sample non-linear earthquake parameters (e.g., hypocenter and origin time). Furthermore, the ability of the modified HMC to estimate non-linear parameters diminishes with respect to the high degree of non-linearity that is inherent to some types of events, such as induced earthquakes. To address this, we adjust the modified HMC to be run in multiple stages, combined with pre-determined initial prior sets. We test this adjustment using synthetic and real data from an induced earthquake event in the Groningen gas field in the Netherlands. We start by obtaining an initial estimate of the prior information and use it to draw multiple initial prior sets. We then run the HMC for each initial prior set in multiple stages where the results from the current stage serve as the prior for the next stage. As the final step, we form the final posterior distributions by selecting results that give the best fit between the observed and modeled data. Within this approach, we estimate ten earthquake parameters those are the six components of a full moment tensor solution, the centroid (three coordinate components), and the earthquake's origin time (including the static time corrections for each recording station). After obtaining the final results, we compare our findings with those of an existing earthquake catalog and several other research results. Given the available fault map of Groningen's subsurface, we found that our results have a higher degree of correlation with respect to the major subsurface faults. ... Read more

Dynamic geomechanical modeling can generate the seismic wavefield caused by a fault rupture. In dynamic fault-rupture modeling, the source is considered to be finite, with a limited extent both in space and in time. This contrasts with the definition of a point source, which is generally assumed to explain the seismic wavefield caused by an earthquake. Most earlier seismic inversion studies, including those of the induced earthquakes caused by depletion of the Groningen gas field, were performed assuming a point source. Still, finding a point-source reference from the seismic wavefield, even when generated by finite faulting, is important in order to calibrate the geomechanical simulation with field-seismic observations. To this end, we have developed a workflow that links geomechanical forward modeling to seismic moment-tensor inversion. We have tested this workflow for the dynamic rupture considering a realistic 3D layered earth model. At first, we simulate the triggering of dynamic fault slip at the center of a fault plane. Next, we invert the seismograms recorded by receivers located on or near the surface to obtain the full moment-tensor point-source representation and the location of the earthquake. The results of inversion show similar waveforms for both the point source and the finite source. The location of the inverted point source is within 400 m from the center of the slip patch. The double-couple components of the inverted moment tensor also match with the strike and the dip of the fault plane. ... Read more

Induced seismicity from a gas-producing region such as Groningen is believed to be caused by reservoir depletion due to long-term gas production. However, because of the complexity and uncertainty regarding the underground structure and composition, it is difficult to quantify the effect on induced seismicity due to gas production. Here we use finite-element modelling to investigate the seismogenic potential of a pre-existing fault reactivated due to fluid depletion, considering different model settings. By applying quasi-static poroelastic loading representing reservoir depletion, the stress and strain fields are derived from the resulting displacement field. The equilibrium of the fault is then evaluated using either rate-and-state or slip-weakening behaviour for friction. When the critical state is reached on the fault, where the shear stress is greater than the friction, the reactivation of the fault takes place. This reactivation is simulated by using a dynamic solver to observe the propagation and the arrest of the dynamic faulting, as well as the resultant wavefield due to seismic slip. By comparing the depletion value at both aseismic and seismic ruptures, and looking at the stress distribution on the fault, the pattern of rupture nucleation, and the resulting seismic wavefield, we are able to evaluate separately the effect of different model settings, including the geometry and material property of both caprock and reservoir, reservoir depletion pattern, and the friction law. Furthermore, by combining our study with the observed seismic wavefield, it is possible to obtain useful insights on the spatial variation in the source region. ... Read more

Induced seismicity from a gas-producing region such as Groningen is believed to be caused by reservoir depletion due to long-term gas production. However, because of the complexity and uncertainty regarding the underground structure and composition, it is difficult to quantify the effect on induced seismicity due to gas production. Here we use finite-element modelling to investigate the seismogenic potential of a pre-existing fault reactivated due to fluid depletion, considering different model settings. By applying quasi-static poroelastic loading representing reservoir depletion, the stress and strain fields are derived from the resulting displacement field. The equilibrium of the fault is then evaluated using either rate-and-state or slip-weakening behaviour for friction. When the critical state is reached on the fault, where the shear stress is greater than the friction, the reactivation of the fault takes place. This reactivation is simulated by using a dynamic solver to observe the propagation and the arrest of the dynamic faulting, as well as the resultant wavefield due to seismic slip. By comparing the depletion value at both aseismic and seismic ruptures, and looking at the stress distribution on the fault, the pattern of rupture nucleation, and the resulting seismic wavefield, we are able to evaluate separately the effect of different model settings, including the geometry and material property of both caprock and reservoir, reservoir depletion pattern, and the friction law. Furthermore, by combining our study with the observed seismic wavefield, it is possible to obtain useful insights on the spatial variation in the source region. ... Read more

Modelling dynamic rupture is essential to correctly describe the process of induced seismicity. Defmod, an open-source finite-element code featuring quasi-static loading, co-seismic volumetric strain, and dynamic rupture, is used to simulate the entire chain of induced seismicity, from pressure evolution due to fluid injection and extraction, building up of stress, and nucleation of dynamic faulting, to wavefield propagation towards the surface. To study induced earthquakes caused by fluid extraction, we modelled the behaviour of a 2-D poroelastic medium including a predefined fault by assigning a fluid source, either constant or varying, in a homogeneous reservoir layer to induce a pressure-field change. For each quasi-static step, the pressure field difference generates a displacement field that in turn affects the pressure through a coupling matrix, depending on Biot's coefficient. The rate of pressure variation is subject to the fluid source as well as the material properties, e.g., porosity and fluid mobility, which affect the speed and distribution of the stress build-up on the fault and thus the pattern of rupture nucleation. In addition, we implemented a predefined pressure profile to simulate the induced rupture in case of a uniform depletion of the reservoir to allow for a comparison with other studies. The results provide useful insights on the causality between reservoir-pressure behaviour and the induced seismicity. ... Read more