# Towards a workflow to evaluate geological layering uncertainty on CO 2 injection simulation

Capucine Legentil, Jeanne Pellerin, Margaux Raguenel, Guillaume Caumon

► To cite this version:

Capucine Legentil, Jeanne Pellerin, Margaux Raguenel, Guillaume Caumon. Towards a workflow to evaluate geological layering uncertainty on CO 2 injection simulation. Applied Computing and Geosciences, 2023, 18, pp.100118. 10.1016/j.acags.2023.100118. hal-04106548

HAL Id: hal-04106548

https://hal.science/hal-04106548

Submitted on 28 May 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>
# Towards a workflow to evaluate geological layering uncertainty on CO₂ injection simulation

Capucine Legentil<sup>a</sup>, Jeanne Pellerin<sup>b</sup>, Margaux Raguenel<sup>b,*</sup>, Guillaume Caumon<sup>a</sup>

## ARTICLE INFO

**Keywords:**

Tetrahedral meshing

CO₂ storage

Geomodeling

Flow simulation

Uncertainties

## ABSTRACT

We propose a workflow for updating 3D geological meshed models to test different layering scenarios and to assess their impact on the simulation of CO₂ injection. This workflow operates on a tetrahedral mesh that encodes rock unit information as well as rock physical properties. The alternative layering meshes are built by modifying the input mesh and inserting a new horizon defined by a scalar field. Modifying consistently a 3D meshed model while keeping its quality is a challenge that we tackle using the advanced capabilities of MMG, an open source remeshing library. CO₂ injection is then simulated with GEOSX, an open-source, multiphysics, and multilevel simulation solver. We demonstrate this workflow for stratigraphic layering uncertainty assessment on a simple synthetic layered reservoir on the flank of a salt diapir. Comparison of simulation results is eased since modifications of the mesh are localized to the area around the inserted horizon. The consistent results highlight the role of stratigraphic unconformities for trap integrity. This work opens a promising path for developing numerical simulation of CO₂ injection on unstructured meshes by combining advanced coupled flow-geomechanical models in geological domains affected by structural uncertainties.

## 1. Introduction

Among the various approaches to mitigate climate change, underground CO₂ sequestration has received significant attention during recent years (see Ajayi et al., 2019; Benson and Cole, 2008 and references therein). The potential storage sites having sufficient capacity consist of depleted oil and gas reservoirs (e.g., Hamza et al., 2021; Oldenburg et al., 2001), and deep saline aquifers (e.g., Celia et al., 2015; Sifuentes and Blunt, 2009). In both cases, it is essential to assess the risk of leakage and migration of the CO₂ towards shallower formations. One potential leakage hazard is the integrity of seal layers that are subject to various uncertainties (Benson and Cole, 2008; Fleury et al., 2011; Li et al., 2006). For example, small-scale sedimentary and stratigraphic heterogeneities can be the source of vertical migration pathways through the cap rock. Such sedimentary heterogeneities may also either increase or decrease the injectivity in the reservoir rocks depending on the spatial connectivity and geometry of rock types (Issautier et al., 2015). Permeable faults and fractures not identified on available data may be another source of CO₂ escape pathways (Bond et al., 2007) as well as fractures induced by injection pressure. Mitigating the risk of fracturing calls for optimizing the injection rate, and therefore for the characterization of hydro-mechanical rock properties

in the reservoir and the overburden, of the in-situ stress state, of fluid pressure and temperature in the geological formations, and for advanced coupled flow-geomechanical simulation codes able to predict the fate of the injected CO₂. It has also been shown that the mineralogy can play a significant role on the fluid-rock interactions and on the trapping mechanisms (e.g., Broseta et al., 2012).

Overall, geological characterization and modeling are essential approaches to understand and to evaluate all relevant aspects of CO₂ sequestration projects. However, uncertainties always exist on geological features away from wells. When seismic data are available, such uncertainty can be significantly reduced, but is not annihilated because of limited seismic resolution (Lallier et al., 2012). Testing the impact of such geological uncertainties on the outcomes of the coupled physical simulation of CO₂ injection, is, therefore, essential to mitigate risks (Wellmann et al., 2014).

In this paper, we focus on the geometrical aspects of geological uncertainty (Wellmann and Caumon, 2018), and propose a new approach to incrementally update a geological model to assess the impact of likely but undetected geological features. As such features may display relatively complex geometry, we choose to represent the geological model as a tetrahedral mesh. Indeed, unstructured meshes offer a good

<footnote><sup>a</sup> Université de Lorraine, CNRS, GeoRessources, 54000 Nancy, France</footnote>

<footnote><sup>b</sup> TotalEnergies, CSTJF, Pau, France</footnote>

<footnote>* Corresponding author.</footnote>

<footnote>E-mail address: margaux.raguenel@totalenergies.com (M. Raguenel).</footnote>

<footnote n="1">¹ https://github.com/GEOSX/GEOSX</footnote>
geometric accuracy for representing complex geological features and their resolution can be adapted spatially. As the creation of such meshes can be challenging and time-consuming, (e.g., Zehner et al., 2015), our strategy is to use a local model updating approach (Caumon et al., 2003; Legentil et al., 2022; Suter et al., 2017). More precisely, we extend the approach of Legentil et al. to three dimensions to locally insert an interface in an existing tetrahedral mesh. The use of tetrahedral meshes also allows to run coupled flow-geomechanical physical simulation on the same representation using either iterative or implicit coupling. Several numerical discretization schemes have been proposed on tetrahedral meshes, such as control volume-finite elements methods e.g. Shao et al. (2021). In this paper, we use GEOSX¹ an open-source, multiphysics simulation tool jointly created by Lawrence Livermore National Laboratory, Stanford University, and TotalEnergies (Gross and Mazuyer, 2021) that operates on unstructured meshes.

To demonstrate our local meshed model updating methodology and how it could be used in storage site assessment, we consider the case of stratigraphic unconformities, which have received relatively limited attention to-date in the area of CO₂ sequestration. Notably, Shariatipour et al. (2016), showed the impact of the possible alteration below an erosional unconformity overlying multi-layered dipping reservoir using a corner-point grid reservoir representation. These stratigraphic unconformities can occur at multiple nested scales in the stratigraphic record (Miall, 2016; Sadler and Jerolmack, 2015), which could impact the hydrodynamic connectivity of reservoirs (Lallier et al., 2012). Considering unconformities in the evaluation of potential sequestration sites is important (Wilkinson et al., 2013). In this paper, we consider a simplified but realistic layered reservoir located on the flanks of a diapir, a setting where the interaction between sediment deposition and salt displacement can typically generate multiple stratigraphic unconformities around the diapir called halokinetic sequences (Giles and Lawton, 2002). Interpreting the number and the geometry of these unconformities can be difficult because of limited seismic bandwidth and subsalt imaging challenges (Clausolles et al., 2019; Rojo et al., 2016).

