Acidizing is a widely used technique for intensification of geothermal wells. At certain conditions, acid injection results in a formation of highly conductive channels, e.g. wormholes. The associated rock matrix dissolution is highly nonlinear and may influence the fluid flow as well as rock properties. Wormhole patterns and breakthrough time are mainly controlled by acid strength and injection rate. Due to high complexity and, therefore, high computational cost, existing conventional simulators usually do not include complex chemical reactions or use weak (sequential) coupling of flow with reactive transport. Although the sequential approach is good for various applications, it is limited by the CFL condition and might suffer from poor convergence. In this project, we investigate the wormhole phenomenon to provide a stable fully implicit method (FIM) to numerically solve the reactive flow and transport problem coupled with equilibrium and kinetic reactions associated with the acid injection. The developed framework is aimed to predict a pattern of dissolution as well as wormhole breakthrough time. For an accurate description of species interactions, we directly connect PHREEQC equilibrium computation with FIM framework for kinetic dissolution using the Element Balance approach. This coupling was achieved with an application of a recently introduced Operator-Based Linearization scheme utilized within a Delft Advanced Research Terra Simulator (DARTS).

Chemical Enhanced Oil Recovery (CEOR) methods can increase the oil recovery of a reservoir to more than 60% of the volumes originally in-situ. Typical oil recoveries range from 20% to 40% of OIIP under traditional primary and secondary recovery stages. The understanding of the mechanisms that control such complex processes is essential considering the continuously growing worldwide energetic demand. History matching of core flooding experiments data through numerical modeling is a powerful tool to understand the physical parameters and mechanisms that control the flow of fluids in novel recovery strategies as Foam-Assisted Chemical Flooding (FACF) in which the effects of oil/water interfacial tension generated by surfactant are combined with mobility control provided by foam.
This report presents the mechanistic modeling study of seven experiments on different core flooding techniques including continuous gas injection, Water-Alternating-Gas injection (WAG), surfactant injection and Foam-Assisted Chemical Flooding (FACF) on Bentheimer sandstone cores. CT scan images were acquired to monitor the evolution of saturation distributions, starting from the primary drainage stage. A 1D model was built for each experiment to history match the saturations revealed by the CT scan imaging through the cores at different pore volumes injected, pressure drops, production fractions, and oil recovery. Good agreement between simulation results and experiments was obtained, with minor mismatches in breakthrough times (0.04 ± 0.02 PVI) for some experiments.
For the history match of primary drainage and forced imbibition, the parameters of the Brooks-Corey equation for calculating relative permeabilities were obtained from a regression analysis on the available observed data. Actual residual saturations were lower than the average saturations reported at the end of each flooding stage.
History-matching of the continuous gas injection at waterflooding residual oil saturation and WAG at connate water saturation experiments demonstrated that the initial water saturation in secondary and tertiary gas injection strategies is a parameter that greatly controls gas mobility. Gas trapping and its effect on the flooding was confirmed for the WAG experiment.
The mechanistic modeling of surfactant injection at under-optimum salinity conditions showed that the geochemical modeling to determine initial ionic concentrations in the aqueous phase can be neglected in absence of naphthenic acids in the oleic phase, even when alkali is present in the solution. The shape and size of the oil bank formed after mobilization of residual oil are controlled by the proper numerical representation of the phase behavior in the oil/microemulsion system, whilst other features as oil bank velocity are controlled by surfactant adsorption. Relative permeability curves at low trapping numbers (forced imbibition) and high trapping numbers (under-optimum salinity conditions) were obtained.
Finally, history match of the drive foam in FACF led to the conclusion that the local equilibrium model for foam model implemented in UTCHEM may not be able to cover the entire range of possible mechanisms controlling the foam formation, specifically, the mechanism by which foam is generated due to increase in interstitial velocity and local pressure drop in zones of reduced effective porosity near a previously formed oil bank.

