Geodetic data spanning different phases of the earthquake cycle offer insights into the spatiotemporal interplay between processes driving surface deformation, such as viscoelastic relaxation, afterslip, and (re)locking. However, quantifying their contributions and explaining pre- and post-earthquake displacements with a single set of rheological parameters is challenging. We set up a 2-D earthquake cycle finite element model that simulates the mantle and a thin low-viscosity shear zone with a temperature-dependent linear Maxwell or nonlinear power-law rheology. We use the ensemble smoother with multiple data assimilation to estimate ensembles of parameters describing the rheological makeup of the subduction zone. We assimilate onshore and offshore displacement time series acquired before and after the 2011 Tohoku-Oki earthquake. Our models provide a unique, robust solution using a temperature-dependent power-law rheology. The estimated creep parameters for the mantle wedge deeper than ∼50 km and sub-slab mantle align with laboratory experiments. However, different creep parameters are required for the shallow part of the mantle wedge than the deeper part to explain the observed postseismic response—highlighting the need for shallow viscoelastic relaxation. The trade-off between water fugacity and activation energy hinders their individual estimation but yields a well-constrained viscosity structure. The spatial distribution of vertical displacements as well as the temporal signature of early postseismic horizontal displacements are required to estimate individual parameters for afterslip and viscoelastic relaxation. Afterslip occurs downdip of the coseismic rupture. Near-trench landward motion during the early postseismic period is driven by elastic stress release beneath the oceanic plate and sub-slab asthenospheric flow.

Forecasting solar radiation is critical for balancing the electricity grid due to increasing production from solar energy. To this end, we need precise simulation of clouds, which is traditionally done by numerical weather prediction. However, these large-scale (LS) models struggle especially with forecasting stratocumulus clouds because their coarse vertical resolution cannot capture the sharp inversion present at stratocumulus cloud top. To address this issue, we employ large eddy simulation (LES), which operates at high resolution and has demonstrated superior accuracy in simulating stratocumulus clouds. However, LES relies on input data from a LS model, which is imperfect. To reduce the uncertainty caused by the LS data, we integrate a single ensemble Kalman filter step at the start of simulation in the LES model, utilizing local observations. Our results show that this approach is computationally feasible, robust, and reduces prediction error at assimilation by 50%. The improvement diminishes after approximately 1 hour of simulation due to the influence of large-scale forcing. Future work will focus on enhancing the LS inflow through nested simulations with realistic lateral boundary conditions to sustain the improvements in forecasting accuracy.

Probabilistic forecasts are regarded as the highest achievable goal when predicting earthquakes, but limited information on stress, strength, and governing parameters of the seismogenic sources affects their accuracy. Ensemble data-assimilation methods, such as the Ensemble Kalman Filter (EnKF), estimate these variables by combining physics-based models and observations. While the EnKF has demonstrated potential in perfect model experiments using earthquake simulators governed by rate-and-state friction (RSF) laws, challenges arise from the non-Gaussian distribution of state variables during seismic cycle transitions. This study investigates the Adaptive Gaussian Mixture Filter (AGMF) and the Particle Flow Filter (PFF) as alternatives for improved stress and velocity estimation in earthquake sequences compared to Gaussian-based methods like the EnKF. We test the AGMF and the PFF's performance using Lorenz 96 and Burridge–Knopoff 1D models which are, respectively, standard simplified atmospheric and earthquake models. This approach, using widely recognized and commonly used testbed models in their fields, makes the methods and findings accessible to both the data assimilation and seismology communities, while supporting comparisons and collaboration. We test these models in periodic, and aperiodic conditions, and analyze the impact of assuming Gaussian priors on the estimates of the ensemble methods. The PFF demonstrated comparable performance in chaotic scenarios, yielding lower RMSE for the estimates of the Lorenz 96 models and stronger resilience to underdispersion for the Burridge–Knopoff 1D models. This is vital given the limited and sparse historical earthquake data, underscoring the PFF's potential in enhancing earthquake forecasting. These results emphasize the need for careful data assimilation method selection in seismological modeling.

