Space-Time Adaptive Methods for Subsurface Flow Simulations – STEERS

In numerical simulations, the result computed by the computer differs from reality; this discrepancy is known as numerical error. It arises because a continuous phenomenon (for example, the underground movement of CO2) is replaced by calculations performed on a spatial grid with discrete time steps. If this error becomes too large, predictions (such as potential leaks) become inaccurate. Such errors can be estimated using error indicators, which rely on the already computed solution. Mathematical theory ensures that the error, i.e. the difference between the exact solution and the computed solution, is bounded from above by the indicator. This indicator then guides the automatic refinement or coarsening of the grid or time step, thereby maintaining simulation reliability while limiting computational cost.

Constructing an error indicator generally involves three key components: selecting a measure of the error (typically a norm in a Hilbert space), deriving an upper bound for this error using computable quantities, and establishing a lower bound. The upper bound is essential to guarantee that the simulation error remains controlled. The lower bound ensures that the error indicator and the actual error are of the same order of magnitude.

To reduce computational costs, the STEERS project proposes:

1. The development of a high-order hybrid discontinuous Galerkin method capable of handling polytopal meshes while preserving mass conservation;

2. The design and analysis of space-time adaptive algorithms driven by error indicators to refine or coarsen the mesh as needed;

3. The implementation of these methods within an open-source software library.

The approach and developments are validated on problems of increasing complexity: starting with linear models (such as Darcy flow), progressing to nonlinear models (multiphase flow), and finally addressing domains that explicitly account for the presence of fractures.

The objective of the project's first task is to simulate advection-diffusion problems in fractured porous media, guiding mesh refinement and time-step selection using error indicators. This ambitious goal is being achieved through the following steps:

1. Extension of Error Indicators to Stationary Advection-Diffusion on Polytopal Meshes:

We have extended error indicators to the stationary advection-diffusion model featuring heterogeneous diffusion coefficients and velocity fields on polytopal meshes (polygons, polyhedra). The model is discretized using a Discontinuous Galerkin (DG) method, which naturally accommodates polytopal elements. Our primary achievement was extending error indicators previously proposed and analyzed in the literature for simplicial meshes (triangles, tetrahedra) to the more general polytopal framework. This required combining various recent mathematical tools. In parallel, a Julia code leveraging the Gridap library is under development. This code implements various high-order hybrid methods, specifically Hybrid High-Order (HHO) and Discontinuous Galerkin (DG), along with their associated error indicators for stationary diffusion models. The incorporation of advection is currently underway. A paper presenting the analysis of these indicators and their performance is being drafted.

2. Fracture Handling and Hybrid DG-HHO Methods:

We are integrating fractures into the adaptation procedure by implementing a method that combines DG discretization within the interior of fracture elements and HHO on the edges representing intersections between fractures. This approach significantly reduces the number of unknowns compared to applying HHO across all elements. To achieve the optimal balance between accuracy and computational cost, the initial triangular mesh of the fracture network undergoes an agglomeration procedure to generate a new mesh composed of a small number of polygons. The adaptation procedure then refines this polygonal mesh until the desired precision level is reached. Preliminary results obtained on moderately sized fracture networks show very promising prospects for reducing computational costs. A paper on this work is in progress, and the findings will be presented at the CMWR 2026 conference.

3. Space-Time Adaptation for Wave Equations:

The study of a space-time adaptation algorithm for the wave equation has enabled the initial incorporation of the time dimension into our simulations. This adaptation algorithm is generic (independent of the specific physical model) and can therefore be reused for ongoing work on flow problems. To handle mesh changes, a projection step of discrete quantities between successive meshes may be required. This work has resulted in two publications.

The next phase of the project aims to complete the work on advection-diffusion models through the following steps: extending error indicators from the stationary case to the time-dependent (unsteady) advection-diffusion problem; incorporating the advection and time-dependent components into the Julia code; validating the space-time adaptation algorithm for both HHO and DG methods; accounting for fractured porous media; and transitioning from the current Julia prototype to a fully-fledged open-source software library.

In the longer term, the project has two additional objectives: extending these methods to handle nonlinear diffusion models and multiphase flow models.

Developing and improving CO2 capture and storage (CCS) technologies are essential steps for meeting the global net zero emissions expected by 2050 (Paris agreement, adopted in 2015). The efficiency and safety of such new technologies could be hindered by the presence of fractures. Fracture are among of the primary risk of CO2 leakage. Their structures and connections play a critical role in the dynamics of governing subsurface processes. Numerical simulation is a key tool in order to perform physical process forecasting and risk assessment studies of these new technologies. Building trustworthy predictions require models that accurately describe the subsurface and the various processes over large space time scales.
Such simulations are not easily tractable as their computational cost would be prohibitive with classical numerical methods based on quasi uniform meshes and uniform time steps. Therefore, we propose a novel approach based on space time adaptive numerical methods. For a given level of accuracy, STEERS will reduce the computational cost thanks to three main ingredients: 1/ a combined Hybrid High-Order / Discontinuous Galerkin (HHO/DG) method for the spatial discretization on agglomerated meshes; 2/ a posteriori error estimates steered space time mesh adaptivity algorithms; 3/ an implementation of the proposed space time algorithms in an open source parallel library. This library will be validated throughout the project on applications of increasing complexity: from linear problems (Darcy type) to non linear ones (multiphase flow), from small scale problem to large scale ones (up to thousands of fractures).