The interest in carbonate lithologies has risen in the last few decades, due to discoveries of giant and supergiant carbonate reservoirs. While studying the properties of this lithology, the existence of cave systems at great depths was highlighted, since their presence can lead to drilling related problems, as well as change the flow direction of fluids in the subsurface. The possibility of predicting their orien- tations, dimensions and mechanical stability in the subsurface would therefore be of great help in the hydrocarbon business, which is why in this study, thanks to data provided by a Lidar laser scanner, a dimensional analysis of caves’ conduits related to their position within the karst system is performed. The origin of caves can be hypogenic or epigenic and they can present different patterns of distribution of conduits. When deformation events cause the generation of orthogonal sets of joints in the rock layers and fractures intersect with each other and with bedding planes, the fluid flow in the subsurface becomes channelized into these features. The combination of this and speleogenetic processes leads to the development of maze caves, whose conduits’ dimensions vary within a system as a consequence of differences of continuity and connectivity along a fracture network. Variations in the shapes of caves’ conduits can be observed as well, depending on the supply of sediment and the position relative to the water table during dissolution. Mechanical deformation of karst systems can be simulated by using the linear elastic deformation mathematical model. In this way it is possible to study the influence that mechanical properties such as Young’s modulus and Poisson’s ratio have on the mechanical stability of caves in the subsurface.

A coupled description of flow and thermal-reactive transport is spanning a wide range of scales in space and time. Subsurface reservoir heterogeneity with complex multi-scale features increases the complexity further. The spatial resolution required during a simulation is dependent on the type and degree of geological heterogeneity, but also on the flowphysics and thewell locations (Karimi-Fard &Durlofsky, 2014). Traditional upscaling or multi-scale techniques are usually focused on the accuracy of the pressure solution and often ignore the transport. Improving the transport solution can however be quite significant for the performance of the simulation, especially in complex applications related to reactive and compositional flow. The use of a method called Adaptive Mesh Refinement (AMR) enables the grid to adapt dynamically during the course of the simulation, which facilitates the efficient use of computational resources (Karimi-Fard & Durlofsky, 2014; Cusini et al., 2016). In this work, the aim was to develop such an AMR framework for general-purpose reservoir simulation. The approach uses a multi-level grid and is applicable to fully unstructured grids. The adaptivity of the grid in the developed AMR procedure is based on a hierarchical grid. The multi-level grid is constructed starting with the static geological model, a fine-scale model that accurately approximates the prospective reservoir. The control volumes present in the model are described by a list of volumes, depths and reservoir properties, and by a connectivity list for all neighbouring blocks, with corresponding transmissibility. The coarser levels are constructed through cell partitioning of the fine-scale model, called level 0. A global flow-based upscaling is applied to redefine the cell properties at the coarser levels. Each coarser level is also described by a list of volumes, depths, reservoir properties and a connectivity list with corresponding transmissibility. An inter-level connectivity list is constructed which includes the connections between control volumes of consecutive levels. The connectivity lists with corresponding transmissibility, along with the cell properties of each control volume are combined to formglobal arrays describing the full hierarchy of levels. The simulation is then started at the coarsest level, while keeping the wells at the finest level of refinement for more accurate representation of the solution. Dynamic adaptivity is performed, based on adaptivity criteria which were developed specific to the used application, resulting in a new grid configuration. The multi-level connectivity list togetherwith the cell properties are redefined at each time step and used as input for the next simulation. The developed AMR framework was implemented in Delft Advanced Research Terra Simulator (DARTS, 2019) which uses the recently developed linearization technique called Operator-Based Linearization, developed by Voskov (2017). The performance of the proposed approach was illustrated for several applications, including geothermal energy extraction and reactive transport. For geothermal applications, several models were tested: a homogeneous reservoir with unstructured grid, a synthetic fluvial model with shale blocks illustrating the role of shale in thermal recharge, a highly heterogeneous fluvial system with low net-to-gross ratio, and several layers of a highly-heterogeneous SPE10 reservoir. The results of these flow examples demonstrate a high level of solution accuracy relative to the reference fine-scale solution. The solution error is considerably decreased as compared to the coarse-scale model. Important features are captured at a higher resolution while keeping the amount of cells to a reasonable number. For the fluvial model, a different approach was adopted for cell aggregation, consisting of grouping the facies together, where shale-regions are kept coarse, and highpermeability facies are kept at a fine level. This significantly improved the solution as compared to the coarse model, where cold front propagation is overestimated. For the application of reactive flow, the dissolution of calcium carbonate through the formation of wormholes was analyzed. Here as well, the solution is considerably improved when using the AMR model. The wormholes are accurately represented, with improved computational performance.