Inga Berre

Eirik Keilegavlen, Alessio Fumagalli, Runar Berge et al.

Fractures are ubiquitous in the subsurface and strongly affect flow and deformation. The physical shape of the fractures, they are long and thin objects, puts strong limitations on how the effect of this dynamics can be incorporated into standard reservoir simulation tools. This paper reports the development of an open-source software framework, termed PorePy, which is aimed at simulation of flow and transport in three-dimensional fractured reservoirs, as well as deformation of the reservoir due to shearing along fracture and fault planes. Starting from a description of fractures as polygons embedded in a 3D domain, PorePy provides semi-automatic gridding to construct a discrete-fracture-matrix model, which forms the basis for subsequent simulations. PorePy allows for flow and transport in all lower-dimensional objects, including planes (2D) representing fractures, and lines (1D) and points (0D), representing fracture intersections. Interaction between processes in neighboring domains of different dimension is implemented as a sequence of couplings of objects one dimension apart. This readily allows for handling of complex fracture geometries compared to capabilities of existing software. In addition to flow and transport, PorePy provides models for rock mechanics, poro-elasticity and coupling with fracture deformation models. The software is fully open, and can serve as a framework for transparency and reproducibility of simulations. We describe the design principles of PorePy from a user perspective, with focus on possibilities within gridding, covered physical processes and available discretizations. The power of the framework is illustrated with two sets of simulations; involving respectively coupled flow and transport in a fractured porous medium, and low-pressure stimulation of a geothermal reservoir.

Bernd Flemisch, Inga Berre, Wietse Boon et al.

This paper presents several test cases intended to be benchmarks for numerical schemes for single-phase fluid flow in fractured porous media. A number of solution strategies are compared, including a vertex and a cell-centered finite volume method, a non-conforming embedded discrete fracture model, a primal and a dual extended finite element formulation, and a mortar discrete fracture model. The proposed benchmarks test the schemes by increasing the difficulties in terms of network geometry, e.g. intersecting fractures, and physical parameters, e.g. low and high fracture-matrix permeability ratio as well as heterogeneous fracture permeabilities. For each problem, the results presented by the participants are the number of unknowns, the approximation errors in the porous matrix and in the fractures with respect to a reference solution, and the sparsity and condition number of the discretized linear system. All data and meshes used in this study are publicly available for further comparisons.

Mats Kirkesæther Brun, Thomas Wick, Inga Berre et al.

This paper concerns the analysis and implementation of a novel iterative staggered scheme for quasi-static brittle fracture propagation models, where the fracture evolution is tracked by a phase field variable. The model we consider is a two-field variational inequality system, with the phase field function and the elastic displacements of the solid material as independent variables. Using a penalization strategy, this variational inequality system is transformed into a variational equality system, which is the formulation we take as the starting point for our algorithmic developments. The proposed scheme involves a partitioning of this model into two subproblems; phase field and mechanics, with added stabilization terms to both subproblems for improved efficiency and robustness. We analyze the convergence of the proposed scheme using a fixed point argument, and find that under a natural condition, the elastic mechanical energy remains bounded, and, if the diffusive zone around crack surfaces is sufficiently thick, monotonic convergence is achieved. Finally, the proposed scheme is validated numerically with several bench-mark problems.

Mats Kirkesæther Brun, Elyes Ahmed, Inga Berre et al.

This paper concerns splitting-based iterative procedures for the coupled nonlinear thermo-poroelasticity model problem. The thermo-poroelastic model problem we consider is formulated as a three-field system of PDE's, consisting of an energy balance equation, a mass balance equation and a momentum balance equation, where the primary variables are temperature, fluid pressure, and elastic displacement. Due to the presence of a nonlinear convective transport term in the energy balance equation, it is convenient to have access to both the pressure and temperature gradients. Hence, we introduce these as two additional variables and extend the original three-field model to a five-field model. For the numerical solution of this five-field formulation, we compare three approaches that differ by how we treat the coupling/decoupling between the flow and/from heat and/from mechanics; these approaches have in common a simultaneous application of the fixed-stress splitting scheme on both the non-linearity and the coupling structure of the problem. More precisely, the derived procedures transform a nonlinear and fully coupled problem into a set of simpler subproblems to be solved sequentially in an iterative fashion. We provide a convergence proof for the derived algorithms, and validate our results through several numerical examples.

