Safe and sustainable exploitation of geo-energy resources requires not only a comprehensive evaluation of the performance and economics of the corresponding projects but also the assessments of the associated risks, including the risk of induced seismicity. Indeed, seismic events may arise from the reactivation of natural faults and fractures due to subsurface engineering activities. Numerous anthropogenic activities including geothermal energy production and CO2 geological storage have been identified as potential triggers of these seismic events. The risks associated with induced seismicity stem from the potential for surface movement, structural damage, and negative impacts on both the environment and human health. The corresponding risk assessments highly rely on geomechanical modeling which is progressively being integrated into the reservoir modeling workflows. This integration demands high levels of integrability, flexibility and performance from the computational engines employed. These requirements, along with the complexities of the underlying physical, numerical and implementation aspects severely constrain the availability of suitable computational capabilities. Increasing societal concerns about induced seismicity amplify the demand for such capabilities, highlighting a lack of efficient and comprehensive solutions in both academia and industry. This thesis contributes to bridging this gap through the development of an innovative modeling framework. Leveraging the ubiquity of Finite Volume Methods (FVM) in traditional reservoir simulations, the newly proposed FVM schemes for coupled fluid mass and momentum balance equations present an opportunity for seamless integration of geomechanical modeling into existing reservoir modeling frameworks. As a result, the proposed approach satisfies the aforementioned requirements and presents an accurate and efficient framework for the investigation of induced seismicity in geo-energy applications. The core innovation of the thesis is represented by fully implicit schemes of FVM for the coupled modeling of faulted poroelastic media implemented in the opensource Delft Advanced Research Terra Simulator (DARTS). The schemes are based on coupled multi-point flux and multi-point stress approximations derived from the local conservation of fluid mass and forces. They support arbitrary material heterogeneity, anisotropy, boundary conditions, fluid properties, and friction laws. To further improve the performance of coupled modeling, block-partitioned preconditioning strategy has been implemented. Besides, first-of-its-kind nonlinear scheme of FVM for the pure elasticity problem has been proposed and implemented in DARTS...@en

Pore pressure fluctuation in subsurface reservoirs and its resulting mechanical response can cause fault reactivation. Numerical simulation of such induced seismicity is important to develop reliable seismic hazard and risk assessments. However, modeling of fault reactivation is quite challenging, especially in the case of displaced faults, i.e., faults with non-zero offset. In this paper, we perform a systematic benchmarking study to validate two recently developed numerical methods for fault slip simulation. Reference solutions are based on a semi-analytical approach that makes use of inclusion theory and Cauchy-type singular integral equations. The two numerical methods both use finite volume discretizations, but they employ different approaches to represent faults. One of them uses a conformal discrete fault model (DFM) while the other employs an embedded (non-conformal) fault model. The semi-analytical test cases cover a vertical frictionless fault, and inclined displaced faults with constant friction and slip-weakening friction. It was found that both numerical methods accurately represent pre-slip stress fields caused by pore pressure changes. Moreover, they also successfully cope with a vertical frictionless fault. However, for the case with an inclined displaced fault with a constant friction coefficient, the embedded method can not converge for the post-slip phase, whereas the DFM successfully coped with both constant and slip-weakening friction coefficients. In its current implementation, the DFM is therefore the model of choice when accurate simulation of local faulted systems is required.@en

In this work, we describe our decisions made to perform the FluidFlower simulation study and discuss various aspects of the benchmark that are different from our usual subsurface simulation practice. We will discuss the impact of various modeling choices on the outcomes of the simulation models, such as gridding, discretization, and solver strategies, and the lessons learned, taking into account the different conditions of the FluidFlower study compared to conditions commonly dealt with in subsurface simulation. We will start with a brief description of the DARTS framework utilized for compositional simulation, the thermodynamic and physical modeling related to the atmospheric CO2 -brine system, and the modeling workflow used in our benchmark submission. Additionally, we describe a custom nonlinear solver developed for the atmospheric benchmark conditions to improve convergence including the linear solver strategy since our default two-stage preconditioner does not perform effectively. To make meaningful comparisons between each of the modeling choices, we define a baseline model which is a simplified version of our setup in the main FluidFlower benchmark. The baseline model is then used to study the effect of Cartesian and unstructured meshes and a two-point flux approximation compared with a multi-point flux approximation for capturing the physics at play. We conclude our work with lessons learned and future recommendations.

@en

Elliptic differential operators describe a wide range of processes in mechanics relevant to geo-energy applications. Extensively used in reservoir modeling, the Finite Volume Method with TPFA can be consistently applied to discretize only a specific type of application under severe assumptions. In this paper, we introduce a positivity preserving Nonlinear Two Point Stress Approximation (NTPSA) based on the recently developed collocated Finite Volume scheme for linear elastic mechanics. The gradient reconstruction is different from the one used in Nonlinear TPFA, but a similar form of weighting scheme is employed to reconstruct the traction vector at each interface. The convergence of the scheme is tested with a homogeneous anisotropic stiffness tensor. The motivation behind the implementation of a new discretization framework in mechanics is to develop a uniform discretization technique preserving monotonicity for generic poromechanics applications.