Bayesian-based data assimilation methods integrate observational data into geophysical forward models to obtain the temporal evolution of an improved state vector, including its uncertainties. We explore the potential of a variant, a particle method, to estimate mechanical parameters of the overriding plate during the interseismic period. Here we assimilate vertical surface displacements into an elementary flexural model to estimate the elastic thickness of the overriding plate, and the locations and magnitudes of line loads acting on the overriding plate to produce flexure. Assimilation of synthetic observations sampled from a different forward model than is used in the particle method, reveal that synthetic seafloor data within 150 km from the trench are required to properly constrain parameters for long wavelength solutions of the upper plate (i.e. wavelength ∼500 km). Assimilation of synthetic observations sampled from the same flexural model used in the particle method shows remarkable convergence towards the true parameters with synthetic on-land data only for short to intermediate wavelength solutions (i.e. wavelengths between ∼100 and 300 km). In real-data assimilation experiments we assign representation errors due to discrepancies between our incorrect or incomplete physical model and the data. When assimilating continental data prior to the 2011 Mw Tohoku-Oki earthquake (1997-2000), an unrealistically low effective elastic plate thickness for Tohoku of ∼5-7 km is estimated. Our synthetic experiments suggest that improvements to the physical forward model, such as the inclusion of a slab, a megathrust interface and viscoelasticity of the mantle, including accurate seafloor data, and additional geodetic observations, may refine our estimates of the effective elastic plate thickness. Overall, we demonstrate the potential of using the particle method to constrain geodynamic parameters by providing constraints on parameters and corresponding uncertainty values. Using the particle method, we provide insights into the data network sensitivity and identify parameter trade-offs.

Research into land subsidence caused by groundwater withdrawal is hindered by the availability of measured heads, subsidence, and forcings. In this paper, a parsimonious, linked data-driven and physics-based approach is introduced to simulate pumping-induced subsidence; the approach is intended to be applied at observation well nests. Time series analysis using response functions is applied to simulate heads in aquifers. The heads in the clay layers are simulated with a one-dimensional diffusion model, using the heads in the aquifers as boundary conditions. Finally, simulated heads in the layers are used to model land subsidence. The developed approach is applied to the city of Bangkok, Thailand, where relatively short time series of head and subsidence measurements are available at or near 23 well nests; an estimate of basin-wide pumping is available for a longer period. Despite the data scarcity, data-driven time series models at observation wells successfully simulate groundwater dynamics in aquifers with an average root mean square error (RMSE) of 2.8 m, relative to an average total range of 21 m. Simulated subsidence matches sparse (and sometimes very noisy) land subsidence measurements reasonably well with an average RMSE of 1.6 cm/year, relative to an average total range of 5.4 cm/year. Performance is not good at eight out of 23 locations, most likely because basin-wide pumping is not representative of localized pumping. Overall, this study demonstrates the potential of a parsimonious, linked data-driven, and physics-based approach to model pumping-induced subsidence in areas with limited data.

Probabilistic forecasts are regarded as the highest achievable goal when predicting earthquakes, but limited information on stress, strength, and governing parameters of the seismogenic sources affects their accuracy. Ensemble data-assimilation methods, such as the Ensemble Kalman Filter (EnKF), estimate these variables by combining physics-based models and observations.While the EnKF has demonstrated potential in perfect model experiments using earthquake simulators governed by rate-and-state friction (RSF) laws, challenges arise from the non-Gaussian distribution of state variables during seismic cycle transitions. This study investigates the Adaptive Gaussian Mixture Filter (AGMF) and the Particle Flow Filter (PFF) as alternatives for improved stress and velocity estimation in earthquake sequences compared to Gaussian-based methods like the EnKF. We test the AGMF and the PFF's performance using Lorenz 96 and Burridge-Knopoff 1D models which are, respectively, standard simplified atmospheric and earthquake models. We test these models in periodic, and aperiodic conditions, and analyze the impact of assuming Gaussian priors on the estimates of the ensemble methods. The PFF demonstrated comparable performance in chaotic scenarios, yielding lower RMSE for the estimates of the Lorenz 96 models and stronger resilience to underdispersion for the Burridge-Knopoff 1D models. This is vital given the limited and sparse historical earthquake data, underscoring the PFF's potential in enhancing earthquake forecasting. These results emphasize the need for careful data assimilation method selection in seismological modeling