Eren Ucar, Eirik Keilegavlen, Inga Berre et al.

Simulating the deformation of fractured media requires the coupling of different models for the deformation of fractures and the formation surrounding them. We consider a cell-centered finite-volume approach, termed the multipoint stress approximation (MPSA) method, which is developed in order to discretize coupled flow and mechanical deformation in the subsurface. Within the MPSA framework, we consider fractures as co-dimension one inclusions in the domain, with the fracture surfaces represented as line pairs in 2D (faces in 3D) that displace relative to each other. Fracture deformation is coupled to that of the surrounding domain through internal boundary conditions. This approach is natural within the finite-volume framework, where tractions are defined on surfaces of the grid. The MPSA method is capable of modeling deformation considering open and closed fractures with complex and nonlinear relationships governing the displacements and tractions at the fracture surfaces. We validate our proposed approach using both problems for which analytical solutions are available and more complex benchmark problems, including comparison with a finite-element discretization.

Antoine Tambue, Inga Berre, Jan M. Nordbotten

Simulation of geothermal systems is challenging due to coupled physical processes in highly heterogeneous media. Combining the exponential Rosenbrock--Euler and Rosenbrock-type methods with control-volume (two-point flux approximation) space discretizations leads to efficient numerical techniques for simulating geothermal systems. In terms of efficiency and accuracy, the exponential Rosenbrock--Euler time integrator has advantages over standard time-dicretization schemes, which suffer from time-step restrictions or excessive numerical diffusion when advection processes are dominating. Based on linearization of the equation at each time step, we make use of matrix exponentials of the Jacobian from the spatial discretization, which provide the exact solution in time for the linearized equations. This is at the expense of computing the matrix exponentials of the stiff Jacobian matrix, together with propagating a linearized system. However, using a Krylov subspace or Leja points techniques make these computations efficient. The Rosenbrock-type methods use the appropriate rational functions of the Jacobian from the spatial discretization. The parameters in these schemes are found in consistency with the required order of convergence in time. As a result, these schemes are A-stable and only a few linear systems are solved at each time step. The efficiency of the methods compared to standard time-discretization techniques are demonstrated in numerical examples.

Eren Ucar, Inga Berre, Eirik Keilegavlen

Shear dilation based hydraulic stimulations enable exploitation of geothermal energy from reservoirs with inadequate initial permeability. While contributing to enhancing the reservoir's permeability, hydraulic stimulation processes may lead to undesired seismic activity. Here, we present a three dimensional numerical model aiming to increase understanding of this mechanism and its consequences. The fractured reservoir is modeled as a network of explicitly represented large scale fractures immersed in a permeable rock matrix. The numerical formulation is constructed by coupling three physical processes: fluid flow, fracture deformation, and rock matrix deformation. For flow simulations, the discrete fracture matrix model is used, which allows the fluid transport from high permeable conductive fractures to the rock matrix and vice versa. The mechanical behavior of the fractures is modeled using a hyperbolic model with reversible and irreversible deformations. Linear elasticity is assumed for the mechanical deformation and stress alteration of the rock matrix. Fractures are modeled as lower dimensional surfaces embodied in the domain, subjected to specific governing equations for their deformation along the tangential and normal directions. Both the fluid flow and momentum balance equations are approximated by finite volume discretizations. The new numerical model is demonstrated considering a three dimensional fractured formation with a network of 20 explicitly represented fractures. The effects of fluid exchange between fractures and rock matrix on the permeability evolution and the generated seismicity are examined for test cases resembling realistic reservoir conditions.

