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.

@en

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.

@en

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.

@en

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.@en

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@en

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.@en

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.

@en

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.@en

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’.@en

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.

@en

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.

@en

Naturally fractured reservoirs can pose challenges for energy operations such as hydrocarbon production, CO2 storage, and geothermal energy production. Fluid flow in these reservoirs is greatly affected by fracture properties such as orientation and aperture, whose magnitude is mainly influenced by the stresses on the reservoir rocks. Simulating fractures and their behavior tends to be computationally intensive, but recent advances in Discrete Fracture Models (DFM) have successfully overcome computational complexity and allow for the explicit inclusion of discrete fractures in reservoir simulations. However, there are still challenges in dealing with uncertainties, including fracture aperture and the effect of in-situ stresses on the fracture surface and their effect on the fluid behavior. This study explores the use of data-assimilation techniques to help quantify these uncertainties. We combine a recent implementation of DFM on the Delft Advanced Research Terra Simulator (DARTS) with both ensemble and gradientbased data-assimilation methods. Our results show that data assimilation can help to understand the dynamic behavior of fluids in fractured reservoirs. Using this technique, we obtain a more accurate representation of the stresses acting on the reservoir and how they affect the fracture aperture. This information is essential for more efficient reservoir management.@en

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.@en

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.

@en

The installed capacity of geothermal systems for direct use of heat is increasing worldwide. As their number and density is increasing, the their interaction with subsurface faults becomes more important as they could lead to safety risks from induced seismicity. Assessment and management of such risks is essential for the further development and extension of geothermal energy for heating. At the same time, the economic output of geothermal systems can be marginal and is hence often supported by subsidy schemes. A combined assessment of fault stability and economic output could help operators to balance economic and safety aspects, but this is currently not common practice. In this study we present a methodology to assess field development plans based on fault stability and Net Present Value (NPV) using reservoir simulations of a fluvial, heterogeneous sandstone representative of the majority of direct-use Dutch geothermal systems. We find that the highest friction coefficient leading to exceedance of the Mohr–Coulomb failure criteria in this sandstone is 0.17; such values could be encountered in clay-rich fault gouges. Similar or lower fault permeability compared to the reservoir results in no changes and an increase respectively of both NPV and fault stability with larger Fault-to-Well Distance (FWD). Fault permeability higher than the reservoir permeability results in a minor increase in NPV with smaller FWD. Our results demonstrate that a combined analysis of thermal, hydraulic, mechanical and economic assessment supports a responsible and viable development of geothermal resources at a large scale. The importance of a high spatial density of supporting stress data will be essential for a better understanding and quantification of economic and fault stability effects of geothermal operations.

@en

A recursive ensemble Kalman filter (EnKF) is used as the data assimilation scheme to estimate strength and stiffness parameters simultaneously for a fully coupled hydro-mechanical slope stability analysis. Two different constitutive models are used in the hydro-mechanical model: the Mohr-Coulomb (MC) model and the Hardening Soil (HS) model. The data assimilation framework allows the investigation of the effect of constitutive behaviour on its ability to estimate the factor of safety using measurements of horizontal nodal displacement at the sloping face. In a synthetic study, close-to-failure and far-from-failure cases of prior property estimations illustrate the effect of initial material property distribution with different material models. The results show that both models provide a reliable factor of safety when the distribution of prior parameters is selected close-to-failure. However, the HS model results in the improved estimation of factor of safety for the far-from-failure case while this is not the case for the MC model. In addition, for the same level of accuracy the computational effort required for the HS model is comparatively less than for the MC model.

@en