More than 30% of the world’s hydrocarbon reserves are located in carbonate reservoirs, and this percentage is likely to increase, as a result of discoveries of new giant oil fields in carbonate rocks, generically named “Pre-salt layers”. However, there are still some problems in understanding karst systems that still unresolved. The karst caves are one of the suitable analogs for karstic reservoirs that also spread all over the world. In this study, the mobile LIDAR (Light Detection and Ranging) data is used to characterize the geometries and to analyze the stability of tunnels under several depths from 5 caves in Bahia, Brazil. The studied caves are representing both of the karstification mechanism (Epigene & Hypogene). In general, there are two tunnel shapes among the caves: horizontal ellipse & Vertical ellipse shapes. Several factors could be controlled the shape origin of the tunnels, but from this study mainly caused by lithology or the geology structure factor. By comparing with the structural data, the conduit orientation generally shows the same trend. Therefore, these conditions suggest a geometrical correlation between the fractures and the caves, and that the observed fractures almost certainly acted as conduits for fluid flow. The stability analysis showed that the vertical ellipse tunnels are more stable compared to the horizontal ellipse. However, all of the tunnels are already unstable in shallow depth (on average less than 2 km depth). The sensitivity tests show several parameters that would affect the stability: Rock Properties (Rock Mass strength, Rock mass elastic, density), number of tunnels in the system and the distances between multiple tunnels.

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.

The decrease in the production of hydrocarbons in combination with a growing urgency to reduce carbon emissions drives a rapid study and application of geothermal energy in the Netherlands. The produced heat from the low enthalpy geothermal energy can be used for heating up the buildings and greenhouses. A doublet system is used in the geothermal scheme, which consists of one production well and one injection well. The reinjection is done to maintain the reservoir pressure and reduce pressure decline due to production, earthquake and subsidence prevention, and environmental safety. Thus, it is crucial to model the aquifer distributions appropriately as the injection and production well should be placed at the communicated sand bodies.
In this study, the architecture of the Alblasserdam Member in the Drechtsteden, West Netherlands Basin is modelled. An investigation of the paleo flow and the distribution of sand prone succession within the Alblasserdam Member is done by integrating seismic interpretation, reservoir sedimentology, and petrophysical evaluation. The data used consist of seismic data, well log dataset and core data from the acquired hydrocarbon exploration data. The Alblasserdam Member was deposited during the Late Jurassic to Early Cretaceous when the high tectonic activities occurred in the West Netherlands Basin and followed by series of inversion impulses made this member bounded by various members at the bottom and the top. The seismic response of the Alblasserdam Member shows a high thickness and depth variations. The sediments of the Alblasserdam Member are highly accumulated in the central part of the study area, bounded by two major faults striking SE-NW.
Several sections of well correlation are made based on cycles of changing accommodation to sediment supply ratio (A/S cycles). The correlation is supported with seismic interpretation, allow to map the distribution of sand-rich interval of the Alblasserdam Member in the Drechtsteden. Two potential aquifer intervals are found. A southwestward shifting of the main fluvial system is observed from the thickness trend of the shallower aquifer intervals. A petrophysical evaluation is made to analyze the properties of the potential aquifer that encompass net to gross thickness (N/G), average porosity, and average permeability. The deeper potential aquifer has a higher clean sand N/G, but lower average permeability, that might be caused by the compaction, and mineral precipitation after the inversion period, especially in the highly inverted area. In the end, a recommendation to place a pair of production and injection wells are proposed based on the high accumulation of sand-rich intervals, sufficient depth to produce heat >70^oC and absences of the fault. The wells are placed toward SE-NW trend which is parallel to the distribution of the main fluvial distribution and the orientation of the faults, with production well placed at the deeper level.