Ivar Stefansson, Inga Berre, Eirik Keilegavlen

Over the last decade, finite volume discretizations for flow in porous media have been extended to handle situations where fractures dominate the flow. These discretizations have successfully been combined with the discrete fracture-matrix models to yield mass conservative methods capable of explicitly incorporating the impact of fractures and their geometry. When combined with a hybrid-dimensional formulation, two central concerns are the restrictions arising from small cell sizes at fracture intersections and the coupling between fractures and matrix. Focusing on these aspects, we demonstrate how finite volume methods effectively can be extended to handle fractures, providing generalizations of previous work. We address the finite volume methods applying a general hierarchical formulation, facilitating implementation with extensive code reuse and providing a natural framework for coupling of different subdomains. Furthermore, we demonstrate how a Schur complement technique may be used to obtain a robust and versatile method for fracture intersection cell elimination. We investigate the accuracy of the proposed elimination method through a series of numerical simulations in 3D and 2D. The simulations, performed on fractured domains containing permeability heterogeneity and anisotropy, also demonstrate the flexibility of the hierarchical framework.

Juan Michael Sargado, Eirik Keilegavlen, Inga Berre et al.

Numerical simulations of brittle fracture using phase-field approaches often employ a discrete approximation framework that applies the same order of interpolation for the displacement and phase-field variables. Most common is to use linear finite elements to discretize the linear momentum and phase-field equations. However the use of $P_1$ Lagrange shape functions to model the phase-field is not optimal, since the latter develops cusps for fully developed cracks that in turn occur at locations correspoding to Gauss points of the associated FE model for the mechanics. Such feature is challenging to reproduce accurately with low order elements, and consequently element sizes must be made very small relative to the phase-field regularization parameter in order to achieve convergence of results with respect to the mesh. In this paper, we combine the standard $P_1$ FE discretization of stress equilibrium with a cell-centered finite volume approximation of the phase-field evolution equation based on the two-point flux approximation that is constructed on the same simplex mesh. Compared to a pure FE formulation utilizing linear elements, the proposed framework results in looser restrictions on mesh refinement with respect to the phase-field length scale. Furthermore, initialization of the history field is straightforward and accomplished through a local procedure. The ability to employ a coarser mesh relative to the traditional implementation is shown for several numerical examples, demonstrating savings in computational cost on the order of 50 to 80 percent for the studied cases.

Inga Berre, Wietse Boon, Bernd Flemisch et al.

This call for participation proposes four benchmark tests to verify and compare numerical schemes to solve single-phase flow in fractured porous media. With this, the two-dimensional suite of benchmark tests presented by Flemisch et al. 2018 is extended to include three-dimensional problems. Moreover, transport simulations are included as a means to compare discretization methods for flow. With this publication, we invite researchers to contribute to the study by providing results to the test cases based on their applied discretization methods.

Eren Ucar, Inga Berre, Eirik Keilegavlen

Understanding the controlling mechanisms underlying injection-induced seismicity is important for optimizing reservoir productivity and addressing seismicity-related concerns related to hydraulic stimulation in Enhanced Geothermal Systems. Hydraulic stimulation enhances permeability through elevated pressures, which cause normal deformations, and the shear slip of pre-existing fractures. Previous experiments indicate that fracture deformation in the normal direction reverses as the pressure decreases, e.g., at the end of stimulation. We hypothesize that this normal closure of fractures enhances pressure propagation away from the injection region and significantly increases the potential for post-injection seismicity. To test this hypothesis, hydraulic stimulation is modeled by numerically coupling fracture deformation, pressure diffusion and stress alterations for a synthetic geothermal reservoir in which the flow and mechanics are strongly affected by a complex three-dimensional fracture network. The role of the normal closure of fractures is verified by comparing simulations conducted with and without the normal closure effect.