This study investigates the integration of machine learning (ML) and data assimilation (DA) techniques, focusing on implementing surrogate models for Geological Carbon Storage (GCS) projects while maintaining the high fidelity physical results in posterior states. Initially, we evaluate the surrogate modeling capability of two distinct machine learning models, Fourier Neural Operators (FNOs) and Transformer UNet (T-UNet), in the context of CO₂ injection simulations within channelized reservoirs. We introduce the Surrogate-based hybrid ESMDA (SH-ESMDA), an adaptation of the traditional Ensemble Smoother with Multiple Data Assimilation (ESMDA). This method uses FNOs and T-UNet as surrogate models and has the potential to make the standard ESMDA process at least 50% faster or more, depending on the number of assimilation steps. Additionally, we introduce Surrogate-based Hybrid RML (SH-RML), a variational data assimilation approach that relies on the randomized maximum likelihood (RML) where both the FNO and the T-UNet enable the computation of gradients for the optimization of the objective function, and a high-fidelity model is employed for the computation of the posterior states. Our comparative analyses show that SH-RML offers a better uncertainty quantification when compared to the conventional ESMDA for the case study.

Featured Application: This study emphasizes the importance of comprehensive monitoring, calibration, and optimization of storage strategies in a saline aquifer. It also highlights the need to manage geomechanical risks and uncertainties. By understanding these risks and employing suitable monitoring techniques, the integrity and safety of GCS can be ensured, contributing to the reduction of CO ₂ emissions. Geological Carbon Storage (GCS) involves storing CO ₂ emissions in geological formations, where safe containment is challenged by structural and stratigraphic trapping and caprock integrity. This study investigates flow and geomechanical responses to CO ₂ injection based on a Brazilian offshore reservoir model, highlighting the critical interplay between rock properties, injection rates, pressure changes, and ground displacements. The findings indicate centimeter-scale ground uplift and question the conventional selection of the wellhead as a monitoring site, as it might not be optimal due to the reservoir’s complexity and the nature of the injection process. This study addresses the importance of comprehensive sensitivity analyses on geomechanical properties and injection rates for advancing GCS by improving monitoring strategies and risk management. Furthermore, this study explores the geomechanical effects of modeling flow in the caprock, highlighting the role of pressure dissipation within the caprock. These insights are vital for advancing the design of monitoring strategies, enhancing the predictive accuracy of models, and effectively managing geomechanical risks, thus ensuring the success of GCS initiatives.

We seek to quantify bulk viscoelastic flow, afterslip, and locking, within a rheological framework that is consistent over the entire earthquake cycle. We address this using an ensemble smoother. We construct a 2D finite element seismic cycle model with a power-law rheology in the asthenosphere. A priori information, such as a realistic temperature field and a coseismic slip distribution, is integrated into the model. Model pre-stresses are initialized during repeated earthquake cycles wherein the accumulated slip deficit is released entirely. We tailor the last earthquake to match the observed co-seismic slip of the 2011 Tohoku earthquake. The heterogeneous rheology structure is derived from the temperature field and experimental flow laws. Additionally, we simulate afterslip using a thin viscoelastic shear zone. We focus on constraining power-law flow parameters for the asthenosphere and the shear zone. We assimilate 3D GEONET GNSS displacement time series acquired before and after the 2011 Tohoku earthquake. Power-law viscosity parameters are successfully retrieved for all domains. The data require separate viscoelastic domains in the mantle wedge above and below ~50 km depth. The sub-slab asthenosphere has viscoelastic properties that are distinctly different from the mantle wedge. The trade-off between the power-law activation energy and water fugacity hinders their individual estimation. The wedge viscosity is >10^19 Pa·s during the interseismic phase. Postseismic afterslip and bulk viscoelastic relaxation can be individually resolved from the surface deformation data. Afterslip is substantial between 40-50 km depth and extends to 80 km depth. Bulk viscoelastic relaxation in the wedge concentrates above 150 km depth with viscosities <10^18 Pa·s. Landward motion of the near-trench region occurs during the early postseismic period without the need for a separate low-viscosity channel below the slab.

