Accurate geothermal and chemical dissolution simulation using adaptive mesh refinement on generic unstructured grids

A coupled description of flow and thermal-reactive transport is spanning a wide range of scales in space and time, which often introduces a significant complexity for the modelling of such processes. Subsurface reservoir heterogeneity with complex multi-scale features increases the modelling complexity even further. Traditional multiscale 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 thermal-compositional flow. The use of an Adaptive Mesh Refinement enables the grid to adapt dynamically during the simulation, which facilitates the efficient use of computational resources. This is especially important in applications with thermal flow and transport where the region requires high-resolution calculations as often localized in space. In this work, the aim is to develop an Adaptive Mesh Refinement framework for geothermal reservoir simulation. The approach uses a multi-level connection list and can be applied to fully unstructured grids. The adaptivity of the grid in the developed framework is based on a hierarchical connectivity list. First, the fine-scale model is constructed, which accurately approximates all reservoir heterogeneity. Next, a global flow-based upscaling is applied, where an unstructured partitioning of the original grid is created. Once the full hierarchy of levels is constructed, the simulation is started at the coarsest grid. Grid space refinement criteria is based on the local changes and can be adjusted for specific models and governing physics. The multi-level connectivity lists are redefined at each timestep and used as an input for the next. The developed Adaptive Mesh Refinement framework was implemented in Delft Advanced Research Terra Simulator which uses the Operator-Based Linearization technique. The performance of the proposed approach is illustrated for several challenging geothermal applications of practical interest.