Simulation of multiphase flow in natural subsurface formations include selection of time-step size, i.e., the discrete snapshots (steps) over which the time-dependent nonlinear process is investigated. The simulation results would naturally depend on the size of the discrete time step, according to the order of accuracy of the implemented numerical scheme. For accurate analyses the required time steps might be needed to be very small, which then leads to very costly simulations, which are often times impractical in real-world simulation applications. Taking coarse-scale time-step sizes can reduce the computational costs, however, it results in lower simulation accuracy. Note that a common simulation practice is to take constant time-step everywhere in the computational domain. To resolve this challenge, this MSc thesis report presents a novel local time-stepping (LTS) method for multiphase flow in porous media. In this method, different time-steps can be applied in one computational domain; allowing for employing small time steps in sub-domains with higher sensitivity to the time step size and bigger time steps elsewhere. The LTS method is developed for all types of sequentially-coupled and fully-coupled simulation strategies. For all of these developments, the interfacial flux continuity is the key requirement to connect the sub-domains of different time zones. The fluxes chosen for the inner stages of the coarse time-step are consistent with the chosen simulation strategy. More precisely, in the case of the implicit-pressure-explicit-saturation (IMPES) method these fluxes are lagged in time (explicit), while for the sequential implicit and fully implicit (FIM) methods they are calculated at the current (implicit) time-step. For implicit transport method, the refined zone partition matrices are developed to algebraically form the refined zone Jacobian matrix and refined zone residual vector. The numerical results show that only a fraction of the domain can be simulated with small time steps; while the rest can be simulated at a bigger one. Local time-stepping method provides accurate simulation results compared with the fine-scale time-step simulation results for all three simulation strategies. Specially, the local time-stepping for FIM simulations, which is the mostly common commercial-grade simulation approach, can preserve the accuracy of the results and also the computational efficiency. For the studied cases, the computational efficiency gain of LTS for sequential implicit method is not as significant as that of the FIM. To exploit full adaptivity in simulation, as presented in Appendix B, one can directly apply the developed method of this thesis to adaptive implicit method (AIM); where sequential and fully implicit coupling approaches are applied in different sub-domains of the reservoir. Overall, LTS casts a promising approach to optimize accuracy-efficiency tradeoff when it comes to time step selection.

Reservoir simulation is a broadly used tool in the oil and gas production industry. It can be used for production history matching, performance forecasting and field development. A common Enhanced Oil Recovery (EOR) method is the injection of miscible gases. With the injection of gas into an oil reservoir, the system has to be considered as multiphase and multicomponent. Compositional simulation gives an accurate representation of the physics of multiphase and multicomponent systems. Solving the system at a high resolution with the compositional description is computationally very intensive. Simulating at a lower resolution with a coarse grid reduces the computational loads, but it becomes less accurate. Upscaling procedures are used to convert the high resolution, fine scale simulation parameters to lower resolution, coarse scale simulation parameters.

In this research, an upscaling procedure is examined that uses a thermodynamic non-equilibrium correction to keep a high accuracy for the coarse scale compositional simulation. This non-equilibrium phenomena is represented by a source term in the thermodynamic equilibrium equation. This source term is obtained from a fine scale simulation and upscaling the phase compositions. The source terms are tabulated against the composition and are used in the coarse scale simulation. An isothermal system with number components and number of phases is considered in this research. The capillary and diffusive effects are neglected as well as the chemical reactions and adsorption. For the fine scale, a 2D grid of 102x20 is used and the coarse scale grid is 1D with 102x1 cells. The fine scale simulation is run using Automatic Differentiation General Purpose Research Simulator (AD-GPRS) and the coarse scale simulation is run with an IMPEC simulator written in Matlab.

To quantify the accuracy of the coarse scale simulation a misfit value is calculated between the coarse scale saturation distribution and the upscaled fine scale saturation distribution.
Comparing the coarse saturation distribution with the upscaled saturation distribution shows that this method allows high accuracy results to be maintained with non-equilibrium upscaling of the compositional flow.

A correlation is made between the source term and composition and used in a coarse simulation instead of the tables created from the fine scale simulation. This reduces the accuracy compared to the fine scale, but it still gives a satisfactory improvement of coarse scale results. To improve the accuracy of the simulation using the correlation, an optimization procedure is performed. The optimization improved the accuracy to the level of the accuracy of the simulation using the full source term tables. With a correlation between the source term and composition, the global fine scale compositional solution is no longer necessary for the coarse simulation.