The particle filter is a data assimilation method based on importance sampling for state and parameter estimation. We apply a particle filter in two different quasi-static experiments with models of subsidence caused by a compacting reservoir. The first model considers uncorrelated model state variables and observations, with observed subsidence resulting from a single source of strain. In the second model, subsidence is a summation of subsidence contributions from multiple sources which causes spatial dependencies and correlations in the observed subsidence field. Assimilating these correlated subsidence fields may trigger weight collapse. With synthetic tests, we show in a model of subsidence with 50 independent state variables and spatially correlated subsidence a minimum of 10¹³ particles are required to have information in the posterior distribution identical to that in a model with 50 independent and spatially uncorrelated observations. Spatial correlations cause an information loss which can be quantified with mutual information. We illustrate how a stronger spatial correlation results in lower information content in the posterior and we empirically derive the required ensemble size for the importance sampling to remain effective. We furthermore illustrate how this loss of information is reflected in the log likelihood, and how this depends on the number of model state variables. Based on these empirical results, we propose criteria to evaluate the required ensemble size in data assimilation of spatially correlated observation fields.

In this comprehensive study, we discuss a novel approach to enhance data assimilation and uncertainty quantification in the field of Geological Carbon Sequestration (GCS). We specifically address the complexities of channelized reservoirs, which pose significant challenges due to non-Gaussian permeability distributions and the intricate non-linear physics of CO2 injection processes. Our innovative method integrates Fourier Neural Operators (FNOs) and Transformer UNet (T-UNet) with advanced data assimilation techniques - the Surrogate-based Hybrid Ensemble Smoother with Multiple Data Assimilation (SH-ESMDA) and the Surrogate-based Hybrid Randomized Maximum Likelihood (SH-RML). These techniques make use of the very efficient computation of gradients that neural networks provide and they not only improves the speed of data processing but also enhances the accuracy of predictions in synthetic data assimilation experiments for GCS applications. A key element of our approach is the use of proxy models alongside high-fidelity simulations, ensuring the consistency and reliability of physical posterior distributions. We utilized Alluvsim for detailed geological modeling and the Delft Advanced Research Terra Simulator (DARTS) for comprehensive fluid flow simulations, providing a comprehensive understanding of reservoir dynamics. A synthetic case study on a channelized reservoir model for CO2 sequestration demonstrates the effectiveness of these methods, with improvements in predicting CO2 plume migration and pressure dynamics within the reservoir. The results of our study show that the integration of FNOs and T-UNet with SH-ESMDA and SH-RML leads to enhanced prediction capabilities, particularly in the challenging context of channelized reservoirs. The SH-ESMDA method proves to be highly efficient in speeding up the data assimilation process without compromising accuracy, while SH-RML demonstrates a more effective history matching compared to standard Ensemble Smoother with Multiple Data Assimilation (ESMDA) techniques, indicating a robust strategy for assimilating complex data. This research not only contributes to the realm of GCS but also presents a novel solution for the integration of artificial intelligence with traditional methodologies that can be applied in various fields where data assimilation and uncertainty quantification are crucial. Our study paves the way for future advancements in this domain, highlighting the potential of AI-driven techniques in enhancing data assimilation and uncertainty quantification for GCS projects.