**Contributions.** In this paper we propose a workflow for updating 3D geological models to test different layering scenarios and to assess their impact on the simulation of CO₂ injection (Fig. 1). We apply this general workflow on a demonstrative synthetic model of a 3D salt diapir. From an input model, we build two alternative caprock layer geometries in the reservoir formation that is penetrated by a salt body (Section 2). A tetrahedral mesh is available for the input model, and we define the geometry of the two alternative models of the horizon to insert by a scalar field on this input mesh. The two model meshes are then directly obtained by inserting the new horizon (Section 3). Modifying consistently a 3D meshed model while keeping its quality is a challenge that we tackle using advanced capabilities of Mmg (2022), the open-source remeshing tool our workflow relies upon. In Section 4, we use these tetrahedral meshes to run simulations of CO₂ injection with GEOSX and to compare the impact of the caprock geometries on the efficiency of the CO₂ trapping. The consistent results show the benefits of our workflow to efficiently screen different geological scenarios while keeping the mesh modifications to a minimum.

## 2. Defining alternative geological models with implicit modeling

In this paper, we demonstrate our workflow on the comparison of geological layering in the vicinity of a salt dome. Numerical simulation is key to understand the pressure front propagation, the CO₂ plume displacement, mechanical constraint, and optimize injection well parameters. In this work, we do not aim at exploring the full capacity and analysis of the coupled simulation on stochastic geological models; rather, our goal is to demonstrate a proof of concept on a relatively simple but realistic case. The first step is to build the synthetic geological models.

<table><caption>**Table 1**

Lithological and petrophysical characterization of the geological model.</caption><thead><tr><th>Color</th><th>Geological Formation</th><th>Lithology</th><th>Porosity</th><th>Permeability (m<fnref n="2" />)</th></tr></thead><tbody><tr><td></td><td>Overburden</td><td>Limestone</td><td>0.2</td><td>1e-14 (10 mD)</td></tr><tr><td></td><td>Reservoir</td><td>Sandstone</td><td>0.2</td><td>1e-11 (10 D)</td></tr><tr><td></td><td>Underburden</td><td>Limestone</td><td>0.2</td><td>1e-14 (10 mD)</td></tr><tr><td></td><td>Salt geobody</td><td>Halite</td><td>0.1</td><td>1e-18 (10 μD)</td></tr><tr><td></td><td colspan="4">Inserted layer</td></tr><tr><td></td><td>Cap rock</td><td>Shale</td><td>0.01</td><td>1e-16 (0.1 mD)</td></tr></tbody></table>

**Input model.** We start with an input synthetic model (Fig. 2) of a diapir salt body crossing three rock units that is built from the available online model<fnref n="2" /> proposed in Pellerin et al. (2015). The layer deformation reproduces typical geometries encountered in salt tectonics (Jackson and Hudec, 2017). From bottom to top, the model of 1.5 km × 1.5 km × 1 km is constituted of the halite base, a limestone formation, the sandstone reservoir formation, and a limestone top formation. For simplicity, we set constant petrophysical values per layer (Table 1).

The input model is a tetrahedral mesh on which the geological features (layers and horizons) are identified by tags on the tetrahedra and on their triangle facets. It was generated with Skua-Gocad<fnref n="3" /> which integrates robust mesh generation functionalities in the Finite Element Mesh Constructor module. Note that structural models and the corresponding surfaces can be generated with open-source software (de la Varga et al., 2019; Grose et al., 2021). To generate tetrahedral meshes, an efficient open-source alternative to the meshing functionalities of Skua-Gocad is Tetgen,<fnref n="4" /> provided that a watertight surface structural model is available (Si, 2015), as used in Pellerin et al. (2017). This initial unstructured mesh is conform to the interfaces separating the different rock formations (Fig. 2). By conform we mean that the tetrahedra faces of the 3D mesh match exactly the triangles of surfaces that were used to build the structural model.

**Caprock layer geometries alternatives.** We consider two possible caprock geometries in the top of the highly permeable sandstone reservoir. Assuming that rock type data is available along a well drilled on the side of the model, this low permeability caprock layer has an uncertain geometry in the vicinity of the diapir. We propose two alternatives (Fig. 3): Model A with a constant thickness caprock layer, corresponding to the simplest geometric scenario extrapolating borehole stratigraphy, and Model B with a unconformity of the caprock layer with the overburden and the reservoir rock. Model B reflects the possibility of the interplay between sediment deposition and erosion during diapir growth which is known to produce halokinetic sequences characterized by stratigraphic unconformities possibly surrounding the diapir (Giles and Rowan, 2012).

**Implicit definition of the caprock geometry.** To define the caprock layer geometries we rely on implicit structural modeling (also called level-set representation). An implicit 3D surface is defined as the isovalue of a scalar field. In our case, the scalar field values are stored at the nodes of a mesh and linearly interpolated in the cells. In geological modeling this representation is the base of the powerful implicit stratigraphic modeling strategy, which is the base of modern geomodeling tools. We follow the approach used in Caumon et al. (2013), Collon and Caumon (2017), Frank et al. (2007) to build the scalar fields on the input tetrahedral mesh of the model. We use the RINGToolKit Gocad

<footnote n="2">² https://www.ring-team.org/ring-dl/public/models/syntheticComplexityModels.zip</footnote>

<footnote n="3">³ https://www.pdgm.com/products/skua-gocad</footnote>

<footnote n="4">⁴ https://www.wias-berlin.de/software/index.jsp?id=TetGen</footnote>
## Section 2. Definition of the initial model and the different scenarios with implicit modeling

<figure><img src="image_2.png" /></figure>

## Section 3. Meshing the models by local modification of the input mesh

<figure><img src="image_3.png" /></figure>

## Section 4. Simulation of CO₂ injection

<figure><img src="image_4.png" /><figcaption>**Fig. 1.** Overview of the proposed workflow as presented in this work. 1a. An input geological model is meshed with tetrahedra, 1b. A layering is defined by a scalar field (Section 2) allowing to create alternative models. 2. The tetrahedral mesh is locally modified by inserting an isovalue of the scalar field using Mmg and Lumos3D (Section 3). 3. CO₂ injection simulations are run with GEOSX. 4. The results are analyzed and visualized with ParaView (Section 4).</figcaption></figure>

plugin<fnref n="5" /> that permits to compute a scalar field on the input mesh nodes from a surface horizon using a least-square discrete optimization (Frank

et al., 2007). Note that an open-source implementation of the method exists in LoopStructural (Grose et al., 2021).

## 3. Meshing the models by local modification of the input mesh

In this section, we generate tetrahedral meshes of model A and model B that are adequate for CO₂ injection simulation with GEOSX.

<footnote n="5">⁵ https://www.ring-team.org/technologies/20-web-site/softwares/gocad-plugins/195-ringtoolkit</footnote>
<figure><img src="image_5.png" /><figcaption>Fig. 2. Input synthetic model of the salt diapir. (Left) The surface geomodel representing the interfaces between rock layers triangle surfaces. (Right) The initial tetrahedral mesh of the 4 geological formations (8,617 vertices and 45,894 tetrahedra). Rock types are identified for each tetrahedron according to Table 1.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Fig. 3. The two scenarios considered for the caprock base geometries in the vicinity of the salt diapir. Model A has a constant thickness caprock layer (deep purple) and Model B has an unconformal caprock layer. The interfaces to insert are defined as an isovalue of a 3D scalar field (bottom).</figcaption></figure>