@en

Quantification of the poromechanical response of subsurface formations due to human-induced pore pressure fluctuations is critical for the performance and stability assessment of many geo-energy systems. In particular, natural faults in the subsurface introduce the hazard of induced seismicity. Numerical modeling of fault reactivation is challenging, while the specific details of induced stresses and fault slip in reservoirs with displaced (i.e. non-zero offset) faults may cause additional challenges depending on the type of numerical formulation employed. To facilitate the systematic development and testing of numerical tools for the simulation of induced seismicity in faulted reservoirs we developed a set of semi-analytical test problems of increasing complexity, based on inclusion theory and Cauchy singular integral equations. With these we investigate the accuracy of two recently developed Finite Volume (FV) schemes with collocated and staggered arrangements of unknowns. One of them employs a conformal discrete fault model (DFM) which can guarantee sufficient accuracy at the cost of adaptive mesh refinement but may suffer from modelling and computational challenges when addressing large-scale realistic geological configurations. The second one employs an embedded (or non-conformal) discrete fault model (EDFM) which avoids the need for excessive mesh refinement, but of which the accuracy and the range of applicability are still to be investigated. We found that both numerical schemes accurately represent the pre-slip Coulomb stresses, but show different degrees of accuracy in representing the resulting depletion-induced fault slip. The semi-analytical benchmark data are available via DOI 10.4121/22240309.

@en

We present a scalable collocated Finite Volume Method (FVM) to simulate induced seismicity as a result of pore pressure changes. A discrete system is obtained based on a fully-implicit fully-coupled description of flow, elastic deformation, and contact mechanics at fault surfaces on a flexible unstructured mesh. The cell-centered collocated scheme leads to a convenient integration of the different physical equations, as the unknowns share the same discrete locations on the mesh. Additionally, a generic multi-point flux approximation is formulated to treat heterogeneity, anisotropy, and cross-derivative terms for both flow and mechanics equations. The resulting system, though flexible and accurate, can lead to excessive computational costs for field-relevant applications. To resolve this limitation, a scalable processing algorithm is developed and presented. Several proof-of-concept numerical tests, including benchmark studies with analytical solutions, are investigated. It is found that the presented method is indeed accurate and efficient; and provides a promising framework for accurate and efficient simulation of induced seismicity in various geoscientific applications.

@en

An increasing number of geo-energy applications require the quantitative prediction of hydromechanical response in subsurface. Integration of mass, momentum, and energy conservation laws becomes essential for performance and risk analysis of enhanced geothermal systems, stability assessment of CO2 sequestration and hydrogen storage, resolving the issue of induced seismicity. The latter problem is of particular interest because it exposes safety risks to people and surface infrastructure.

Implicit coupling of conservation laws is computationally demanding and the solution procedure often uses different numerical methods for different laws that complicates simulation. Recently developed Finite Volume (FV) schemes for poromechanics present a unified approach for the modeling of conservation laws in geo-energy applications. Contact mechanics at faults requires special attention due to the inequality constraints it imposes and nonlinear friction laws that strongly affect the occurrence of seismicity.

We develop a cell-centered FV scheme for the purpose of integrated simulation in Delft Advanced Research Terra Simulator (DARTS) platform. The scheme proposes a unified numerical framework capable to resolve conservation laws in a fully implicit manner using a single collocated grid. Coupled multi-point flux and multi-point stress approximations provide mass, momentum, and heat fluxes at the faces of the computational grid. We use a conformal discrete fracture model to incorporate faults, where the multi-point approximations of fluxes respect the discontinuity in displacements. The block-partitioned preconditioner that takes the advantage of linear structure of the coupled problem is developed to facilitate the performance of the simulation.

The proposed numerical scheme are validated against analytical and numerical solutions in a number of test cases. The convergence and stability of the schemes are investigated. It is found that the developed scheme is indeed accurate, stable, and efficient. Thereafter, we demonstrate the applicability of the approach to model fault reactivation at the laboratory scale. In a core injection test, we validate the results of simulation against experimental measurements. Next, we investigate the performance of the different preconditioning strategies. The proposed block-partitioned preconditioning strategy demonstrates the scalability and efficiency of the numerical framework.@en

We develop a collocated Finite Volume Method (FVM) to study induced seismicity as a result of pore pressure fluctuations. A discrete system is obtained based on a fully-implicit coupled description of flow, elastic deformation, and contact mechanics at fault surfaces on a fully unstructured mesh. The cell-centered collocated scheme leads to convenient integration of the different physical equations, as the unknowns share the same discrete locations on the mesh. Additionally, a multi-point flux approximation is formulated in a general procedure to treat heterogeneity, anisotropy, and cross-derivative terms for both flow and mechanics equations. The resulting system, though flexible and accurate, can lead to excessive computational costs for field-relevant applications. To resolve this limitation, a scalable parallel solution algorithm is developed and presented. Several proof-of-concept numerical tests, including benchmark studies with analytical solutions, are investigated. It is found that the presented method is indeed accurate, stable and efficient; and as such promising for accurate and efficient simulation of induced seismicity.

@en