This study presents a method to address the significant uncertainties in subsurface modeling that impact the efficiency of energy transition applications such as geothermal energy extraction and CO₂ geological sequetsration. The approach combines a physics-based geomechanical proxy model with an ensemble smoother with multiple data assimilation (ES-MDA), aimed at enhancing uncertainty quantification through the integration of vertical displacement measurements from fluid production and injection. The data from wells is limited in spatial coverage, while these measurements offer extensive spatial information, improving the understanding of subsurface behavior by reflecting changes in reservoir pressure and temperature. The open-DARTS simulator for fluid flow and a geomechanical proxy are used to perform data assimilation with ES-MDA. By generating an ensemble of model realizations with varied permeability, calculating vertical displacements at the surface, and applying ES-MDA, we effectively identify the probability distribution of the vertical displacement of the model conditioned to observed subsidence data. Entropy is used as a statistical measure to quantify the reduction of uncertainty of subsurface models based on observations. Our approach was tested on a 2D conceptual and 3D realistic datasets, demonstrating its capability to provide data assimilation. This workflow represents an advancement in subsurface modeling, supporting informed decision-making in geothermal energy production and CO₂ sequestration by offering an improved alternative for data assimilation and enhancing tools for uncertainty quantification.

Digital twins of the Earth are digital representations of the Earth system, spanning scales and domains. Their purpose is to monitor, forecast and assess the Earth system and the consequences of human interventions on the Earth system. Providing users with the capability to interact with and interrogate the system, digital twins of the Earth are decision support systems for addressing environmental challenges. By informing humans of their impact on the Earth system, digital twins aspire to promote new pathways moving forward. By answering causal queries through intervention analysis, they can enhance evidence-based policy making. Existing digital twins of the Earth are primarily technological information systems that represent the physical world. However, as the social and physical worlds are intrinsically interconnected, we argue that humans must be accounted for both within and outside digital twins of the Earth: Within twins to represent human impacts and responses that are integral to the Earth system; and outside twins to govern access and development and to guide responsible use of information acquired from twins. Incorporating human interactions in digital twins of the Earth represents a transformative frontier, promising unparalleled insights into Earth system dynamics and empower humans for action.

This paper identifies and explains particular differences and properties of adjoint-free iterative ensemble methods initially developed for parameter estimation in petroleum models. The aim is to demonstrate the methods’ potential for sequential data assimilation in coupled and multiscale unstable dynamical systems. For this study, we have introduced a new nonlinear and coupled multiscale model based on two Kuramoto–Sivashinsky equations operating on different scales where a coupling term relaxes the two model variables toward each other. This model provides a convenient testbed for studying data assimilation in highly nonlinear and coupled multiscale systems. We show that the model coupling leads to cross covariance between the two models’ variables, allowing for a combined update of both models. The measurements of one model’s variable will also influence the other and contribute to a more consistent estimate. Second, the new model allows us to examine the properties of iterative ensemble smoothers and assimilation updates over finite-length assimilation windows. We discuss the impact of varying the assimilation windows’ length relative to the model’s predictability time scale. Furthermore, we show that iterative ensemble smoothers significantly improve the solution’s accuracy compared to the standard ensemble Kalman filter update. Results and discussion provide an enhanced understanding of the ensemble methods’ potential implementation and use in operational weather- and climate-prediction systems.

The feasibility of physics-based forecasting of earthquakes depends on how well models can be calibrated to represent earthquake scenarios given uncertainties in both models and data. We investigate whether data assimilation can estimate current and future fault states, i.e., slip rate and shear stress, in the presence of a bias in the friction parameter. We perform state estimation as well as combined state-parameter estimation using a sequential-importance resampling particle filter in a zero-dimensional (0D) generalization of the Burridge–Knopoff spring–block model with rate-and-state friction. Minor changes in the friction parameter ϵ can lead to different state trajectories and earthquake characteristics. The performance of data assimilation with respect to estimating the fault state in the presence of a parameter bias in ϵ depends on the magnitude of the bias. A small parameter bias in ϵ (+3 %) can be compensated for very well using state estimation (R2 = 0.99), whereas an intermediate bias (−14 %) can only be partly compensated for using state estimation (R2 = 0.47). When increasing particle spread by accounting for model error and an additional resampling step, R2 increases to 0.61. However, when there is a large bias (−43 %) in ϵ, only state-parameter estimation can fully account for the parameter bias (R2 = 0.97). Thus, simultaneous state and parameter estimation effectively separates the error contributions from friction and shear stress to correctly estimate the current and future shear stress and slip rate. This illustrates the potential of data assimilation for the estimation of earthquake sequences and provides insight into its application in other nonlinear processes with uncertain parameters.