Since implementing robust and efficient mesh generation or mesh modification algorithms is extremely challenging in 3D, some might even say a life-time achievement (Loseille, 2019), we take advantage of existing open-source implementations for tetrahedral meshes: (Mmg, 2022; Dapogny et al., 2014), Gmsh (Geuzaine and Remacle, 2009), Tetgen (Si, 2015), and focus on their consistent integration in geomodeling workflows. We rely in particular on the Mmg platform, that provides functionalities to insert level-set interfaces in tetrahedral meshes. This allows us to extend to 3D the approach of Legentil et al. (2022) who proposed to insert a line modeling a fluid contact in a 2D model (a line in a triangle mesh) to evaluate the impact of its position on wave propagation. In 3D, we modify locally the tetrahedral mesh to insert a surface implicitly defined. The main advantage of this incremental approach is that we avoid two time-consuming steps: (i) the update of the surface structural model and (ii) the generation of a new 3D mesh from scratch.

## 3.1. Constraints for meshing

Before diving in mesh operations, let us review the requirements that our global workflow (Fig. 1) imposes on the meshing step. Our objective is to simulate CO₂ injection and migration with GEOSX. To run the simulations efficiently and to obtain consistent results, the generated mesh should be valid. Valid means that all cells have a strictly positive volume and do not contain any degenerate edge or

GEOSX operates on unstructured conformal meshes. Conformal means that the intersection of two cells (tetrahedra) is either empty or is a face (triangle), edge, or point shared by the two cells. Generating valid conformal meshes is a necessary condition, but is not a sufficient one to obtain consistent results. A valid mesh with a bad quality can prevent simulation convergence or produce non-physical results. (Knupp, 2007) defines mesh quality as “*the characteristics of a mesh that permit a particular numerical PDE simulation to be efficiently performed, with fidelity to the underlying physics, and with the accuracy required for the problem*”. He further details the related inherent difficulties and challenges of defining reliable mesh quality metrics. As GEOSX does not yet provide mesh quality metrics for tetrahedral meshes, we limit our objective to generate meshes that do not contain defective elements with very small or large angles, and that have a consistent size.

Our objective of comparing simulations on similar structural models adds another constraint. The insertion of the new horizon modifies the input mesh, however these modifications should be kept local. This is a crucial point for geological modeling at basin scale as typically considered for CO₂ storage.

## 3.2. Implicit surface insertion

In this section, we insert the level-set interfaces defined in Section 2 in the input tetrahedral mesh (Fig. 2) to generate meshes of Model A and Model B. To restrict the set of tetrahedra that are remeshed, we first
<figure><img src="image_7.png" /><figcaption>Fig. 4. The intersection of good quality tetrahedra with a level set typically generates several very small as well as badly shaped tetrahedra. The level set (blue) is inserted near the apex of two adjacent tetrahedra. The red tetrahedron is subdivided into 4 tetrahedra, the green tetrahedron into 6 tetrahedra. These bad quality tetrahedra will be replaced by better shaped elements during the mesh optimization of Mmg.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Fig. 5. Processing of scars at the intersection of the new interface with the input model interfaces. (a) The red tetrahedra are incident to edges created at those intersections (tag Mmg: 10). (b) Close view on the scar left by the interface insertion in the mesh since the modification of tetrahedra outside the reservoir (inside the salt dome) is prevented. (c) The remeshing of the buffer zones that encompass the intersections permit to recover good quality tetrahedra.</figcaption></figure>

define a modification zone, then the level-set iso-value is discretized
before being post-processed to optimize mesh quality.

**Restricting the modification zone.** The modification region is a set of tetrahedra of the input mesh where the level-set interface will be inserted. We build this set of tetrahedra as the intersection of the set of tetrahedra that are contained in the reservoir layer and the set of tetrahedra that are at a distance to the level-set inferior to a given threshold. Additional constraints could be easily added provided that the region where the level-set is to be inserted is large enough to perform mesh modification operations and to optimize mesh quality.

**Interface discretization with mmg.** We insert the surface in the mesh with Mmg. The first step is to cut each tetrahedron intersected by the level-set into pieces ensuring the conformity between the newly created tetrahedra. Note that this intersection is much more complex in 3D than in 2D: in 2D, a triangle is either cut in two triangles if the level-set intersects one vertex and one edge, or in three triangles if the level-set intersects two edges. In 3D, the number of configurations is larger, and the consistency of the subdivision between adjacent tetrahedra should be guaranteed. An example is given in Fig. 4: two adjacent tetrahedra are cut respectively into 4 and 6 tetrahedra which are typically defective elements in the view of the mesh quality discussed in the previous section. The closer the level-set is to the nodes and edges of a tetrahedron, the lower the quality of the elements resulting from

the intersection. The second step is, therefore, to optimize the quality of
tetrahedra of the modification region with Mmg, which performs local
edge flips, collapses or splits (Dapogny et al., 2014).

**Scar removal and quality improvement.** Mmg remeshing functionality drastically improves the quality of the worst tetrahedra created during the surface insertion. However, a few ill-shaped tetrahedra may remain. Indeed, some tetrahedra cut into pieces when inserting the interface are outside the modified region. They were split to keep the mesh conformity and should be post-processed. These tetrahedra have a face on one of the input interface and are intersected by the level-set (red tetrahedra on Fig. 5-a). They are subdivided into lower quality tetrahedra that may be problematic for running the simulations (Fig. 5-b).

Therefore, we propose a post-processing step to optimize the quality
of these tetrahedra after the level set discretization. We create a buffer
zone around them to remove the remaining scars of the intersection
to remesh. The buffer zone initialized to the red tetrahedra of Fig. 5,
is then enlarged using a greedy iterative algorithm based on node
adjacencies. At each iteration, tetrahedra containing at least one node
of a tetrahedron in the buffer at the previous iteration are added to the
buffer. The buffer zone should be large enough to accommodate mesh
improvements, but should remain small enough to limit the impact
on the complete mesh. For Model A and Model B we used only two
iterations. The resulting meshes are obtained in 10 s (machine with
<figure><img src="image_9.png" /><figcaption>Fig. 6. Tetrahedral meshes obtained for Model A (48,935 tetrahedra) and Model B (48,888 tetrahedra) with Mmg parameters used: hmin = hmax =75 m and hausd = 10 m.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Fig. 7. Comparison of general purpose a priori mesh quality metrics computed on Model A and Model B. The obtained values for aspect ratio, dihedral angles, and edge lengths are similar, the green color indicating overlapping values. Note that the difference of maximal values of the aspect ratios is due to the lower angle of the contact between the cap rock and the reservoir layer in the models.</figcaption></figure>

processor Intel(R) Xeon(R) CPU E5-1620 v2 @3.70 GHz, 8 Cores)
(Fig. 6).

Resulting meshes for simulation. The quality of the obtained mesh is checked to prevent well-known defective elements for simulations. We compute general purpose metrics using the Verdict library in Par-aView (Stimpson et al., 2007) (Fig. 7). We consider the distributions of three geometrical measures on tetrahedra: the aspect ratio, measured as $\frac{L_{max}}{2\sqrt{6r}}$ (with $L_{max}$ is the maximal edge length and $r$ the radius of the inscribed sphere), the minimal dihedral angle and the edge lengths. The small angle introduced by the caprock unconformity of Model B has a clear impact on the minimal value of the minimal angle as well as the maximal value of the aspect ratio.

