Source-linked AI summary
PorePy: An Open-Source Software for Simulation of Multiphysics Processes in Fractured Porous Media
Eirik Keilegavlen, Runar Berge, Alessio Fumagalli, Michele Starnoni, Ivar Stefansson, Jhabriel Varela, Inga Berre
TL;DR
Fractured-rock simulation must handle high-aspect-ratio fractures, complex networks, and coupled processes across fracture and matrix. This paper develops and implements mixed-dimensional modeling and discretization principles in the open-source PorePy software. PorePy performs well on flow, poroelasticity, and fracture-deformation benchmarks, and demonstrates coupled flow-transport and deformation simulations.
Problem
Fractured-rock simulation frameworks must handle thin, extensive fractures, contrasting fracture and matrix properties, multiphysics coupling, and complex network intersections.
Source: §1 · pp. 1–2Method
The paper defines a mixed-dimensional framework and implements it in PorePy, decomposing geometry into matrix, fracture, and intersection objects with equations on objects and their connections.
Source: §1 · p. 3Results
PorePy performs well on flow, poroelastic deformation, and fracture-deformation benchmarks, and demonstrates fully coupled flow with nonlinear transport and fracture deformation under poromechanical loading.
Source: §6 · p. 26Takeaways & Limitations
The framework supports automatic meshing of complex fracture networks in 2D and 3D and simulations spanning flow, transport, poroelasticity, and fracture deformation.
Source: §6 · p. 26The framework permits coupling only between subdomains exactly one dimension apart, ruling out direct 3D–1D coupling.
Source: §2 · p. 6
Abstract
from arXiv · showhide
Development of models and dedicated numerical methods for dynamics in fractured rocks is an active research field, with research moving towards increasingly advanced process couplings and complex fracture networks. The inclusion of coupled processes in simulation models is challenged by the high aspect ratio of the fractures, the complex geometry of fracture networks and the crucial impact of processes that completely change characteristics on the fracture-rock interface. This paper provides a general discussion of design principles for introducing fractures in simulators, and defines a framework for integrated modeling, discretization and computer implementation. The framework is implemented in the simulation software PorePy, which can serve as a flexible prototyping tool or multiphysics problems in fractured rocks. Based on a representation of the fractures and their intersections as lower-dimensional objects, we discuss data structures for mixed-dimensional meshes, formulation of multiphysics problems and discretizations that utilize existing software. We further present the implementation of these concepts in the PorePy open-source software tool, which is aimed at coupled simulation of flow and transport in three-dimensional fractured reservoirs as well as deformation of fractures and the reservoir in general. We present validation by benchmarks for flow, poroelasticity and fracture deformation in fractured porous media. The flexibility of the framework is then illustrated by simulations of fully coupled flow and transport and of injection driven deformation of fractures. All results reported herein can be reproduced by openly available simulation scripts.
1. Introduction
Fractured-rock simulation is difficult because fractures are thin yet extensive, properties and governing physics differ across fractures and matrix, and networks have complex intersections. The paper presents PorePy, an open-source framework for mixed-dimensional multiphysics modeling, and positions it between existing modeling approaches while supporting flexible prototyping.
1. Introduction
- Fractured-rock simulation must address high fracture aspect ratios, contrasting fracture and matrix properties and physics, multiphysics coupling, and complex intersections.The paper assumes averaged fracture apertures and similar quantities are available at the discretization scale, rather than modeling fracture-wall geometry in detail.Source: §1 · pp. 1–2
- Upscaled models are efficient and support multiphysics, but their accuracy depends on how well the fractured domain resembles a continuum and on the quality of upscaling.Upscaling can be analytically tractable for simple geometries and dynamics, but extremely challenging for complex geometries and multiphysics couplings.Source: §1 · p. 2
- Discrete fracture network models explicitly represent fractures while omitting the rock matrix, offering efficient fracture-dynamics simulations but excluding processes outside the network unless coupled to a continuum model.The passage cites developed DFN models for coupled flow and transport –.Source: §1 · p. 2
- Discrete fracture matrix models combine explicit representation of important fractures with upscaled treatment of others and the host rock, balancing explicit and upscaled approaches.They can represent equations in the matrix, fractures, and intersections; fractures and intersections are often modeled as lower-dimensional objects embedded in the 3D matrix [18],.Source: §1 · p. 2
- Conforming DFM methods are the paper’s focus, while non-conforming EDFM and XFEM avoid conforming-mesh complexity but require other fracture-interaction treatments.The paper names DuMuX, CSMP, MOOSE-FALCON, OpenGeoSys, and Flow123d as simulators incorporating DFM principles –.Source: §1 · pp. 2–3
- A flexible prototyping framework is motivated by the need to study strongly coupled processes in mixed-dimensional geometries, paralleling the role of MRST in non-fractured porous media.The paper presents PorePy as both design principles for mixed-dimensional multiphysics frameworks and an open-source implementation.Source: §1 · p. 3
- PorePy combines mixed-dimensional geometry, separately defined equations for geometric objects and their connections, and reusable fixed-dimensional discretizations in an open-source framework.The framework supports linear and nonlinear problems, automatic meshing of fractured domains in 2D and 3D, and models of flow, transport, and deformation.Source: §1 · p. 3
2 Design principles for mixed-dimensional simulation tools
The framework represents fractures, intersections, and matrix as geometric objects connected through interfaces, with explicit rules for coupling and communication. It organizes meshes and data as a graph, uses projections and mortar grids for inter-object transfer, and separates subdomain discretization from interface coupling and global assembly.
2 Design principles for mixed-dimensional simulation tools
- The framework represents matrix, fractures, and intersections as objects, allowing equations within objects and on their connections; lower-dimensional intersections may be explicit or replaced by a model for fracture-fracture interaction.Explicitly representing all objects is the preferred solution because it permits flexible dynamics and dimension-independent coupling; omitting intersections can make interaction modeling challenging, especially for complex dynamics and 3D networks.Source: §2 · pp. 4–5
- Couplings are restricted to neighboring subdomains exactly one dimension apart, expressed through interface variables, and interface laws depend only on the interface and its immediate neighbors.The framework therefore excludes direct 3D–1D and same-dimensional couplings; communication between same-dimensional objects must pass through a lower- or higher-dimensional object.Source: §2 · p. 6
- The mixed-dimensional geometry is organized as a graph with subdomains as nodes and interfaces as edges, providing a location for grids, parameters, variables, equations, and discretization choices.The graph structure follows from subdomains seeing only neighboring interfaces and interfaces connecting their two neighboring subdomains.Source: §2 · p. 6
- Conforming mesh construction identifies fracture intersections, meshes objects of all dimensions, and identifies neighboring cells and faces; PorePy uses Gmsh and standard fixed-dimensional meshes for the subdomains.Automatic decisions about whether nearby objects are distinct depend on prescribed mesh size, and mesh-size tuning can be delicate; PorePy supplies minimum and target sizes, though Gmsh may override them.Source: §2 · pp. 6–7
- Separate interface meshes act as mortar grids, enabling non-matching subdomain grids and potential speedups by pairing fine fracture grids with coarser matrix grids.The mortar space must be designed carefully when combining non-matching grids, non-uniform discretizations, or different physics to avoid introducing instabilities.Source: §2 · p. 7
- Projection operators transfer variables between interfaces and neighboring subdomains, with their construction depending on whether variables are intensive or extensive.The implementation considers lowest-order projections constructed from overlapping areas between interface cells and neighboring mesh faces or cells.Source: §2 · pp. 7–8
- Subdomain equations can often be defined and discretized as fixed-dimensional problems with interface terms as boundary conditions or sources, while interface laws handle coupling and projections.This separation supports reuse of fixed-dimensional simulation code; its usefulness depends on the physical process being modeled.Source: §2 · p. 8
- The global multiphysics system has blocks for subdomains and interfaces and, within each, blocks for primary variables, a structure useful for solver design and result analysis.PorePy’s global degree-of-freedom manager tracks the block structure and numbering of individual degrees of freedom.Source: §2 · pp. 8–9
3 Modeling, discretization and implementation
This section presents modeling, discretization, and implementation of flow, transport, poroelastic fracture deformation, and frictional contact in mixed-dimensional frameworks. It describes the governing equations and coupled variables, alongside a modular grid and implementation structure that supports model variations.
3.1 Flow in fractured media
- Variational and integral flow formulations retain the structure of fixed-dimensional problems, enabling software reuse and allowing different discretizations on neighboring subdomains.From the interface perspective, each subdomain discretization acts as a Neumann-to-Dirichlet map from interface fluxes to pressures projected onto the interface.Source: §3 · pp. 11–12
- The flow implementation couples subdomain discretizations through interface laws on a graph, without direct coupling between neighboring subdomains or requiring them to use the same numerical method.PorePy separates the coupling structure from subdomain discretizations and offers RT0-P0, Mvem, Tpfa, and Mpfa schemes for flow.Source: §3 · p. 13
3.2 Fully coupled flow and transport
- In coupled flow and transport, concentration-dependent viscosity couples the equations, while Darcy flux drives advective transport; finite-volume discretization uses single-point upstreaming and MPFA diffusion.The interface advective flux uses an upstream-like operator based on the Darcy interface flux.Source: §3 · pp. 14–15
3.3 Poromechanical fracture deformation by contact mechanics
- The fracture-deformation model describes opening and slip through displacement jumps, with displacement variables on matrix-fracture interfaces rather than displacements inside fractures.The exposition considers non-intersecting fractures and restricts aperture changes to the normal displacement jump.Source: §3 · pp. 15–17
3.3.2 Implementation in mixed-dimensional framework
- The contact displacement remains stationary until tangential force exceeds frictional resistance, then slides parallel to the tangential force; contact conditions use stress without fluid pressure.The contact law distinguishes open, sticking, and slipping states.Source: §3 · p. 17
- Poroelastic contact implementation extends the Mpsa finite-volume method –, adding fracture pressure and interface flux variables while treating the matrix-fracture boundary as Neumann for flow.These flow couplings are standard across discretization schemes for single-phase flow.Source: §3 · p. 18
- Displacement is imposed through a Dirichlet boundary condition and projected between matrix and mortar grids; stress continuity and nonlinear contact conditions are projected and solved with a semi-smooth Newton method [61],.The contact relation depends on whether the fracture is open, sticking, or slipping.Source: §3 · p. 18
- The mixed-dimensional grid and modular implementation divide the complex coupled system into manageable parts and allow model variations with minimal adjustments.The system contains different governing equations across subdomains and couplings between variables on different grids.Source: §3 · p. 18
4 Validation
The validation benchmarks PorePy’s flow, poroelasticity, and fracture-deformation capabilities against established test cases. The reported results include agreement with reference solutions and convergence behavior, with scripts available to reproduce the simulations.
4 Validation
- The validation uses three established test cases spanning flow, poroelasticity, and fracture deformation; scripts for reproducing reported results are openly available.The cases probe discretization schemes, multiphysics, and time-dependent problems.Source: §4 · p. 18
4.1 Flow in 2d fractured porous media
- The flow benchmark tests intersecting and isolated conductive and blocking fractures using several homogeneous, heterogeneous, and nonconforming discretization and grid configurations.The benchmark includes Tpfa, Mpfa, RT0-P0, Mvem, and combinations of schemes across subdomains.Source: §4 · p. 19
- PorePy flow solutions agree well with the fine-grid reference; matrix pressure converges first-order across the considered methods, while fracture-subdomain convergence is sublinear.The benchmark compares homogeneous, heterogeneous, and nonconforming discretizations.Source: §4 · pp. 19–20, Figure 6 · p. 20
4.2 Mandel’s problem in poroelasticity
- The poroelastic Mandel simulation gives pressure and displacement results in good agreement with, using a 622-triangle mesh.The case uses quasi-static Biot equations, Mpsa/Mpfa for mechanics and flow, and implicit Euler time discretization.Source: §4 · pp. 20–21, Figure 8 · p. 21
4.3 Sneddon’s problem of fracture deformation
- Sneddon’s fracture-deformation study finds average first-order convergence in relative normal displacement across fracture resolutions.The study averages 140 simulations per resolution across seven fracture angles and 20 grid realizations; six refinement levels comprise 840 simulations.Source: Figure 9 · p. 21, §4 · p. 22
5 Applications: Multiphysics simulations
The applications demonstrate PorePy’s flexibility on fully coupled flow and transport with non-matching meshes, and on injection-driven fracture deformation in a poromechanical setting. Together, they illustrate simulations of complex coupled processes on mixed-dimensional fracture networks.
5.1 Fully coupled flow and transport
- The flow-and-transport example applies automatic differentiation and Newton’s method to a nonlinear problem on matching and non-matching grids in a complex fracture network.The network includes 15 fracture domains, 62 one-dimensional domains, and 9 zero-dimensional domains; transport is advection-dominated in fractures and diffusion-dominated in the matrix.Source: §5 · pp. 22–23
- Matching and non-matching meshes produce average fracture-network concentration profiles with no significant differences, and the non-matching solution shows no irregularities.At t=20, concentration has increased in considerable parts of the rock matrix through diffusion and advection from the fractures.Source: Figure 10 · p. 23, §5 · p. 23, Figure 11 · p. 24
5.2 Poroelasticity and fracture deformation
- The fracture-deformation example models three fractures under 25 days of injection and a subsequent 25-day relaxation phase in an anisotropic background stress regime.The setup illustrates interactions between processes in the rock matrix and critical events on individual fractures.Source: §5 · p. 24
- Injection produces tangential displacement jumps on all three fractures, while normal jumps appear on fractures 1 and 3 but not the more distant fracture 2.During relaxation, normal jumps vanish while tangential jumps persist due to friction; tangential jumps on fractures 2 and 3 increase slightly at shut-in.Source: §5 · p. 24, Figure 12 · p. 25, Figure 13 · p. 25
6 Conclusions
The paper presents design principles for simulation software for fractured porous media and implements them in the open-source tool PorePy. The framework uses lower-dimensional representations of fractures and intersections, and PorePy demonstrates its capabilities through benchmarks and multiphysics simulations.
6 Conclusions
- The paper presents design principles for flexible simulation software to support rapid prototyping of models and discretization methods for fractured porous media.The framework is implemented in the open-source simulation tool PorePy.Source: §6 · p. 25
- The framework represents fractures and their intersections as separate lower-dimensional objects within a Discrete Fracture Matrix model.It supports multiphysics flexibility and reuse of code written for non-fractured domains, making it suitable for extending other software to mixed-dimensional problems.Source: §6 · p. 26
- PorePy provides automatic meshing of complex two- and three-dimensional fracture networks and methods for flow, transport, poroelastic deformation, and fracture contact mechanics.The implementation performs well on flow, poroelastic deformation, and fracture deformation benchmarks, and demonstrates fully coupled flow and nonlinear transport simulations.Source: §6 · p. 26