Different data assimilation schemes such as the ensemble Kalman filter (EnKF), ensemble smoother (ES) and ensemble smoother with multiple data assimilation (ESMDA) are implemented in a hydro-mechanical slope stability analysis. For a synthetic case, these schemes assimilate displacements at the crest and the slope to estimate strength and stiffness parameters. These estimated parameters are then used to estimate the system's state and factor of safety (FoS). The results show that EnKF provides an FoS estimation with a mean close to the truth and with the smallest standard deviation, with ESMDA using the largest amount of assimilation steps also providing a mean close to the truth but with less confidence. The ES and ESMDA with fewer assimilation steps underestimate the FoS approximation and have low confidence. Assimilating measurements over a longer period provides a more accurate parameter, state and FoS estimation. ES has the best computational performance, with ESMDA performing worse, with its performance dependent on the number of assimilation steps. The computational performance of the EnKF is better than ESMDA but around 50% worse than the ES. Non-linearity of the underlying problem is a key cause of the multi-step assimilation processes having a better performance.

Data assimilation methods have been implemented on a slope stability problem, and the performance of different constitutive models and data assimilation schemes has been investigated. In the first part, a data assimilation scheme called the ensemble Kalman filter (EnKF) is implemented using a finite element model (FEM) and its performance with different constitutive models (the Mohr-Coulomb (MC) and Hardening Soil (HS) material models) is investigated to study their effect on the parameter and the factor of safety (FoS) estimation. Measurements of horizontal displacement are assimilated. The results from a synthetic example show that the HS model can generally be used to get reliable results for parameter and FoS estimation. However, using the MC model does not always output reliable parameter and FoS estimation. In the second part, the performance of different data assimilation schemes, i.e., the EnKF and ensemble smoother with multiple data assimilation (ESMDA), is studied with the preferred constitutive material model (the HS model). The results of a synthetic case show that the EnKF results in a narrower distribution for the FoS than the ESMDA method, while the latter results in FoS estimation which is closer to the ‘truth’.

Our ability to forecast earthquakes and slow slip events is hampered by limited information on the current state of stress on faults. Ensemble data assimilation methods permit estimating the state by combining physics-based models and observations, while considering their uncertainties. We use an ensemble Kalman filter (EnKF) to estimate shear stresses, slip rates and the state θ acting on a fault point governed by rate-and-state friction embedded in a 1-D elastic medium. We test the effectiveness of data assimilation by conducting perfect model experiments. We assimilate noised shear-stress and velocity synthetic values acquired at a small distance to the fault. The assimilation of uncertain shear stress observations improves in particular the estimates of shear stress on fault segments hosting slow slip events, while assimilating observations of velocity improves their slip-rate estimation. Both types of observations help equally well to better estimate the state θ. For earthquakes, the shear stress observations improve the estimation of shear stress, slip rates and the state θ, whereas the velocity observations improve in particular the slip-rate estimation. Data assimilation significantly improves the estimates of the temporal occurrence of slow slip events and to a large extent also of earthquakes. Rapid and abrupt changes in velocity and shear stress during earthquakes lead to non-Gaussian priors for subsequent assimilation steps, which breaks the assumption of Gaussian priors of the EnKF. In spite of this, the EnKF still provides estimates that are unexpectedly close to the true evolution. In fact, the forecastability for earthquakes for the same alarm duration is very similar to slow slip events, having a very low miss rate with an alarm duration of just 10 per cent of the recurrence interval of the events. These results confirm that data assimilation is a promising approach for the combination of uncertain physics and indirect, noisy observations for the forecasting of both slow slip events and earthquakes.