3.3. *Implementation*

In this section we provide details on our implementation to ensure
the reproducibility of our results.

Mmg platform. (Mmg, 2022) is an open source software for triangle and tetrahedral remeshing. It provides tools (1) to generate, adapt and optimize 2D triangulations, (2) to adapt and optimize 3D surface triangulation and isovalue discretization, and (3) to adapt and optimize tetrahedral meshes and implicit domain meshes. Mmg has been mostly developed for computational fluid dynamics applications and for shape optimization applications (Dobrzynski, 2005; Dapogny et al., 2014). Available as a C library or as an executable, Mmg can be easily integrated in a workflow and its functionalities were demonstrated on simulation of the subsonic turbulent nozzle jet flow, to optimize an elastic mechanical system, to accelerate the full waveform inversion, see e.g. Balarac et al. (2021), Jacquet (2021). Mmg is an open source software jointly developed by Inria,<fnref n="6" /> the University of

Bordeaux, Bordeaux INP,<fnref n="7" /> the UPMC<fnref n="8" /> and the CNRS.<fnref n="9" /> Mmg is multi-platform (Windows, Linux and MacOS), distributed under LGPL license and available on Github.<fnref n="10" />

Lumos3D. To generate automatically the meshes for a potentially large number of alternative geological models, we implemented an extension of the library Lumos2D (Legentil et al., 2022): the Lumos3D library,<fnref n="11" /> which is written in C++ and depends on the C mmg application programming interface (API). From an input meshed model, it allows to restrict the discretization and the optimization to the modified region. To prevent remeshing of the entire mesh, the options `-nomove`, `-noinsert`, `-noswap` are used by calling the Mmg API for level-set discretization. The `-nomove` option blocks the movement of vertices and the `-noinsert` one prevents the creation of new vertices. The `-noswap` option prevents the edge swapping between two adjacent tetrahedra. For the optimization step, the options `hmin` and `hmax` are also used to set the mesh size to adapt the mesh resolution. To control the maximal distance between the piecewise linear representation boundary representation and the reconstructed boundary, the `-hausd` option is set. A decreasing Hausdorff distance, decreases the distance between the level-set and the triangle surface, high curvature areas are then refined.

Mmg uses the Medit file format (`.mesh`). We encode geological model entities with integer tags defined on the tetrahedra. The scalar field and/or metric for remeshing values are defined on the vertices of the mesh and are stored in another file format (`.sol`). The local parameters for multi-material models are defined in a Mmg file (`.mmg3d`).

4. Numerical simulation of CO₂ injection

Having built tetrahedral meshes for Model A and Model B, we
run simulations of CO₂ injection in the vicinity of the salt dome with

<footnote n="6">⁶ https://www.inria.fr/en</footnote>

<footnote n="7">⁷ https://www.bordeaux-inp.fr/en</footnote>

<footnote n="8">⁸ https://www.sorbonne-universite.fr/</footnote>

<footnote n="9">⁹ https://www.cnrs.fr/en</footnote>

<footnote n="10">¹⁰ https://github.com/MmgTools/mmg</footnote>

<footnote n="11">¹¹ https://www.ring-team.org/technologies/420-lumos</footnote>
GEOSX. The results are then compared to assess the role of stratigraphic unconformities on trap integrity.

## 4.1. GEOSX

GEOSX is an open-source, multiphysics, and multilevel physics simulation tool jointly created by Lawrence Livermore National Laboratory, Stanford University, and TotalEnergies (Gross and Mazuyer, 2021). GEOSX<fnref n="12" /> is a code actively developed, it is distributed under LGPL 2.1 license terms. Designed for scalability on multiple CPUs and multiple GPUs, GEOSX offers a suite of physical solvers that rely on various numerical schemes such as Finite Volume methods and Finite Element methods. These different schemes are implemented for unstructured grids to accommodate the complex geometries of the sub-surface. For flow simulations, GEOSX offers the choice to use the classical two-point flux approximation (TPFA) or a more advanced hybrid mimetic discretization scheme (MFD), see e.g. Lipnikov et al. (2014). MFD scheme is more accurate on unstructured grids but more computationally intensive. GEOSX has the flexibility of assigning solvers to specific regions (flow, or coupled flow and poromechanics, or mechanics) of the field or specific times of the simulation. This is an advantage to be efficient when handling large geological models, such as basins and saline aquifers, over long simulation periods (Gross and Mazuyer, 2021).

## 4.2. Flow simulation in porous media

The injection of CO2 in water is described as a multiphase multi-component system. We summarize the equations governing this physical process following the documentation of GEOSX.<fnref n="13" />

*Mass conservation equations.* For a component c, the mass conservation is expressed as:

$$ \phi \frac{\partial}{\partial t} \left( \sum_{\ell} \rho_{\ell} y_{c\ell} S_{\ell} \right) + \nabla \cdot \left( \sum_{\ell} \rho_{\ell} y_{c\ell} \mathbf{u}_{\ell} \right) - \sum_{\ell} \rho_{\ell} y_{c\ell} q_{\ell} = 0, \quad (1) $$

where $\phi$ is the porosity of the medium, $S_{\ell}$ the saturation of phase $\ell$, $y_{c\ell}$ the mass fraction of component c in phase $\ell$, $\rho_{\ell}$ the phase density, $q_{\ell}$ the volumetric source term, and t the time. The formulation used is isothermal.

*Darcy's law.* The multiphase extension of Darcy's law is used, with the phase velocity $\mathbf{u}_{\ell}$ written as a function of the phase potential gradient $\nabla\Phi_{\ell}$:

$$ \mathbf{u}_{\ell} := -k\lambda_{\ell}\nabla\Phi_{\ell} = -k\lambda_{\ell}(\nabla(p - P_{c,\ell}) - \rho_{\ell}g\nabla z), \quad (2) $$

with k the rock permeability, $\lambda_{\ell} = k_{\ell}/\mu_{\ell}$ the phase mobility, defined as the phase relative permeability divided by the phase viscosity, p the reference pressure, $P_{c,\ell}$ the capillary pressure, g the gravitational acceleration, and z the depth. The computation of the component properties is described in the following paragraph.

The combination of the mass conservation equations with Darcy's law gives a set of $n_c$ equations written as:

$$ \phi \frac{\partial}{\partial t} \left( \sum_{\ell} \rho_{\ell} y_{c\ell} S_{\ell} \right) - \nabla \cdot \mathbf{k} \left( \sum_{\ell} \rho_{\ell} y_{c\ell} \lambda_{\ell} \nabla \Phi_{\ell} \right) - \sum_{\ell} \rho_{\ell} y_{c\ell} q_{\ell} = 0. \quad (3) $$

*Constraints and thermodynamic equilibrium.* The pore space is completely filled by the phases, which can be put as a volume constraint such as: $\sum_{\ell} S_{\ell} = 1$. The system is also closed by the thermodynamic equilibrium constraint: $f_{c\ell} - f_{cm} = 0$, with $f_{c\ell}$ the fugacity of component c in phase $\ell$ and $f_{cm}$ the fugacity of component c in phase