Introduction
Production development of prospective reservoirs includes the use of various technologies that provide information at many different scales. These scales range from core plugs being a few centimeters in size to well logs detecting properties a few meters around the well, and to seismic imaging covering a significant volume with limited resolution (few meters vertically and 10's of meters horizontally). However, time and capital limitations result in sparse direct sampling of reservoir rock and fluid properties. This is why the construction of reservoir models, through integration of these data using geostatistical reservoir description algorithms, has become a crucial step in resource development ( Branets et al., 2009 ). These algorithms conventionally result in fine-scale descriptions of reservoir properties (porosity, permeability) on grids of tens of millions of cells ( Christie, 1996 ). tion and often ignore the transport. However, in complex applications related to chemical and compositional EOR (Enhanced Oil Recovery), reservoir storage and geothermal industry, the number of conserved chemical species can be large which makes any improvement in transport solution quite significant for the performance and robustness of the simulation. A technique called Adaptive Mesh Refinement (AMR) provides an effective means for adapting the resolution of a model to solution requirements. This method is well developed in many areas of computational physics (e.g. fluid dynamics and solid mechanics) but is however not widely used for practical reservoir simulation ( Karimi-Fard and Durlofsky, 2014 ).
In today's literature, several researchers have developed and proposed AMR procedures to capture the local nature of transport processes. Bahrainian and Dezfuli (2014) have developed a novel unstructured grid generation algorithm which considers the effect of geological features and well locations in the grid resolution. This strategy involves the definition and construction of the initial grid based on the geological model, geometry adaptation of geological features and grid resolution control. Trangenstein (2002) used the combination of highresolution discretization methods with dynamically adaptive mesh refinement for a two-component single-phase model for miscible flooding. Pau et al. (2012) proposed an AMR algorithm for compressible two-phase flow in porous media. The method is implemented within a block structured adaptive mesh refinement framework which allows the grids to dynamically adapt to flow features and enables efficient parallelization of the algorithm. The coarse-scale permeability was obtained by averaging the fine-scale permeability. Similar techniques have been developed for compositional simulation ( Sammon et al., 2003 ), thermal problems ( Christensen et al., 2004 ), improved/enhanced oil recovery processes ( Van Batenburg et al., 2011 ), Discrete Fracture Networks ( Berrone et al., 2019 ), and many more applications.
In this work, the aim was to develop a dynamic AMR scheme using an unstructured multi-level gridding framework, for geothermal simulation in complex reservoirs. The focus lied particularly on thermalreactive flow and transport formulation which are required for a wide range of subsurface applications relevant to the energy transition including geothermal. Notice that heterogeneity plays a very important role in geothermal applications ( Shetty et al., 2018;Babaei and Nick, 2019 ). The geothermal doublet lifetime and heat recovery rate usually vary a lot with both reservoir parameters and operational management where uncertainties due to heterogeneity are dominating ( Willems and M. Nick, 2019 ). Heterogeneity in flow path and shale facies play an important role in water heat recharge which directly affects doublet performance at low net-to-gross ratio ( Crooijmans et al., 2016 ). Besides, complex heat extraction process and corresponding chemical interactions can also amplify the effect of heterogeneity ( Cui et al., 2016;Kala and Voskov, 2020 ).
As a starting point of our framework, a fine-scale geological model has to be constructed accurately approximating all reservoir heterogeneity. In reservoir simulation, this model is often represented by an array of volumes, depths and a connectivity list ( Lim et al., 1995 ) describing each control volume. Next, a global flow-based upscaling was applied and an unstructured partitioning of the original grid was constructed as suggested in ( Karimi-Fard and Durlofsky, 2014 ). This partitioning provides coarser levels of the original model which is also described by an array of volumes, depths and a connectivity list. A coarser connectivity list includes connections between control volumes at the given level as well as interconnections between the levels. Once the full hierarchy of levels is constructed, the simulation is started at the coarsest grid. Grid space refinement criterion is developed for particular applications. The multi-level connection list is reconstructed at each time step and used for the simulation. The proposed approach was implemented in Delft Advanced Research Terra Simulator (DARTS) ( Kala and Voskov, 2020;Wang et al., 2020 ).

Governing equations
General-purpose reservoir simulation is based on the solution of governing equations which describe mass and energy transfer of various species in the subsurface. The flow of mass and energy in a system with phases and components are described in this section. For this general-purpose thermal-compositional model, component mass conservation equations and a single energy conservation equation need to be solved ( Khait and Voskov, 2018b ). When chemical reactions occur in the system, an additional term describing kinetic reactions is added to the mass conservation equation ( Kala and Voskov, 2020 ). These governing relations are described as: where is the time, is the stoichiometric coefficient associated with kinetic reaction , is the rate of kinetic reaction . The right-hand side of the mass conservation Eq. (1) is the kinetic term which describes reactions. It is set to zero when no chemical processes are involved in the system. The rest of the terms in the system can be described as functions of spatial coordinate and/or physical state : ( , ) porosity, ( ) the mole fraction of component in phase , ( ) phase saturation, ( ) phase molar density, ⃖⃖⃖ ⃗ ( , ) phase velocity, ̃ ( , , ) source of phase , U ( ) phase internal energy, U ( ) rock internal energy, ℎ ( ) phase enthalpy, ( , ) thermal conduction.
An exception is the phase source term, which is also dependent on -well control variables.
The rock internal energy and thermal conduction are assumed to be spatially homogeneous for simplification of the problem, meaning that they are characterized as functions of the spatial coordinate only. Phase flow velocity ⃖⃖⃖ ⃗ is assumed to follow Darcy's law, expressed as: where ( ) permeability tensor, ( ) relative permeability of phase , ( ) phase viscosity, ( ) pressure in phase , ⃖ ⃗ ( ) gravity vector, ( ) depth (backward oriented).
The nonlinear unknowns in this system of equations are the pressure , the overall compositions of each component and the enthalpy ℎ .

Modeling approach
In order to solve the governing Eqs.
(1) and 2 , we apply a finitevolume discretization on a general unstructured mesh and perform a backward Euler approximation in time to both equations, where the phase velocities ⃖⃖⃖ ⃗ are substituted by the Darcy relation (3) : Here is the control volume for which the system is being solved, = ̃ is a source of phase , is the previous time step whereas + 1 is the time step we want to solve for. Capillarity and gravity are neglected in these equations, and a Two-Point Flux Approximation (TPFA) with an upstream weighting is applied. Δ , the phase potential, therefore simply becomes the difference in pressure between blocks connected via interface , while Δ is the temperature difference between these blocks; Γ = Γ ∕ is a phase transmissibility, where Γ is the constant geometrical part of the transmissibility (involving permeability and geometry of the control volume). Finally Γ = Γ is the thermal (conductive) transmissibility ( Khait and Voskov, 2018b ). This system of equations is solved for each mesh element in time, where the unknowns are the composition of the components and the pressure for the mass conservation Eq. (4) , and the pressure and enthalpy for the energy Eq. (5) .
In general-purpose reservoir simulation, the solving process requires the linearization of strongly nonlinear governing equations. In conventional reservoir simulators, a Newton-Raphson based method is typically used for the linearization, which solves on each nonlinear iteration a linear system of equations in the following form: where is the residual and is the Jacobian, which is the derivative of the residual with respect to primary nonlinear unknowns, defined at a nonlinear iteration . In this work, we use a recently-developed approach called Operator Based Linearization (OBL). The main idea of OBL is to transform the discretized mass and energy conservation Eqs. (4) and (5) to an operator form, where space-dependent and state-dependent properties of governing equations are separated. This provides the opportunity to approximate the representation of the exact physics of a problem through the discretization of the state-dependent properties. The underlying methodology of OBL is explained in more details in Voskov (2017) and Khait and Voskov (2018a,b) .

Connectivity list
The proposed AMR technique uses the Finite Volume Method (FVM) for discretization. The implementation of the finite volume discretization method to the mass conservation Eq. (1) requires the evaluation of the flow between two adjacent control volumes in terms of the cell pressures. Using a Two-Point Flux Approximation (TPFA), the flow rate is defined as: where: flow rate at interface of cells and , Γ phase transmissibility at interface of cells and , pressure of cell , pressure of cell . Indexing is based on a Cartesian structured mesh for simplicity.

Table 1
Connectivity list of the example grid from Fig. 1 .
To evaluate the flux between two adjacent control volumes, a socalled connectivity list is constructed, i.e. for each interface between two neighbouring control volumes, the indices of these cells are listed together with the transmissibility ( Lim et al., 1995 ). The result is a list with all connection pairs present in the grid. A few important points to be noted are: • Each connection consists of only two elements, • The connection pairs are not repetitive, • No-flow boundaries imply the absence of connections and are hence not listed in the connectivity list. Fig. 1 shows a simple example of a 2D Cartesian structured grid, with corresponding cell indexing. Table 1 shows its connectivity list. The list is expressed as two arrays, cell and cell , where each column represent a connection pair. Each pair has an associated interface transmissibility stored in the connectivity list.

Multi-level grid generation
The adaptivity of the grid in the developed AMR scheme is based on a hierarchical representation of connectivity list. The simulation grid is composed of several predefined levels representing the same geological properties at different resolutions. We start with a fine-scale model ( static geological model ) which accurately represents all reservoir heterogeneity. This grid is defined as level 0 and represents our finest level. The modeling grid is defined by a list of control volumes, depths, reservoir properties (all spatially distributed properties required to solve the discretized relations 4 and 5 ) for each mesh element, and a list of connectivity with corresponding transmissibility between neighbouring cells.
Next, level 1 is defined, where control volumes are constructed by aggregating fine grid cells. Upscaling is applied to redefine volume, depth and reservoir properties at a coarser level. A connectivity list, with corresponding transmissibility, is constructed for this level and inter-level connections are defined in addition. Similarly, more levels of coarsening can be constructed. A control volume in grid-level always consists Fig. 2. 2D Multi-level grid with three preconstructed grids (levels) with an example simulation grid which is constructed by aggregating control volumes from different levels. of cells from grid-level ( − 1) , resulting in a hierarchical relationship ( Karimi-Fard and Durlofsky, 2014 ). The simulation grid is then obtained by combining control volumes from grids of different levels. A schematic representation of this procedure is illustrated in Fig. 2 .

Cell aggregation
A mesh consists of a set of finite control volumes, each having vertices with allocated coordinates. To conduct cell aggregation, the centroid is first computed for each mesh element within the grid. Fig. 3 shows an example 2D unstructured grid to illustrate how cell aggregation is conducted. As can be seen, in this particular example, each cell has 3 vertices, and a centroid (represented in red) with coordinates and defined as ( , where and are the coordinates of the vertices. Each mesh element has an assigned index number. Cell aggregation is then carried out by dividing the grid in the -and -direction (and in the -direction for 3D models) into equidistant intervals Δ and Δ using a predefined scaling factor. Each interval has coordinates [ , + Δ ] in the -direction and [ , + Δ ] in the -direction. Centroids of cells whose coordinates are within a given -area are aggregated to form one coarse cell. To check whether a fine cell is within a given plane which will form coarse cell , the following algorithm is implemented for the coordinates and of the centroid of fine cell For the given 2D unstructured grid example in Fig. 3 , the so-called level 1 -i.e. the next level of coarsening -is shown in Fig. 4 . The numbers represent the assigned indices of the newly constructed coarse cells. If one wants to construct an additional level, the same procedure can be followed with a larger -and -range partitioning, where grid cells of level 1 are aggregated to form level 2.
For further steps into the generation of the levels, a list -"fines in coarse " -is constructed where the corresponding indices of the aggregated fine cells are listed for each coarse cell. Table 2 tabulates this list for the example above ( Figs. 3 to 4 ). This type of list is generated for each coarse level (level > 0 ) in the hierarchical grid. These lists are stored for the construction of the cell properties (e.g. volume, porosity) Fig. 3. 2D unstructured grid with centroids and with range partitioning (represented by the white lines) in the x-and y-direction with Δ and Δ spacing respectively. Aggregation is carried out for cells whose centroid fall within a given x-( [ ∶ + Δ ] ) and y-range ( [ ∶ + Δ ] ). E.g., all cells whose centroids are found within the yellow-highlighted 2D range are aggregated to form one coarse cell.  of the coarse levels, where the cell data from the fine level is needed during upscaling. Note that cell aggregation can also be conducted while taking care of highlighting geological features (e.g. fractures) and different facies in the model. For example, cell aggregation can be conducted by grouping domains with the same facies together into one coarse cell, or, in fractured reservoirs, by aggregating cells by isobar contours similar to Karimi-Fard and Durlofsky (2014) .
After cell aggregation is conducted, the connectivity list is then constructed describing all connections within each level and the inter-level connections. To illustrate the methodology, we use the simple structured grid from Fig. 1 , where cell aggregation was performed to form one coarse level.
In the proposed AMR scheme, the connectivity list of each level is determined systematically. Each mesh element consists of a set of vertices . E.g. a triangular mesh element comprises 3 vertices, and a Cartesian grid comprises 4 vertices. These vertices are numbered uniquely. The vertices comprised in a cell are stored in a list; this is done for each mesh element in level 0. To determine whether two control volumes and are adjacent, we take the intersection of both sets of vertices. That is: Each geometry has a different criterion. For 2D shaped mesh elements, the interface is a line; for 3D shaped cells, the interface is a plane. Hence the criterion is that the intersection length should equal 2 for 2D shapes and 3 or more for 3D shapes. This methodology is applied to the finest Table 3 Interfaces contained in each cell for level 0. The interfaces are expressed as connection pairs, describing the faces for each cell i .
Coarse cell I Faces
Next, the common faces between each coarse cell are determined. This is implemented by evaluating the intersection between the set of faces belonging to coarse cell and the set of faces forming coarse cell . This is expressed as: If a given coarse cell has one or multiple common faces with another coarse cell , these two cells form neighbouring blocks. For transmissibility computation in further steps, the area of the connecting interface is stored, which is here expressed as the sum of the intersecting fine grid faces.
For inter-level connections, a similar method is implemented. For each coarse cell in level , the intersection of its set of faces with the set of faces of a given fine cell is determined. This operation is conducted for every fine cell in level ( − 1) , except for the fine cells comprised in the evaluated coarse cell ( ∈ ). This is expressed mathematically as follows: Similarly, if a given coarse cell has a common face with a fine cell , the two cells are connected. This procedure is applied between all levels and ( − 1) . The result is a list of connections within level 0, a list of connections for each level , and an inter-level connectivity list, which describe the full hierarchical grid.

Transmissibility and upscaling
In this work, the AMR method is implemented for unstructured grids of any geometry. The definition of the transmissibility for unstructured grids is expressed as: where: Γ 12 transmissibility between cells 1 and 2, Γ 12 constant geometrical part of the transmissibility, mobility of a given phase , interface area, permeability of cell , distance between centroid of cell to interface area , unit vector normal to the interface, unit vector along the line joining centroid of cell to the center of interface A.
Here, the directional permeability of each cell is expressed as the magnitude of the cell's [ , , ] coordinates multiplied by the unit vector .
To solve the mass conservation Eq. (1) , the flow rate must be computed for the interface of every neighbouring cells. It is therefore necessary to compute the transmissibility for each dual connection listed in the connectivity list. The result is a list consisting of all connections, with their corresponding transmissibility. This methodology is applied at the finest level of refinement, level 0.
For thermal problems, another type of transmissibility Γ must be computed to approximate thermal conductive flux in the energy Eq. (2) .
Since thermal rock conduction is not as heterogeneous as permeability, the thermal transmissibility is defined as the geometric coefficient, that is, the area of the interface divided by the sum of the distances 1 and 2 from centroids to interface , multiplied by the average conduction 12 : As mentioned earlier, level 0 is represented by an array of volumes, depths and reservoir properties which are derived from the static geological model. Once the hierarchical grid is constructed, all cell properties must be redefined for the coarser levels (level > 0 ). This is done by upscaling the properties of the corresponding fine grid cells. The volume is upscaled by simply summing the volumes of the aggregated fine grid cells ; Depth upscaling is done by taking the average of the fine scale depths. The porosity, thermal conductivity, and rock heat capacity are upscaled using a volumetric averaging. For example, the sum of the porosity multiplied by the corresponding cell volume of each fine cell is taken over the total volume of the coarse cell ; In this study, for the upscaling of permeability, we use the flowbased upscaling technique developed by Karimi-Fard et al. (2006) ; Gong et al. (2008) ; Karimi-Fard and Durlofsky (2012) . This technique uses the pressure solution when the system has reached steady-state to compute the flow across each interface. The transmissibility can then be derived by rearranging the flow Eq. (7) . These approaches can be applied to unstructured coarse grids with generally-shaped control volumes ( Karimi-Fard and Durlofsky, 2014 ). The coarsening technique defines the coarse transmissibility Γ between two adjacent control volumes and . This is expressed as: The coarse-grid average pressures and , and the coarse-grid flow rate , are computed using a fine-grid pressure solution. These quantities are given by: where and define the fine-scale pressures in the corresponding coarse blocks. In the flow rate expression , indicates the interface between fine cells and and Γ denotes the transmissibility for this interface. This interface comprises a portion of the interface between coarse blocks and . For inter-level connections, a similar approach is used. For a given fine cell and coarse cell with interface , Eq. (18) is used with = , the pressure of fine cell , and the pressure of coarse cell . This procedure is conducted for each inter-level connection found within the hierarchical grid.
For thermal problems, a similar method can be implemented, but is however not computationally efficient as temperature takes significantly longer to reach a steady state. We therefore use Eq. (14) to compute the upscaled thermal transmissibility, where the area is expressed as the sum of the fine-scale faces which compose interface , and the distances and represent the distances between the cell centroid and the centroid of the coarse interface.

Global indexing
Once all needed parameters at every hierarchical level are evaluated, which include cell properties and a connectivity list with associated transmissibility for each level and between levels, it is necessary to combine the levels in order to form a global hierarchical set of grids. To combine the levels, it is however necessary to assign a unique indexing to each and every mesh element contained in the multi-level grid. This is where global indexing plays a role. For convenience, indexing is ordered starting from the coarsest level. An example of global indexing is shown in Fig. 6 .
The procedure used to assign global indexing is to simply add the number of cells of the previous level(s) to the current level. E.g., for level 1 from Fig. 6 , numbering starts at the total number of cells of level 2 2 ; for level 0, numbering starts at 2 + 1 . This procedure is applied to the bookkeeping lists such as "fines in coarse " and to the connectivity list of the corresponding levels. After the global indexing is applied to the connectivity lists, the connectivity list at each level and the interlevel connection lists are combined into one list. This is conducted by concatenating these lists to form one list describing all existing connections within the hierarchical grid. Regarding the array of cell properties, global indexing is applied by simply concatenating the arrays together in the right order, i.e. from the coarsest level to the finest level. This way, indexing is done in the same order as the global indexing. The result is a global array of volumes, depths and relevant reservoir properties describing each mesh element within the hierarchical grid. Having constructed the hierarchical grid and assigned it global indexing, the pre-processing stage is complete and the simulation with dynamic adaptivity can be performed.

Dynamic adaptivity framework
To determine whether grid adaptivity is necessary, we define refinement and coarsening criteria, which are dependent on the application used. In this study, we adopted an approach where the difference in solution variable is analysed between neighbouring blocks. Therefore, the difference in the solution variable of interest is computed between each pair of cells active in the simulation grid. If this difference is higher than a given threshold, both neighbouring blocks are refined. For the coarsening of a set of fine cells, belonging to a given coarse cell, the difference between all the corresponding fine cells and their neighbouring cells is computed; if each and every one of these connections have a difference in solution variable below a given threshold, the fine cells are coarsened to the next consecutive level.
For cells marked for refinement, the corresponding fine cells from the level below are added to the array of active blocks, which is used for implementation of the next time step, while the indices of the coarse cells in question are suppressed. Similarly, the cells marked for coarsening are suppressed from the active cells, and the corresponding coarse blocks are added. Fig. 7 shows an example of a two-level hierarchical grid. The current time step simulation grid is represented on the bottom left. After a check for adaptivity was conducted, cells 1 and 2 were marked for refinement. Hence as explained above, the cell indices 1 and 2 are suppressed from the array of active blocks, and their corresponding fine cell indices are added (6, 7, 10 and 11 for coarse cell 1, and 12, 13, 16 and 17 for coarse cell 2). The scheme at the bottom right of the figure shows the simulation grid which will be used for the next time step. As can be seen cell adaptivity results in an unstructured indexing..
Once the simulation grid is redefined and the array of active cells is updated, the connectivity list and corresponding transmissibility must be redefined. This is done by copying the list of connections for the whole hierarchical grid, where only the connections and corresponding transmissibility involving the active cells are kept, while connections involving non-active cells and their corresponding transmissibility are suppressed. Similarly, the same holds for the array of volume, depth For computation of the next time step solution +1 , the solution of the previous time step is required (see Eqs. (4) and (5) ). However, doesn't have the same grid configuration as the next time step +1 . It is therefore necessary to convert the grid of solution to the same configuration as the simulation grid at +1 . To do so, we use simple mapping techniques. A prolongation operator is firstly used to redefine the solution variable at each cell of the finest level of refinement (level 0). A so-called constant prolongation is implemented; i.e., all sub-domain values are set to the coarse value solution variable : Subsequently, restriction to the new simulation grid is conducted on the prolongated solution; i.e., for cells already at the finest level, the solution stays the same; when several control volumes are grouped into a single coarser control volume, the coarse value is set to the volumeweighted average of all sub-domain values ( Karimi-Fard and Durlofsky, 2014 ): A schematic representation of this procedure for the 2-level hierarchical grid and for the new simulation grid of Fig. 7 ( +1 ) is shown in Fig. 8 . The model, however, necessitates sequential numbering for mesh generation. It can be seen in Fig. 7 that indexing is non-consecutive when grid adaptivity is applied. This is where local indexing comes in play. That is, the active blocks indices are re-numbered in a sequential order to prevent undefined indices in the mesh. The global indexing is stored in a so-called Global to Local array for conversion back to the global indices for adaptivity check in the next time step. The described procedure, which redefines the grid configuration for the next simulation, is repeated at each time . It is also important to note that all previously computed operators in the OBL method are re-used after each successive timestep. This is possible since the parameter space for each state dependent operator in the OBL method is decoupled from any spatial property or discretization. This provides a significant speeds-up of the computation especially when simulation property is expensive to evaluate.
In the synthetic examples used to illustrate the performance of the AMR framework, the first time step simulation is started at the coarsest level. For improved accuracy, the cells containing the wells are kept at the finest level of refinement, level 0.

Applications for geothermal reservoirs
Geothermal technology has recently received substantial attention as an alternative source of energy. However, geothermal production systems have a relatively low return on investment, where uncertainties related to lack of detailed information about subsurface formations can significantly affect the quantification of the economic planning and feasibility of geothermal projects ( Willems, 2017 ). It is therefore important to reduce the uncertainty and produce a high accuracy solution while keeping the computational costs low. Geothermal systems therefore represent a good candidate for implementation of our AMR framework since it keeps the accuracy of simulation process close to the fine-scale while the performance is close to coarse-scale models.
Simulation of geothermal reservoirs implicates the solution of both mass (1) and energy (2) conservation equations where pressure and enthalpy are the solution variables. We are mostly interested in the accurate prediction of the temperature displacement front and resulting thermal breakthrough time. Dynamic adaptivity will be illustrated for 2 synthetic geothermal examples: • A homogeneous reservoir with unstructured meshing, • A heterogeneous fluvial system from Shetty et al. (2018) with low net-to-gross ratio.
In DARTS , the enthalpy is used as nonlinear unknown instead of the temperature. The adaptivity criteria are therefore applied to the enthalpy solution where the difference in enthalpy between two adjacent control volumes is analysed. This is done for each pair of connection within the simulation grid. Here, we applied the following adaptivity criteria: { if Δℎ > 70 , mark cells i and j for refinement, if Δℎ < 40 , ∀ ∈ , mark cells {∀ ∈ } for coarsening.
This adaptivity criteria is a simple heuristic and serves a practical purpose in this work. The proposed AMR method would greatly benefit from a more sophisticated criteria, for example, a criteria based on a posteriori error estimates similar to Vohralík and Wheeler (2013) ; Vohralík and Yousef (2018) .
The geothermal examples are illustrated by showing the fine-scale solution at different time steps versus the AMR solution and the coarsescale solution. Each synthetic example was analyzed quantitatively by conducting an error analysis where the error of both AMR and coarse solution are computed relative to the fine-scale solution. Both the L2 norm and L-infinity norm were calculated for each time step throughout the simulation. Moreover, to define the performance of the AMR method in terms of computational resources, the percentage of grid cells utilized in the simulation using the AMR grid, relative to the total number of cells in the fine-scale model was plotted for each example.

Case 1: homogeneous model
The first model is a simple 2D homogeneous reservoir (constant permeability) with unstructured triangular mesh. We consider a single injector (I) and a single producer (P) configuration. A two-level hierarchical grid is used, with 1420 cells in level 0 and 75 cells in level 1. Fig. 9 illustrates both levels, along with the permeability field (constant permeability of 2000 mD), and the well locations. The simulation parameters for this model are specified in Tables 6 and 7 in the Appendix.
The level 1 is illustrated in Fig. 9 where each color represents a coarse cell. As can be seen, cell aggregation was conducted by dividing the xand y-axes into 5 and 15 equidistant intervals. The cells at the well locations are kept fine at all times. The simulation was conducted for a period of 5500 days. The temperature solution at three different times is shown in Fig. 10 . Figure (a)    coarse-scale solution, and figure (d) shows the node distribution for the AMR simulation run.
The solution on the AMR grid demonstrates a particularly good match with the fine-scale solution. The node distribution shows high concentration along the front and at the well locations, and low concentration behind and ahead of the front, where no significant changes are observed. This considerably lowers the computational time as compared to running the fine-scale model. The coarse-scale solution differs notably from the AMR and fine-scale solution, with a faster cold front propagation at the coarse grid which is more pronounced in comparison at late times = 0 . 3 and 1. The relative error of the AMR solution is significantly lower than the coarse solution in both the L2 and L-infinity norm ( Fig. 11 ). Moreover, the number of cells is considerably reduced (see Fig. 12 ), ranging from 8 to 60%. The trend shows an overall increase as the front propagates, and a decrease when the cold front has reached the producing well, which results in coarsening at locations where no more thermal variations are detected. This considerably improves the performance of simulation since the AMR approach is much more favourable in terms of efficient use of computational resources (see Table 5 in section 6).

Case 2: sugar-cube shale model
Shales are often neglected in conventional reservoir simulation as the convective flow is never acquired in shales due to low permeability. For geothermal applications, they represent an important source of heat for thermal recharge of the cold water front. Modelling of the shales however significantly increases computational time since shales often occupy a significant amount of computational grid. Here, we test an application of our AMR approach to a sugar-cube model where cubes represents shale bodies and space between them fluvial channels. We use a simple 2D setup, with in total, a 5 by 6 shale block configuration. Shale blocks have a permeability of 10 −2 mD while the sand bodies have a permeability of 10 3 mD. The injector and producer are placed at opposite corners of the reservoir as shown in the Fig. 13 . Level 0 consists of 4588 fine cells. Level 1 is constructed differently from the conventional AMR approach with the sand channels kept at fine level, and only the shale blocks are coarsened by a ratio of 100 ( 10 × 10 ). The coarse grid, level 1, contains 1618 cells. The simulation parameters for this model are specified in Tables 6 and 7 in the Appendix. As can be seen on the AMR figure (b), the grid refines as soon as the cold front arrives at proximity to a shale body. The cold front is accurately represented on the AMR grid and there are no differences compared to the fine grid. On the coarse grid however, the cold front propagates further than for the fine and AMR model, which is clearly visible at the late time recording = 0 . 7 . When the cold front passes part of the shales blocks and these have cooled down, coarsening occurs as observed at = 0 . 7 . Fig. 15 depicts the error distribution through time of both the AMR and coarse model relative to the fine model.
As can be seen in Fig. 15 , the error of the coarse model is significantly larger than for the AMR model, where the error is close to zero. The high frequency changes in the error, especially observed in the ∞ norm, seem to correlate with refinement and coarsening of the mesh in between timesteps, similar to what was observed in Berrone (2010) . The percentage of cells used in the AMR grid relative to the number of cells used in the fine grid is shown in Fig. 16 . As can be seen, the percentage of cells ranges from roughly 35% to around 90% halfway through the simulation, when the cold front reaches the producing well, and then lowers to 65% when shale blocks proximal to the injector wells have cooled down to injection temperature, and hence coarsening occurs.
As observed, the computational time and effort was considerably reduced throughout the simulation, and the AMR solution outcome shows a very accurate representation of the fine-scale model (see Table 5 in section 6).

Case 3: fluvial heterogeneous model
Our AMR framework was tested for a heterogeneous reservoir with a low net-to-gross ratio (N/G = 35%). The permeability field ranges from 5 to 3400 mD with a significant amount of shale regions present. The hierarchical grid for this example is a structured grid and it comprises two levels. The finest grid, level 0, consists of 2400 grid cells with 40 cells in the x-direction and 60 cells in the y-direction. Level 1 was reduced to 150 mesh elements, where aggregation was done using 4x4 fine mesh elements, resulting in 10 grid cells in the x-direction and 15 grid cells in the y-direction. The permeability field along with the hierarchical grid for this example is shown in Fig. 17 . The location of the injector (I) and producer (P) are depicted in yellow on the permeability distribution. The simulation parameters for this model are specified in Tables 6 and 7 in the Appendix.
The simulation was conducted until cold water breakthrough reached the producing well. Fig. 18  The AMR mesh exhibits a significant improvement in temperature solution compared to the solution on the coarse grid. Refinement is mainly   focused at the front and slightly beyond the front, while areas where insignificant changes occur remain coarse. Important details, such as fingering effects at the cold water front, which are neglected on the coarse grid, are clearly visible in both fine and AMR solutions, which results in a more accurate representation of this physical phenomenon.
The relative error throughout the simulation run was recorded, where the fine model is taken as reference solution, for comparison be-tween the coarse and AMR model. Fig. 19 shows the L2 norm and the L-infinity norm error in time.
As can be seen, the marked improvement is also recorded in the error analysis, where the error between the coarse and fine model is notably larger than the error between the AMR and fine model. The L2 norm remains relatively constant for the AMR solution whereas it increases slightly in time for the coarse solution. The number of grid-cells used in the simulation ranges from 8 to 70% throughout the simulation (see Fig. 20 ). This represents a significant improvement in computational effort and time, while still capturing important features (see Table 5 in section 6).

Case 4: reactive transport
Carbonate reservoirs host a major part of the world's hydrocarbon reserves. But besides hydrocarbon reserves, the ongoing energy transition has resulted in an increase interest in geothermal systems where many are hosted by carbonate rocks. These reservoirs can have heavily fractured and karstified intervals, resulting in unforeseen hazards during drilling. Furthermore, naturally fractured carbonate reservoirs contain a large uncertainty in flow response due to the poor ability to predict the spatial distribution of discontinuity networks at reservoir-scale. Another important process related to dissolution is well acidization used to increase the production. This process involves the dissolution of reservoir rock to stimulate flow towards the wells. These chemical reactions are localized and form important features for accurate representation of the flow response. Furthermore, reaction rates which occur during dissolution are high, resulting in a sharp front in the flow response.
Moreover, during dissolution, formation and development of an unstable dissolution front with multiple wormholes can occur and its modeling is quite sensitive to the resolution ( Shaik et al., 2018 ). In nearwell acidization processes, the regime which forms a single dominating wormhole is the most preferable. It is therefore important to accurately predict this unstable dissolution while keeping the computational time reasonable. AMR is therefore a good solution to model these reservoirs and chemical processes to solution requirements.
In the flow example analyzed in this study, dissolution involves the following simple reaction where carbonate is dissolved: The kinetic reaction rate for this reaction is where is the mineral surface area, is the kinetic reaction constant, is the ion activity product, is the equilibrium product, and is the solid saturation. Permeability is updated using a power-law relationship defined as follows where 0 and 0 are the initial permeability and porosity, and is the power-law exponent. The proposed model simulates the phenomenon of unstable wormhole formation triggered by small perturbations in permeability. On one side of the reservoir, an injector well is placed which is perforated throughout the whole thickness. On the other side, the producer well is placed, also spanning the entire thickness of the reservoir. The model described in this example has dimensions of 100 by 100 meters. A constant permeability of 1 mD is used with 5% of random noise. The left illustration in Fig. 21 shows the well locations, along with the permeability of the reservoir. The hierarchical grid consists of two levels, where level 0 is an unstructured grid containing 2194 triangular cells. Cell aggregation was conducted to construct level 1, where the x-and y-axes were divided in 10 equidistant intervals, resulting in a grid with only 100 cells. Level 0 and 1 are shown in Fig. 21 . The simulation parameters for this model are specified in Tables 8 and 9 in the Appendix.
The AMR simulation was started at the coarse level, while keeping the well cells at the finest level throughout the entire simulation run. For this application, the adaptivity criterion is based on the solid composition,    As can be seen, the AMR solution is considerably more accurate than the coarse-scale solution. The far-propagating wormhole (at = 1), which is present in both the fine-scale solution and the AMR solution, is not well represented on the coarse-scale solution, where two extensive wormholes are present. The AMR solution however shows a very good representation of the fine-scale solution throughout time. The most extensive wormhole exhibits slight differences in thickness and some mi-  nor variations are observed at the other smaller wormholes. The node distribution follows the front, which is in this example quite dispersed, resulting in refinement spanning a wide area, especially at the last time step. However, considerable computational resources are saved at the early stage of the simulation.
To quantify the differences between fine-scale , AMR and coarsescale solutions, an error analysis was conducted. Here again, both the L2 norm and the L-infinity norm were computed for the AMR and coarse model, relative to the fine-scale solution. The graphs in Fig. 23 depict the outcome. As can be seen, the AMR error is once more significantly less than the coarse model, for both norms. For the L2 norm, the coarse-fine relative error is three times greater than the AMR-fine error at the final time step. The L-infinity norm of the coarse-fine error starts low at the first time step, where no extensive propagation is observed and where the model is close to the initial conditions, but then rapidly increases to 0.8 and remains more or less constant throughout. The L-infinity norm of the AMR-fine error seems to increase in time. This is due to the propagation of initially small errors in the solution. Note however that the relatively big error for both the AMR and coarse-scale model are not representative for this example and are related to another type of instability in the solution.
Similarly to the previous example, we have analyzed the total number of cells used in the AMR model, relative to the total number of cells contained in the fine-scale grid. The graph in Fig. 24 shows this quantity expressed in percentage.
As can be seen, the number of cells used during the AMR simulation is overall less than the number of cells present in the fine-scale model. Initially, the number of cells starts at 20%, which represents the use of  the coarsest level, with both left and right boundaries kept at the finest level. It then increases, fairly steeply, to around 95% due to the high injection velocity which corresponds to the dominating wormhole regime. Around the end of the simulation, the number of cells starts to decrease, which indicates coarsening at some locations. See Table 5 in the next section for the actual computational time. Although almost 100% of the cells is used at two thirds of the simulation, which is computationally expensive, considerable resources are saved in the beginning. Moreover, this problem is sensitive to the resolution which requires refinement at many locations in order to accurately capture the wormhole propagation.

Discussion
In the results section, all computational speed-up was indicated in terms of % of active cells w.r.t. fine-scale. In this section, the actual CPU times are highlighted for all cases. In order to have a fair comparison, the same nonlinear solver (Newton's based update with a fixed number of iterations), linear solver (direct one) and time-stepping strategy is used for the fine, coarse, and AMR runs. All these results are shown in Table 5 .
Even though the overhead of the AMR method is substantial, this can easily be explained by the non-vectorized Python implementation of the AMR procedure vs highly optimized Python and C++ implementation of the conventional simulation used in the coarse and fine simulation. Since the scope of this work is a proof of concept of the proposed AMR procedure prototyped outside of the simulation loop, our AMR code has not been optimized yet. This can be solved by either an application of Numba (just-in-time compiler for Python) or rewriting the procedure in C++. An expected speed-up, in our experience, is around two orders of magnitude when compared to the original Python implementation, thereby reducing the overhead to around 1.5% of the runtime of the AMR method and making it a viable strategy for geoscience applications. The framework presented in this paper is developed in the DARTS platform which can be used for a more general set of applications related to the energy transition. However, it is important to note the major differences in various energy applications. For example, in two types of applications shown in this study (geothermal and chemical dissolution cases), the coarser representation is still capable of accurately capturing important features of the geothermal dynamic behaviour. The coarsescale simulation in chemical dissolution, however, completely fails to represent the same dynamic behaviour (dissolution pattern) and effective characterization of the process (e.g. effective rock dissolution). It is therefore evident, as is also pointed out in the literature, that problems contained localized sharp gradients can greatly benefit from the AMR technique.

Conclusions
This study aimed at developing an Adaptive Mesh Refinement (AMR) technique in Delft Advanced Research Terra Simulator (DARTS) for general-purpose reservoir simulation. The developed AMR framework consists of a multi-level hierarchical grid, where levels are constructed through a mesh partitioning of the fine-scale model -the static geological model -which is represented by an array of properties (e.g. volume and porosity). The framework consists of the construction of the coarse levels through cell aggregation of the next consecutive fine level at the pre-processing stage. The method used to aggregate fine cells includes the grouping of subdomains whose centroids are found within a predefined 3D domain. In this study, domains are grouped by the partitioning of the x-, y-and z-axes into equidistant intervals. However, this strategy can easily be changed and improved.
The aggregation of the subdomains to form a coarser level is stored as an array of indices for the next stages, which consists of the indices of the fine cells comprised in its coarse control volume for each coarse cell.
Next, in order to solve the relevant governing equations, the flow must be computed at each interface present in the mesh. We, therefore, generate a list -called a connectivity list -describing all neighbouring cells within each level and between levels. The fine-scale transmissibility is then computed using the permeability field. Hereafter, a flow-based upscaling is applied in order to acquire the transmissibility of coarser levels and the inter-level transmissibility. Each control volume has defined parameters that are relevant for solving the system (volume, porosity, depth etc).
Once the hierarchy of levels is complete, the simulation can be started. Adaptivity check is performed at every time step, using criteria specific to the application. Once the regions for coarsening and refinement are defined, the solution is prolongated to the finest meshing level, and subsequently restricted from fine to the adaptive simulation grid. A new connection list and grid properties are constructed for the new coarsened schema. Once it is completed, the simulation runs for the next time step using the constructed simulation model.
The accuracy of the method was demonstrated for geothermal applications. Two models were tested, including a homogeneous model with unstructured gridding, a synthetic sugar-cube-like model with high permeability channels surrounded by shale blocks and a heterogeneous fluvial system model with a low net-to-gross ratio. High levels of solution accuracy relative to the reference fine-scale results are observed for both cases. An error analysis was conducted to record the differences between the AMR and the coarse solution relative to the reference finescale solution. The error resulting from the AMR model is significantly lower than for the coarse model, for all tested problems. The overall percentage of grid cells used in the AMR model relative to the fine-scale model is considerably decreased for most problems.
To conclude, the developed AMR method shows high levels of accuracy for both homogeneous and heterogeneous models and can be used for geothermal applications as well as for other applications implemented in DARTS. The number of cells in the AMR simulation, relative to the total number of cells of the finest level, is considerably reduced, which is very favourable in terms of efficient use of computational resources. The framework is applicable to two-and three-dimensional models and for unstructured as well as structured meshes. The applicability of the method to unstructured grids provides an effective means for solving complex geological systems.

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.

CRediT authorship contribution statement
Stephan