Geodetic observations of vertical land motion following a megathrust earthquake are key to a better understanding of processes and parameters controlling the dynamics at subduction margins. The relative contributions of dominant drivers during the postseismic phase, such as viscoelastic relaxation, afterslip and relocking, remain difficult to estimate individually and are often derived at the end of an observation period. Data assimilation can provide the means to estimate model parameters by combining physical models with observations whilst taking their uncertainties into account. We use Bayesian inference in the form of an ensemble smoother to estimate geodynamic parameters during the postseismic phase of the megathrust earthquake cycle. The ensemble smoother uses a Monte Carlo approach to represent the probability density distribution (pdf) of model states with a finite number of realizations. Prior estimates of the imperfect physical model are combined with the likelihood of noisy observations to estimate the posterior pdf of model parameters. We discuss a synthetic data experiment where observations are sampled from a 3D finite element model with noise added to represent errors in the data. With a smoother, observations at all time steps are assimilated in one go, to ensure consistency between estimated parameters and model outputs. We assimilate vertical and horizontal surface displacements into a 2D finite element viscoelastic earthquake cycle model with a power-law rheology. We incorporate heterogeneity into the viscoelastic structure and estimate the viscosity field based on a heterogeneous temperature field and experimental flow laws. Preliminary results show that model parameters, such as the extent of the cold nose, maximum depth of afterslip and flow law parameters can be recovered remarkably well by assimilating synthetic on- and offshore surface observations. We will apply the smoother to postseismic surface displacements of the Tohoku 2011 earthquake to estimate model parameters and driving processes across spatio-temporal domains.

Gas production in the Groningen field in the Netherlands reactivates faults and induces earthquakes, triggering severe societal unrest. Laboratory experiments indicate the relevant lithologies are velocity-strengthening [1], which favors stable creeping rather than instabilities and thus in theory prohibits earthquake nucleation. This makes the observed seismicity hard to explain. However, slide-hold-slide experiments on the same rocks reveal healing: stress drop and energy release of the subsequent slip after reactivation increase as healing time increases [2]. Given that the Groningen region has been tectonically inactive for millions of years, we will show that the healing history is a key factor for the nucleation of induced earthquakes.

We simulate earthquake sequences on a normal fault governed by rate-and-state friction, crosscutting a depleting gas reservoir. We reproduced the healing behavior observed in the laboratory and extended healing time up to millions of years. Our simulations suggest induced seismicity occur after fault reactivation, no matter whether the frictional behavior is velocity-weakening (VW) or -strengthening (VS). Unlike under VW friction where larger earthquakes with longer recurrence intervals occur as production continues, no subsequent earthquakes are expected under VS friction. No theoretically predictable nucleation size is measured: nucleation zone keeps expanding simultaneously when shear stress is dropping. Rupture propagation is mostly restricted to within the reservoir. Its reduced possibility of propagating into the under- and overburden reduces energy release and thus moment magnitude. The more VS the fault is, the lower the maximum slip rate will be, if under the same healing time. In other words, for a larger (a-b), a longer healing time is required before a large earthquake can happen. Therefore, if the healing history is well constrained, we can have an estimate of the seismic hazard of the coming event before it occurs. However, in order to have an accurate hazard estimate, dynamic weakening mechanisms activated at slip rate of cm/s may not be ignored. Adding flash heating shortens the required healing time by at least eight orders of magnitude in one of our representative models. Parameter space is explored systematically to build analytical expressions of earthquake magnitude as a function of healing time and frictional parameters. This study will help to understand in which tectonic histories, induced seismicity can still occur under VS friction. This is important for forecasting induced seismicity and minimizing the hazard in similar lithologies across northwestern Europe and in subsurface production and storage sites worldwide.