m. The thermodynamic equilibrium is ensured through a constitutive fluid model described in the next paragraph. To summarize, the compositional multiphase flow solver assembles a set of $n_c + 1$ equations in each element, i.e., $n_c$ mass conservation equations and one volume constraint equation. The thermodynamic equilibrium is ensured at each nonlinear iteration by the fluid constitutive model.

*CO₂ brine model.* The fluid model includes two components that are transported by one or two fluid phases. The update of the fluid properties is done in two steps:

* The phase fractions ($v_p$) and phase component fractions ($y_{c,p}$) are computed as a function of pressure (P), temperature (T), component fractions ($z_c$), and a constant salinity.

* The phase densities ($\rho_p$) and phase viscosities ($\mu_p$) are computed as a function of pressure, temperature, the updated phase component fractions, and a constant salinity.

The CO₂ phase density and viscosity are computed using the relation defined by Span and Wagner (1996), and the brine density and viscosity by Phillips correlation (Phillips et al., 1981).

## 4.3. Simulation parameterization

To simulate the injection of CO₂ in both models, we run multiphase flow simulations using a Two-Point Flux-Approximation (TPFA) scheme to compute the transmissibilities between cells. In addition to the meshes where geological features are identified and the rock properties defined, we build the two parameter files for GEOSX (these xml files are provided in the supplementary material). We set the following parameters for the simulation:

* The mesh is provided using .msh format. The injection well is reduced to one tetrahedron defined by its center that is located at depth $z = -516$ m. To ensure a fixed position for the well between model A and model B, the injection cell is chosen outside the modified region (the red tetrahedron in Fig. 8) and specified following the documentation.<fnref n="14" />

* The multiphase flow solver is chosen to be TPFA following the model described in Section 4.2. The well solver is described in the same xml entry, using the well constraint with an injection rate defined in an additional section, we inject $0.03 \text{ m}^3/\text{s}$ during $1.6 \times 10^8 \text{ s}$ (about 5 years). The simulation is run for $5 \times 10^8 \text{ s}$ (about 16 years). Injection and post-injection periods as well as the output periodicity are defined as events in GEOSX terms.

* The different regions of the models (Section 2) are defined in the Element Regions section from the tags on tetrahedra. On each region, a constitutive model is defined, i.e., rock properties and fluid models.

* The model is initialized using the field specifications entry. This initialization is done at a constant pressure and constant composition (reservoir filled with water). The equilibrium is done at the first calculation step.

* The user requests the desired outputs, here we used 3D outputs in vtk file format.

## 4.4. Results

The results of the CO₂ injection with GEOSX are visualized with the open-source scientific visualization tool ParaView<fnref n="15" /> (Fig. 8). Only cells where CO₂ concentration is non-zero are visible on top of the background mesh presented in Fig. 6. For both models, the CO₂ plume shape

<footnote n="12">¹² https://github.com/GEOSX/GEOSX</footnote>

<footnote n="13">¹³ https://geosx-geosx.readthedocs-hosted.com/</footnote>

<footnote n="14">¹⁴ https://geosx-geosx.readthedocs-hosted.com/en/latest/docs/sphinx/basicExamples/multiphaseFlowWithWells/Example.html</footnote>

<footnote n="15">¹⁵ https://www.paraview.org/</footnote>
<figure><img src="image_11.png" /><figcaption>**Fig. 8.** Simulation results on Model A (left) and Model B (right) after 3 years of injection. The stratigraphic trap is efficient for Model A, whereas a gas leakage is visible when the caprock layer is affected by erosion for Model B. The tetrahedral cells where CO₂ is injected are highlighted in red.</figcaption></figure>

<table><caption>**Table 2**
CO₂ mass for each region in the different models after three years of injection.</caption><thead><tr><th>mass in kg</th><th></th><th>Model A</th><th>Model B</th></tr></thead><tbody><tr><td rowspan="3">Overburden</td><td>Free</td><td>1.51e<sup>8</sup></td><td>8.55e<sup>8</sup></td></tr><tr><td>Dissolved</td><td>1.42e<sup>8</sup></td><td>2.83e<sup>8</sup></td></tr><tr><td><i>All</i></td><td><b>2.92e<sup>8</sup></b></td><td><b>1.14e<sup>9</sup></b></td></tr><tr><td rowspan="3">Reservoir</td><td>Free</td><td>1.52e<sup>9</sup></td><td>7.46e<sup>8</sup></td></tr><tr><td>Dissolved</td><td>2.44e<sup>8</sup></td><td>2.34e<sup>8</sup></td></tr><tr><td><i>All</i></td><td><b>1.76e<sup>9</sup></b></td><td><b>9.80e<sup>8</sup></b></td></tr><tr><td rowspan="3">Caprock</td><td>Free</td><td>5.36e<sup>7</sup></td><td>1.06e<sup>6</sup></td></tr><tr><td>Dissolved</td><td>5.93e<sup>6</sup></td><td>4.48e<sup>5</sup></td></tr><tr><td><i>All</i></td><td>5.96e<sup>7</sup></td><td>1.51e<sup>6</sup></td></tr><tr><td rowspan="3">Underburden</td><td>Free</td><td>6.35e<sup>7</sup></td><td>6.23e<sup>7</sup></td></tr><tr><td>Dissolved</td><td>1.70e<sup>7</sup></td><td>1.68e<sup>7</sup></td></tr><tr><td><i>All</i></td><td>8.06e<sup>7</sup></td><td>7.91e<sup>7</sup></td></tr><tr><td rowspan="3">Salt</td><td>Free</td><td>9.58e<sup>3</sup></td><td>0.00e<sup>0</sup></td></tr><tr><td>Dissolved</td><td>9.59e<sup>6</sup></td><td>5.78e<sup>6</sup></td></tr><tr><td><i>All</i></td><td>9.60e<sup>6</sup></td><td>5.78e<sup>6</sup></td></tr><tr><td>Full Model</td><td></td><td><b>2.2055e<sup>9</sup></b></td><td><b>2.2042e<sup>9</sup></b></td></tr></tbody></table>

is consistent: CO₂ flows upwards and CO₂ concentration increases just below the caprock. For Model A, we observe that the caprock prevents the CO₂ migration towards the overburden. For Model B, the CO₂ concentration in the overburden is higher than for Model A, the plume goes through the caprock. As expected, the structural trap is not perfect, there is a leakage of CO₂. These results are confirmed by the mass of CO₂ in each region of the model after 3 years of injection (Table 2).

**Reproducibility.** To ensure the reproducibility of the results presented in this paper, we provide all input, intermediate results and output files. They are available in the supplementary materials of the paper, as described in the list below:

* Input mesh of the salt diapir (inputmesh.vtk)

* Level sets for Model A and Model B (levelsetA.sol and levelsetB.sol)

* Tetrahedral meshes of Model A and Model B (meshA.msh and meshB.msh)

* GEOSX input parameter files (simulationA.xml and simulationB.xml)

* Fluid models (fluid.txt)

* Outputs of the simulation after 3 years (outputA.vtk and outputB.vtk)

We provide details about the software used in Table 3.

## 5. Discussion

In this paper, we set the ground for and incremental geometric model updating workflow and showed how it can be used for a relatively simple CO₂ sequestration case study. Nonetheless, several aspects can be discussed and extended in future studies.

