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.