Geothermal operations are expanding and increasingly contributing to the current energy supply. Assessing the long-term operable lifetime of these projects is complicated as the reservoirs they produce from are often deep and subsurface properties are uncertain and spatially variable. The minimum lifetime of a geothermal project usually considers the heat in place in a geothermal reservoir, but not the heat flow from the confining layers into the reservoir during operation. For the economic feasibility and optimal design of a geothermal project it is the key to capture this process, as this allows more accurate prediction of the long-term extraction temperature. Previous studies are not conclusive on the contribution of vertical recharge from the confining layers. This research evaluates the contribution of recharge to the heat production for geothermal projects. In a simulation study with an idealised, homogeneous geothermal system, the reservoir thickness, well placement and production rate are varied to investigate their respective influence on thermal recharge. The results show that the recharge from the vertical confining layers may contribute considerably to the total energy output of a geothermal well. Due to recharge, geothermal systems produce more heat than the heat in place and under specific conditions the produced heat can be more than five times as large. The largest contribution of the recharge of vertically confining layers can be expected under conditions of thin reservoirs and a long interaction time between the injection water and confining layers (i.e. large well spacing, low production rate). The insights of this study may help to optimise the design and operation of individual geothermal projects for optimal utilisation of a geothermal reservoir for energy supply.@en

This work demonstrates the efficiency of using iterative ensemble smoothers to estimate the parameters of an SEIR model. We have extended a standard SEIR model with age-classes and compartments of sick, hospitalized, and dead. The data conditioned on are the daily numbers of accumulated deaths and the number of hospitalized. Also, it is possible to condition the model on the number of cases obtained from testing. We start from a wide prior distribution for the model parameters; then, the ensemble conditioning leads to a posterior ensemble of estimated parameters yielding model predictions in close agreement with the observations. The updated ensemble of model simulations has predictive capabilities and include uncertainty estimates. In particular, we estimate the effective reproductive number as a function of time, and we can assess the impact of different intervention measures. By starting from the updated set of model parameters, we can make accurate short-term predictions of the epidemic development assuming knowledge of the future effective reproductive number. Also, the model system allows for the computation of long-term scenarios of the epidemic under different assumptions. We have applied the model system on data sets from several countries, i.e., the four European countries Norway, England, The Netherlands, and France; the province of Quebec in Canada; the South American countries Argentina and Brazil; and the four US states Alabama, North Carolina, California, and New York. These countries and states all have vastly different developments of the epidemic, and we could accurately model the SARS-CoV-2 outbreak in all of them. We realize that more complex models, e.g., with regional compartments, may be desirable, and we suggest that the approach used here should be applicable also for these models.@en

This article investigates bypassing the inversion steps involved in a standard litho-type classification pipeline and performing the litho-type classification directly from imaged seismic data. We consider a set of deep learning methods that map the seismic data directly into litho-type classes, trained on two variants of synthetic seismic data: (i) one in which we image the seismic data using a local Radon transform to obtain angle gathers, (ii) and another in which we start from the subsurface-offset gathers, based on correlations over the seismic data. Our results indicate that this single-step approach provides a faster alternative to the established pipeline while being convincingly accurate. We observe that adding the background model as input to the deep network optimization is essential in correctly categorizing litho-types. Also, starting from the angle gathers obtained by imaging in the Radon domain is more informative than using the subsurface offset gathers as input.@en

An ensemble-based history-matching framework is proposed to enhance the characterization of petroleum reservoirs through the assimilation of crosswell electromagnetic (EM) data. As an advanced technology in reservoir surveillance, crosswell EM tomography can be used to estimate a cross-sectional conductivity map and associated saturation profile at an interwell scale by exploiting the sharp contrast in conductivity between hydrocarbons and saline water. Incorporating this information into reservoir simulation in combination with other available observations is expected to enhance the forecasting capability of reservoir models and to lead to better quantification of uncertainty.

The proposed approach applies ensemble-based data-assimilation methods to build a robust and flexible framework in which various sources of available measurements can be integrated. A comparative study evaluates two different implementations of the assimilation of crosswell EM data. The first approach integrates the crosswell EM field components in their original form, which entails forward simulation of the observed EM responses from the simulated reservoir state. In the second approach, formation conductivity is derived from the EM data through inversion and is subsequently assimilated into the reservoir model. An image-oriented distance parameterization of the fluid front assimilates the conductivity field in an efficient and robust manner and overcomes issues with data size, errors, and their correlation. Numerical experiments for different test cases with increasing complexity provide insight into the performance of the two proposed integration schemes. The results demonstrate the efficiency of the developed history-matching workflow and the added value of crosswell EM data in enhancing the reservoir characterization and reliability of dynamic reservoir forecasts.@en