For demonstration purposes, we chose constant permeability and porosity per layer to demonstrate the workflow up to simulation. When working with real data models, property geostatistical filling is a key step to represent geological features. On unstructured grids, classical algorithms operating on structured grids cannot be used, but recent geostatistical algorithms, e.g. Biver et al. (2019), Mourlanette et al. (2020) have demonstrated their capabilities on tetrahedral meshes to fill heterogeneous porosity, permeability, and facies values consistent with available data.

Recently (Osman et al., 2021) have shown some examples where petrophysical heterogeneity can be neglected when key geological interfaces are correctly represented in the subsurface model. The proposed surface insertion approach could be used to adaptively test which heterogeneities are best represented as surfaces and to what extent it complements existing volumetric representations of geological heterogeneities in a systematic way.

Therefore, a remaining challenge is to account for a larger range of geological settings. The management of faults is key for this objective. Faults are discontinuities, and fault slips mean discontinuities in the level set inserted to model a horizon. From a geometrical standpoint, there is another important difference: faults may end in the model while horizons never do. To account for faults cutting the model in two blocks, inserting iteratively the horizons in the fault blocks is a possibility with the current workflow. The fault blocks would then be considered independently while maintaining conformity. The additional challenge is to account for faults that are ending in the model.

Remeshing also raises practical challenges for automation. For instance, we chose to inject the CO₂ right below the reservoir in order to keep the same parameter set. Some adaptation would be needed to specify boundary conspirions in the simulator, even if the mesh changes. More significantly, exploring effectively a larger range of scenarios would call for a tighter integration of the model construction, perturbation and mesh updating. As mentioned before, the combination of LUMOS and Mmg with the open-source implicit modeling code LoopStructural could be a valuable option to explore this in the future. Indeed, alternative level set geometries could be generated with LoopStructural before updating the mesh with LUMOS and Mmg.
<table><caption>**Table 3**
Software list (OS stands for Open-source).</caption><thead><tr><th>Software</th><th>Purpose</th><th>Version</th><th>URL</th><th>OS</th><th>Availability</th></tr></thead><tbody><tr><td>Skua-Gocad</td><td>Geological modeling</td><td>2021</td><td>https://pdgm.com/products/skua-gocad/</td><td>No</td><td>Commercial</td></tr><tr><td>RINGToolKit</td><td>Level set building</td><td>2021</td><td>https://www.ring-team.org/technologies/</td><td>No</td><td>Ring Consortium</td></tr><tr><td>Lumos3D</td><td>Workflow</td><td>2022</td><td>https://www.ring-team.org/technologies/</td><td>No</td><td>Ring Consortium</td></tr><tr><td>Mmg</td><td>Remeshing platform</td><td>5.6</td><td>https://www.mmgtools.org/</td><td>Yes</td><td>LGPL3</td></tr><tr><td>GEOSX</td><td>Multiphysics solver</td><td>-</td><td>https://github.com/GEOSX/GEOSX</td><td>Yes</td><td>LGPL2.1</td></tr><tr><td>ParaView</td><td>Scientific visualization</td><td>5.10</td><td>https://www.paraview.org/</td><td>Yes</td><td>BSD</td></tr></tbody></table>

Finally, one major challenge is the mesh quality, a topic that is at the intersection of geometric modeling and numerical simulation. To our best knowledge, no mesh generation method provides guarantees on its output mesh for complex geometries such as those of geological models. However, one benefit of the level-set-based mesh updating is that the surface, which representation is stored directly on the mesh, can be modified more easily, simplifying mesh quality optimization strategies by snap-rounding (see Legentil et al., 2022 for the 2D example and further discussion).

## 6. Conclusion

We propose a workflow for locally updating 3D geological meshed models and demonstrate it to compare the impact of two layering scenarios on the simulation of CO₂ injection. This automatic workflow extends to 3D the method proposed in Legentil et al. (2022). We use open-source software to generate the modified meshed models, to run the CO₂ injection simulation as well as to visualize and to compare the results. Although we made several simplifying assumptions, it opens a very promising path for developing numerical simulations of CO₂ injection on unstructured meshes and incremental approaches for geological model building and uncertainty assessment.

### CRediT authorship contribution statement

**Capucine Legentil:** Conceptualization, Investigation, Methodology, Software, Writing – original draft. **Jeanne Pellerin:** Conceptualization, Methodology, Supervision, Funding acquisition, Validation, Writing – original draft. **Margaux Raguenel:** Conceptualization, Investigation, Methodology, Validation, Writing – review & editing. **Guillaume Caumon:** Conceptualization, Methodology, Supervision, Funding acquisition, Validation, Writing – review & editing.

### Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

### Data availability

Data are available as supplementary materials, and software list with links is provided in the article.

### Acknowledgments

This work was performed in the frame of the RING project (http://ring.georessources.univ-lorraine.fr) supported by the RING Consortium managed by ASGA. This work is additionally supported by TotalEnergies. We acknowledge Inria and the Mmg Consortium for Mmg libraries. The authors would also like to thank Gilles Darche and Francois Hamon for their help, and Florian Wellmann, an anonymous reviewer and the editor Boyan Brodaric for their constructive comments.

## Appendix A. Supplementary data

Supplementary material related to this article can be found online at https://doi.org/10.1016/j.acags.2023.100118.

## References

Ajayi, T., Gomes, J.S., Bera, A., 2019. A review of CO₂ storage in geological formations emphasizing modeling, monitoring and capacity estimation approaches. Pet. Sci. 16 (5), 1028–1063. http://dx.doi.org/10.1007/s12182-019-0340-8, URL: http://link.springer.com/10.1007/s12182-019-0340-8.

Balarac, G., Basile, F., Bénard, P., Bordeu, F., Chapelier, J.-B., Cirrottola, L., Caumon, G., Dapogny, C., Frey, P., Froehly, A., et al., 2021. Tetrahedral remeshing in the context of large-scale numerical simulation and high performance computing.

Benson, S.M., Cole, D.R., 2008. CO₂ sequestration in deep sedimentary formations. Elements 4 (5), 325–331. http://dx.doi.org/10.2113/gselements.4.5.325, URL: https://pubs.geoscienceworld.org/elements/article/4/5/325-331/137789.

Biver, P., Fuet, S., Allard, D., 2019. Direct geostatistical simulation on unstructured grids I: Recent improvements for additive variables. In: Petroleum Geostatistics 2019. European Association of Geoscientists & Engineers, Florence, Italy, pp. 1–5. http://dx.doi.org/10.3997/2214-4609.201902231, URL: https://www.earthdoc.org/content/papers/10.3997/2214-4609.201902231.

Bond, C.E., Gibbs, A.D., Shipton, Z.K., Jones, S., et al., 2007. What do you think this is? "Conceptual uncertainty" in geoscience interpretation. GSA Today 17 (11), 4, Publisher: THE GEOLOGICAL SOCIETY OF AMERICA, INC.

Broseta, D., Tonnet, N., Shah, V., 2012. Are rocks still water-wet in the presence of dense CO₂ or H₂s? Geofluids 12 (4), 280–294. http://dx.doi.org/10.1111/j.1468-8123.2012.00369.x, arXiv:https://onlinelibrary.wiley.com/doi/pdf/10.1111/j.1468-8123.2012.00369.x, URL: https://onlinelibrary.wiley.com/doi/abs/10.1111/j.1468-8123.2012.00369.x.

Caumon, G., Gray, G., Antoine, C., Titeux, M.O., 2013. Three-dimensional implicit stratigraphic model building from remote sensing data on tetrahedral meshes: Theory and application to a regional model of la Popa Basin, NE Mexico. IEEE Trans. Geosci. Remote Sens. 51 (3), 1613–1621. http://dx.doi.org/10.1109/TGRS.2012.2207727, Publisher: IEEE.

Caumon, G., Sword, C.H., Mallet, J.-L., 2003. Constrained modifications of non-manifold B-reps. In: Proceedings of the Eighth ACM Symposium on Solid Modeling and Applications. SM '03, ACM Press, New York, New York, USA, p. 310. http://dx.doi.org/10.1145/781650.781657.

Celia, M.A., Bachu, S., Nordbotten, J.M., Bandilla, K.W., 2015. Status of CO₂ storage in deep saline aquifers with emphasis on modeling approaches and practical simulations: Status of CO₂ storage in deep saline aquifers. Water Resour. Res. 51 (9), 6846–6892. http://dx.doi.org/10.1002/2015WR017609, URL: http://doi.wiley.com/10.1002/2015WR017609.

Clausolles, N., Collon, P., Caumon, G., 2019. Generating variable shapes of salt geobodies from seismic images and prior geological knowledge. Interpretation 7 (4), T829–T841, Publisher: Society of Exploration Geophysicists and American Association of Petroleum.

Collon, P., Caumon, G., 2017. 3D geomodelling in structurally complex areas - implicit vs. explicit representations. In: 79th EAGE Conference and Exhibition 2017. http://dx.doi.org/10.3997/2214-4609.201701144.

Dapogny, C., Dobrzynski, C., Frey, P., 2014. Three-dimensional adaptive domain remeshing, implicit domain meshing, and applications to free and moving boundary problems. J. Comput. Phys. 262, 358–378. http://dx.doi.org/10.1016/j.jcp.2014.01.005.

de la Varga, M., Schaaf, A., Wellmann, F., 2019. Gempy 1.0: Open-source stochastic geological modeling and inversion. Geosci. Model Dev. 12 (1), 1–32. http://dx.doi.org/10.5194/gmd-12-1-2019, URL: https://gmd.copernicus.org/articles/12/1/2019/.

Dobrzynski, C., 2005. Adaptation de Maillage anisotrope 3D et application à l'aéro-thermique des bâtiments (Ph.D. thesis). Université Pierre et Marie Curie - Paris VI, URL: https://tel.archives-ouvertes.fr/tel-00120327.

Fleury, M., Pironon, J., Le Nindre, Y., Bildstein, O., Berne, P., Lagneau, V., Broseta, D., Pichery, T., Fillacier, S., Lescanne, M., Vidal, O., 2011. Evaluating sealing efficiency of caprocks for CO₂ storage: An overview of the geocarbon integrity program and results. Energy Procedia 4, 5227–5234. http://dx.doi.org/10.1016/j.egypro.2011.02.501, URL: https://www.sciencedirect.com/science/article/pii/S1876610211007806, 10th International Conference on Greenhouse Gas Control Technologies.

Frank, T., Tertois, A.L., Mallet, J.L., 2007. 3D-reconstruction of complex geological interfaces from irregularly distributed and noisy point data. Comput. Geosci. 33 (7), 932–943. http://dx.doi.org/10.1016/j.cageo.2006.11.014.
Geuzaine, C., Remacle, J.F., 2009. Gmsh: A 3-D finite element mesh generator with built-in pre- and post-processing facilities. Internat. J. Numer. Methods Engrg. 79 (11), 1309–1331. http://dx.doi.org/10.1002/nme.2579.

Giles, K.A., Lawton, T.F., 2002. Halokinetic sequence stratigraphy adjacent to the El Papalote diapir, Northeastern Mexico. AAPG Bull. 86 (5), 823–840, Publisher: American Association of Petroleum Geologists (AAPG).

Giles, K.A., Rowan, M.G., 2012. Concepts in halokinetic-sequence deformation and stratigraphy. Geol. Soc. Lond. Special Publ. 363 (1), 7–31, Publisher: Geological Society of London.

Grose, L., Ailleres, L., Laurent, G., Jessell, M., 2021. LoopStructural 1.0: time-aware geological modelling. Geosci. Model Dev. 14 (6), 3915–3937. http://dx.doi.org/10.5194/gmd-14-3915-2021, URL: https://gmd.copernicus.org/articles/14/3915/2021/.

Gross, H., Mazuyer, A., 2021. GEOSX: A multiphysics, multilevel simulator designed for exascale computing. In: SPE Reservoir Simulation Conference, SPE, On-Demand, http://dx.doi.org/10.2118/203932-MS, D011S010R007, URL: https://onepetro.org/spersc/proceedings/21RSC/1-21RSC/D011S010R007/470761.

Hamza, A., Hussein, I.A., Al-Marri, M.J., Mahmoud, M., Shawabkeh, R., Aparicio, S., 2021. CO2 enhanced gas recovery and sequestration in depleted gas reservoirs: A review. J. Pet. Sci. Eng. 196, 107685. http://dx.doi.org/10.1016/j.petrol.2020.107685, URL: https://linkinghub.elsevier.com/retrieve/pii/S0920410520307518.

Issautier, B.t., Viseur, S., Audigane, P., Chiaberge, C., Le Nindre, Y.-M., 2015. A new approach for evaluating the impact of fluvial type heterogeneity in CO 2 storage reservoir modeling. C. R. Geosci.

Jackson, M.P., Hudec, M.R., 2017. Salt Tectonics: Principles and Practice. Cambridge University Press.

Jacquet, P., 2021. Inversion par forme d'ondes complète en domaine temporel utilisant des méthodes de Galerkiné Discontinues avancées (Ph.D. thesis). Université de Pau et des pays de l'Adour.

Knupp, P.M., 2007. Remarks on mesh quality. In: 45th AIAA Aerospace Sciences Meeting and Exhibit. Reno, NV, URL: https://www.osti.gov/servlets/purl/1146104.

Lallier, F., Caumon, G., Borgomano, J., Viseur, S., Fournier, F., Antoine, C., Gentilhomme, T., 2012. Relevance of the stochastic stratigraphic well correlation approach for the study of complex carbonate settings: Application to the Malampaya buildup (offshore Palawan, Philippines). Geol. Soc. Lond. Special Publ. 370 (1), 265–275. http://dx.doi.org/10.1144/SP370.12, arXiv:https://www.lyellcollection.org/doi/pdf/10.1144/SP370.12, URL: https://www.lyellcollection.org/doi/abs/10.1144/SP370.12.

Legentil, C., Pellerin, J., Cupillard, P., Froehly, A., Caumon, G., 2022. Testing scenarios on geological models: Local interface insertion in a 2D mesh and its impact on seismic wave simulation. Comput. Geosci. 159, 105013, Publisher: Elsevier.

Li, Z., Dong, M., Li, S., Huang, S., 2006. CO2 sequestration in depleted oil and gas reservoirs—caprock characterization and storage capacity. Energy Convers. Manage. 47 (11–12), 1372–1382. http://dx.doi.org/10.1016/j.enconman.2005.08.023, URL: https://linkinghub.elsevier.com/retrieve/pii/S0196890405002098.

Lipnikov, K., Manzini, G., Shashkov, M., 2014. Mimetic finite difference method. J. Comput. Phys. 257, 1163–1227. http://dx.doi.org/10.1016/j.jcp.2013.07.031, URL: https://linkinghub.elsevier.com/retrieve/pii/S0021999113005135.

Loseille, A., 2019. Some challenges in industrialising mesh adaptation : A focus on accuracy, stability, geometry and parallelism. In: Sixth Workshop on Grid Generation for Numerical Computations. INRIA Saclay Ile-de-France, Palaiseau, France, URL: https://pages.saclay.inria.fr/frederic.alauzet/TetrahedronVI/Slides_Tetrahedron_VI/Friday/Adrien_Loseille-Tetrahedron-2019.pdf.

Miall, A.D., 2016. The valuation of unconformities. Earth-Sci. Rev. 163, 22–71. http://dx.doi.org/10.1016/j.earscirev.2016.09.011, URL: https://www.sciencedirect.com/science/article/pii/S0012825216303257.

Mmg, 2022. Mmg version 5.6.0. URL: https://github.com/MmgTools/mmg.

Mourlanette, P., Biver, P., Renard, P., Nøtinger, B., Caumon, G., Perrier, Y.A., 2020. Direct simulation of non-additive properties on unstructured grids. Adv. Water Resour. 143, 103665. http://dx.doi.org/10.1016/j.adwatres.2020.103665, URL: https://linkinghub.elsevier.com/retrieve/pii/S030917082030213X.

Oldenburg, C.M., Pruess, K., Benson, S.M., 2001. Process modeling of CO₂ injection into natural gas reservoirs for carbon sequestration and enhanced gas recovery. Energy Fuels 15 (2), 293–298. http://dx.doi.org/10.1021/ef000247h, URL: https://pubs.acs.org/doi/10.1021/ef000247h.

Osman, H., Graham, G.H., Moncorge, A., Jacquemyn, C., Jackson, M.D., 2021. Is cell-to-cell scale variability necessary in reservoir models? Math. Geosci. 53 (4), 571–596. http://dx.doi.org/10.1007/s11004-020-09877-y.

Pellerin, J., Botella, A., Bonneau, F., Mazuyer, A., Chauvin, B., Lévy, B., Caumon, G., 2017. RINGMesh: A programming library for developing mesh-based geomet modeling applications. Comput. Geosci. 104, 93–100. http://dx.doi.org/10.1016/j.cageo.2017.03.005.

Pellerin, J., Caumon, G., Julio, C., Mejia-Herrera, P., Botella, A., 2015. Elements for measuring the complexity of 3D structural models: Connectivity and geometry. Comput. Geosci. 76, 130–140, Publisher: Elsevier.

Phillips, S., Igbene, A., Fair, J., Ozbek, H., Tavana, M., 1981. Technical Databook for Geothermal Energy Utilization. Technical Report LBL-12810, 6301274, http://dx.doi.org/10.2172/6301274, URL: http://www.osti.gov/servlets/purl/6301274-us7CIF/native/.

Rojo, L.A., Escalona, A., Schulte, L., 2016. The use of seismic attributes to enhance imaging of salt structures in the barents sea. First Break 34 (11), http://dx.doi.org/10.3997/1365-2397.2016014, URL: https://www.earthdoc.org/content/journals/10.3997/1365-2397.2016014.

Sadler, P.M., Jerolmack, D.J., 2015. Scaling laws for aggradation, denudation and progradation rates: the case for time-scale invariance at sediment sources and sinks. Geol. Soc. Lond. Special Publ. 404 (1), 69–88. http://dx.doi.org/10.1144/SP404.7, arXiv:https://www.lyellcollection.org/doi/pdf/10.1144/SP404.7, URL: https://www.lyellcollection.org/doi/abs/10.1144/SP404.7.

Shao, Q., Matthai, S., Driesner, T., Gross, L., 2021. Predicting plume spreading during CO₂ geo-sequestration: benchmarking a new hybrid finite element–finite volume compositional simulator with asynchronous time marching. Comput. Geosci. 25 (1), 299–323, Publisher: Springer.

Shariatipour, S.M., Pickup, G.E., Mackay, E.J., 2016. Investigation of CO₂ storage in a saline formation with an angular unconformity at the caprock interface. Petrol. Geosci. 22 (2), 203–210. http://dx.doi.org/10.1144/petgeo2015-039, URL: https://www.lyellcollection.org/doi/10.1144/petgeo2015-039.

Si, H., 2015. TetGen, a delaunay-based quality tetrahedral mesh generator. ACM Trans. Math. Software 41 (2), http://dx.doi.org/10.1145/2629697.

Sifuentes, W., Blunt, M.J., 2009. Modeling CO₂ storage in aquifers: Assessing the key contributors to uncertainty. In: Offshore Europe Oil & Gas Conference & Exhibition. Aberdeen, UK, p. 13. http://dx.doi.org/10.2118/123582-MS.

Span, R., Wagner, W., 1996. A new equation of state for carbon dioxide covering the fluid region from the triple-point temperature to 1100 K at pressures up to 800 MPa. J. Phys. Chem. Ref. Data 25 (6), 1509–1596. http://dx.doi.org/10.1063/1.555991, URL: http://aip.scitation.org/doi/10.1063/1.555991.

Stimpson, C.J., Ernst, C.D., Knupp, P., Pébay, P.P., Thompson, D.C., 2007. The Verdict Library Reference Manual. Technical Report, Sandia National Laboratories, URL: www.vtk.org/Wiki/images/6/6b/VerdictManual-revA.pdf.

Suter, E., Cayeux, E., Fris, H.A., Kårstad, T., Escalona, A., Vefring, E., 2017. A novel method for locally updating an earth model while geosteering. Int. J. Geosci. 08 (02), 237–264. http://dx.doi.org/10.4236/ijg.2017.82010.

Wellmann, F., Caumon, G., 2018. 3D structural geological models: Concepts, methods, and uncertainties. In: Advances in Geophysics, Vol. 59. Elsevier, pp. 1–121.

Wellmann, J.F., Finsterle, S., Croucher, A., 2014. Integrating structural geological data into the inverse modelling framework of iTOUGH2. Comput. Geosci. 65, 95–109. http://dx.doi.org/10.1016/j.cageo.2013.10.014, URL: https://linkinghub.elsevier.com/retrieve/pii/S0098300413002781.

Wilkinson, M., Haszeldine, R.S., Mackay, E., Smith, K., Sargeant, S., 2013. A new stratigraphic trap for CO₂ in the UK North Sea: Appraisal using legacy information. Int. J. Greenh. Gas Control 12, 310–322. http://dx.doi.org/10.1016/j.ijggc.2012.09.013, URL: https://www.sciencedirect.com/science/article/pii/S1750583612002277.

Zehner, B., Börner, J.H., Görz, I., Spitzer, K., 2015. Workflows for generating tetrahedral meshes for finite element simulations on complex geological structures. Comput. Geosci. 79, 105–117, Publisher: Elsevier.