An efficient numerical simulator for geothermal simulation A benchmark study

(cid:129) An e ﬃ cient numerical simulator (DARTS) for geothermal simulations and applications is proposed. (cid:129) A set of benchmark tests are performed in DARTS compared to state-of-the-art numerical simulators. (cid:129) The good matches with other simulators verify the capability of DARTS for geothermal simulation. (cid:129) Higher performance is achieved in DARTS owing to the Operator-Based Linearization (OBL) approach. Accurate prediction of temperature and pressure distribution is essential for geothermal reservoir exploitation with cold water re-injection. Depending on our knowledge about the heterogeneous structure of the subsurface, the reservoir development scheme can be optimized and the overall lifetime of the geothermal ﬁ eld can be extended. In this study, we present Delft Advanced Research Terra Simulator (DARTS), which provides fast and accurate energy production evaluation for geothermal applications. This simulation framework is suitable for uncertainty analysis with a large ensemble of models. In DARTS, we select the molar formulation with pressure and enthalpy as primary variables. Besides, the fully-coupled fully-implicit two-point ﬂ ux approximation on unstructured grids is implemented to solve the mass and energy conservation equations. For the nonlinear solution, we employ the recently developed Operator-Based Linearization (OBL) approach. In our work, DARTS is compared with the state-of-the-art simulation frameworks using a set of benchmark tests. We demonstrate that DARTS achieves a good match for both low- and high-enthalpy conditions in comparison to other simulators. At the same time, DARTS provides high performance and ﬂ exibility of the code due to the OBL approach, which makes it particularly useful for uncertainty quanti ﬁ cation in processes involving complex physics.

difficulties for the Cartesian grids to depict the geological structures accurately and require unstructured grids [5] to characterize the complexities.
Several simulators have been used in geothermal applications [6]. TOUGH2 [7] is the state-of-the-art simulator for general-purpose numerical simulation of multiphase fluid and heat transport in porous media. It is widely used for geothermal projects [8][9][10]. In TOUGH2, the natural formulation is implemented with pressure and temperature (or saturation) as primary variables. The IAPWS-IF97 of the International Association for the Properties of Water and Steam [11] is used to calculate water thermodynamic properties. AD-GPRS (Automatic Differentiation General Purpose Research Simulator) [12,13] is a powerful research simulation framework that also provides geothermal capabilities [14]. In AD-GPRS, both natural and molar formulations are implemented within the unified simulation framework, while the formulation used to calculate water and steam properties is the same as [15]. These frameworks provide capabilities for the prediction of geothermal development.
However, the complexity of physics and a large number of grid blocks in high-resolution geothermal models often challenge conventional simulation techniques. The complex physical processes (i.e., multi-phase flow, multi-component reactive transport) encountered in geothermal applications require robust, flexible, and efficient solutions. The governing properties can become highly nonlinear due to the complex behavior of fluid thermodynamic properties with respect to changes of pressure and enthalpy, especially when several phases exist in the system [16]. To accurately delineate the physical process happening underground, an advanced simulation strategy is necessary to improve the convergence of the nonlinear solution. Besides, large-scale reservoir simulation with multi-million control volumes is often needed to characterize and predict the behavior of a geothermal reservoir which slowing down the modelling process. Furthermore, to quantify uncertainties and optimize development strategies, a large ensemble of models are necessary to cover the wide range of parameter settings, which requires high-performance and reliability of forward simulation.
In this study, we present Delft Advanced Research Terra Simulator (DARTS) developed for various industrial applications [17,18]. DARTS includes capabilities for the solution of forward and inverse problems for subsurface fluid and heat transport. For the solution of highly nonlinear problems, the Operator-Based Linearization (OBL) approach is employed in DARTS. The OBL approach was proposed recently for generalized complex multi-phase flow and transport applications and aims to improve the simulation performance [19,20]. For spatial discretization, a finite-volume fully-implicit method in combination with two-point flux approximation on unstructured grids is implemented in DARTS. Besides conventional discretization in temporal and spatial space, DARTS also utilizes discretization in physical space using the OBL approach.
In the OBL approach, the nonlinear terms (i.e., accumulation, flux) in governing partial differential equations are discretized and written in the operator form depending on the physical state. State-dependent operators are translated into multi-dimensional tables in the parameter space. During the simulation process, state-dependent operators are evaluated at the required supporting points of the parameter space. Multi-linear interpolation is then applied to create a continuous description. This representation simplifies the construction of the Jacobian matrix and residuals since the complex physics calculations are translated into generic interpolation between supporting points, which are calculated adaptively [21]. As a result, the programming implementation is significantly simplified preserving high flexibility and performance of the code. Furthermore, the design of the simulation framework supports a further extension to the advanced parallel architectures, e.g., GPU [18,22]. In DARTS, the molar formulation is implemented with pressure and enthalpy [18,23] as primary variables for geothermal simulation.
To keep the high performance, essential cores in DARTS (e.g., linear solver, well controls, OBL interpolation, etc.) are programmed in C++. Different simulation engines for various physical processes (e.g., geothermal and compositional simulation) are implemented in a unified framework. To make the simulator flexible, C++ classes are exposed via a Python interface, which enables users to manipulate DARTS and easily control the simulation process. In this way, DARTS possesses both the performance of C++ and flexibility of the scripting language.
Here, we take a geothermal case as an example to demonstrate how the compatibility is reflected in practice. The Python interface provides DARTS with the capability to embrace complex properties describing specific physical phenomena. Besides the set of integrated geothermal properties implemented in C++, an open-source IAPWS-IF97 formulation is incorporated into DARTS by designing a wrapper around the open-source package in Python. Taking advantage of OBL, the incorporated physics from other sources are used to calculate supporting points while the derivatives are evaluated automatically during interpolation. Therefore, the flexibility of the Python interface provides DARTS with the extended capability to model various physical processes. At the same time, the main C++ routine guarantees the efficiency of the simulation.
The primary objective of this work is to validate DARTS with stateof-the-art simulators for geothermal applications. In this study, we assume that chemical interactions are not affecting the flow of mass and energy in the reservoir. Notice that DARTS framework is already extended for various chemical reactions [24] and their generalization for geothermal applications is ongoing. In this paper, we first briefly introduce the governing partial differential equations and the basics of the OBL approach. Next, single-component (water), single-and multiphase flow is incorporated with different models for benchmark tests. The solution and performance of DARTS are compared with TOUGH2 and AD-GPRS individually.

Methodology
In general, aqueous brine is used as the fluid for thermal circulation in geothermal development. For some applications, CO 2 [25,26] has been proposed as a heat carrier. In addition, minerals can be dissolved by the brine with a number of chemical reactions [24], making the fluid chemistry even more complicated, and hydrocarbon components can be mixed with brine and co-produced [19]. Such type of models requires a complicated equation-of-state (EoS) to describe realistic phase behavior.
To simplify the benchmark comparison, we start with the basic situation where only the water component exists in the studied system. Although only a single component is involved, both liquid and gaseous phases are present in high-enthalpy systems. In this case, the complex EoS of water is required for accurate characterization, as described in [11]. The large contrast in thermodynamic properties between liquid water and saturated steam should also be taken into consideration for the efficient simulation in high-enthalpy systems.
Here, we consider the governing equations and nonlinear formulations for two-phase thermal simulation with water, which can be described by mass and energy conservation equations: where: ϕ is the porosity, s p is the phase saturation, ρ p is the phase The saturation constraint is required to close the system: In addition, Darcy's law is used to describe the fluid flow in the reservoir, where: K is the permeability tensor k [mD], rp is the relative permeability of phase p μ , p is the viscosity of phase p p [Pa·s], p is the pressure of phase p [bars], γ p is the specific weight The rock is compressible, which is reflected by the change of porosity through: where: ϕ 0 is the initial porosity, c r is the rock compressibility [1/bars] and p ref is the reference pressure [bars].
Molar formulation [14,15] is taken as the system nonlinear formulation, in which pressure and enthalpy are chosen as the primary variables. The Newton-Raphson method, as a conventional approach, is usually adopted to linearize the nonlinear system of equations. The resulting linear system of equations on each nonlinear iteration can be expressed in the following form: where: J ω ( ) k is the Jacobian matrix defined at the k th nonlinear iteration.
The conventional linearization approach involves the Jacobian assembly with accurate evaluation of property values and their derivatives with respect to the nonlinear unknowns. The properties and their derivatives are usually based on either piece-wise approximations (such as relative permeabilities) or solutions of highly nonlinear systems using chain rule and inverse theorem [27]. Therefore, the nonlinear solver has to perform extra iterations to capture small variations in properties which are sometimes negligible because of the uncertainties in property evaluation and numerical nature of their representation. For this reason, we use the OBL approach, which helps to improve this behavior.

Operator-Based Linearization
In the OBL approach, the elements in separate terms (e.g., accumulation, flux, etc.) of the governing Eqs. (1) and (2) fully defined by the physical state ω can be grouped and represented by state-dependent operators. Pressure and enthalpy are taken as the unified state variables of a given control volume in geothermal simulation. Upstream weighting of the physical state is used to determine the flux-related fluid properties determined at the interface l. The discretized mass conservation equation in operator form reads: where ω n is the physical state of block i at the previous timestep, ω is the physical state of block i at the new timestep, Γ l is the fluid transmissibility. State-dependent operators are defined as Here, the phase-potential-upwinding (PPU) strategy [21] is applied in DARTS to model the gravity effect. The potential difference of phase p between block i and j can be written as: where: ω j is the physical state of block j at the new timestep, δ ω ( ) p is the density operator for phase p.
The discretized energy conservation equation in operator form can be written as: where: This agglomeration of different physical terms into a single nonlinear operator simplifies the implementation of simulation framework. Instead of performing complex evaluations of each property and its derivatives with respect to nonlinear unknowns, we can parameterize operators in physical space either at the pre-processing stage or adaptively with a limited number of supporting points. The evaluation of operators during the simulation is based on bi-linear interpolation, which improves the performance of the linearization stage. Besides, due to the piece-wise representation of operators, the nonlinearity of the system is reduced, which improves the nonlinear behavior [19,21]. However, to delineate the nonlinear behavior in the system, especially strong nonlinearity (e.g., at high-enthalpy conditions), it is necessary to select a reasonable OBL resolution to characterize the physical space. Too coarse OBL resolution may lead to large error in the solution [20].
A connection-based multi-segment well is used to simulate the flow in the wellbore [18]. The communication between well blocks and reservoir blocks is treated in the same way as between reservoir blocks. Besides, the top well block is connected with a ghost control volume, which is selected as a placeholder for the well control equations. The bottom hole pressure (BHP), volumetric and mass rate controls are available in DARTS to model various well conditions.
As for the BHP well control, the injector and/or producer will operate under fixed bottom hole pressure. A pressure constraint is defined at the ghost well block: The volumetric rate control in DARTS is implemented through volumetric rate operator ζ ω ( ) Similarly, the mass rate control is defined as: (14) where Q p mass is the calculated mass rate ⎡ ⎣ ⎢ ⎤ ⎦ ⎥ ζ ω kg/day , ( ) p mass is the mass rate operator. Any of the described well controls can be coupled with energy boundary conditions, defined by temperature or enthalpy of the injected fluid at the injection well. Since temperature is the function of pressure and enthalpy, it is expressed in operator form and the temperature well control reads: where T ω ( ) is the temperature operator, T target is the target temperature of the injected fluid. Alternatively, enthalpy of the injected fluid can be defined as: where h is the enthalpy of the well control block, h target is the target enthalpy of the injected fluid. For the production well control, enthalpy is taken equal to that of the upstream well block.

Benchmark test
In this section, we perform a set of benchmark tests and compare the simulation results of DARTS with state-of-the-art reservoir simulators TOUGH2 and AD-GPRS. The comparisons are performed with one-, two-and three-dimensional models. At the beginning of each comparison, the selected model with initial and boundary conditions is described, after which the simulation is performed and the result comparison is shown. Finally, we display the performance of different simulators.

One-dimensional case
In a geothermal reservoir, fluid is mainly transported convectively from injection well or influx boundary to production well. At the same time, heat is transported through convective and conductive flow, where convection usually dominates. However, conduction also plays an important role in the development of a geothermal reservoir as the main mechanism of re-charging cold re-injected fluid. Besides, the properties of reservoir fluid can be significantly affected by phase changes. For example, the difference in densities of liquid water and steam has a great impact on heat transport and fluid distribution. Accurate simulation of these mechanisms is necessary for both forward and inverse modelling.
Here, we validate the solution (i.e., pressure, temperature, saturation, etc.) of DARTS with solutions obtained using the two state-of-theart simulators TOUGH2 an AD-GPRS. Two simulation models (one horizontal and one vertical) are selected as benchmarks for one-dimensional comparison. Table 1 lists the parameter settings used in these two models. Fig. 1 shows the initial and boundary conditions of the horizontal model. This model is initialized with hot steam to mimic high-enthalpy geothermal reservoir conditions. The horizontal boundary in the X direction is set with free-flow condition while a no flow condition is assigned to the rest of the boundaries. By influx of cold water, the reservoir block containing hot steam will be cooled down and the steam condensation will be coupled with the flow. Three simulators are set with identical parameters and run with a similar simulation strategy (i.e., time-step, convergence tolerance, etc.).

Horizontal case
The solutions generated by different simulators are shown in Fig. 2. DARTS achieves a perfect match with TOUGH2 in pressure, temperature, and saturation solution. It accurately captures the thermal propagation with a sharp saturation front, behind which a two-phase transition zone reflecting the interaction between cold water and hot steam is observed. Fig. 2b displays the 'staircase' shape, which can be interpreted as the reflection of phase transition on the temperature profile. The first stair represents the two-phase transition zone where pressure and temperature are independent. Although there is a slight difference between the solution of DARTS and TOUGH2 vs. AD-GPRS, the mismatch is minor.

OBL convergence of 1D horizontal model
By performing the interpolation in physical space, OBL significantly reduces the computational resources needed for property calculation. However, reasonable resolution in physical space is necessary to accurately capture the nonlinearity of rock and fluid properties [19,20]. After a good match shown in the 1D horizontal case, the same model is used for sensitivity analysis of the influence of OBL resolution on the accuracy of the solutions. Since the model is initialized with high-enthalpy conditions and presents a complex phase transition process, it is expected to be more challenging for OBL to match the reference solution with a limited resolution due to the high nonlinearity of governing physics. Here, simulations with different OBL resolutions were performed and the corresponding solutions are compared in Fig. 3.
In Fig. 3, the red solid line represents the solution with 512-points of OBL resolution, which is the reference solution. The lowest 8-points resolution introduces the largest error in all of the solution profiles, which is because the coarsest resolution can barely capture the nonlinearity of physics. With the increase of OBL resolution, the solution approaches the reference results gradually.
Already with an OBL resolution of 16 points, DARTS closely matches the reference solution for the pressure profile, as shown in Fig. 9a. However, the temperature and saturation profiles in Fig. 9b and Fig. 9c demonstrate that it is still difficult to obtain accurate results with merely 16 points, and at least 128 points are needed to accurately capture the saturation shock, which indicates that the physical nonlinearity at high-enthalpy conditions heavily relies on the thermodynamic properties, and an accurate thermal solution (i.e., temperature or enthalpy) is essential for the representation of the full physical process in the geothermal reservoir. With an OBL resolution of 256 points, the solution overlaps with the reference line, which demonstrates that sufficient accuracy has been achieved. Fig. 4 shows the nonlinear iteration performed by DARTS at different resolutions and the relative linearization cost per nonlinear iteration in comparison with the reference solution. Here, the relative linearization cost represents the ratio of the CPU time between OBL and reference solution per nonlinear iteration. With the coarsening of OBL resolution, the total number of nonlinear iterations decreases: the coarser the resolution is, the more linear is the physical description and, hence, easier for the simulator to converge. Besides, the linearization cost per Newton step does not decrease much with the coarsening of the resolution. This can be explained by the fact that the time consumption for calculating the supporting points in the physical space only takes a small portion of the linearization process. However, the accuracy of the solution decreases with the resolution as it is shown in Fig. 3. A       reasonable OBL resolution should be selected to keep both accuracy and efficiency when dealing with highly nonlinear physics which can be easily tested for a simple 1D model as shown here. Fig. 5 shows the initial and boundary conditions of the vertical model. This model is initialized with cold water at the top grid cell while hot steam for the rest cells. The initial pressure is set as uniform throughout all grid cells. The top boundary is set with no flow condition while the free-flow condition is applied at the bottom. Due to the large contrast of thermodynamic properties between water and steam, the fluids will redistribute and reach equilibrium under the effect of gravity. The intention is to validate the capability of DARTS in dealing with gravity compared with other simulators.

Vertical case with buoyancy
Since liquid water is much heavier than steam, water flows downwards while steam rises up following the buoyancy effect. During the equilibrium process, heat residing in different phases is exchanged with the transport of fluids. Resulting from this thermal transport, phase transitions take place due to the large variation of enthalpy in liquid water and the steam phases. Besides, the pressure will be redistributed under the gravity effect, which influences the counterbalance of water and steam phase as well.
Figs. 6 and 7 display the evolution of water saturation and enthalpy profiles for different time steps as the simulation proceeds. Water   In the end, reversed enthalpy distribution is obtained with higher enthalpy at the top compared to the initial distribution. Fig. 8 shows the comparison of the result between DARTS, TOUGH2, and AD-GPRS. We can see that DARTS achieves a good match with TOUGH2 in all solution profiles. In Fig. 8a, the pressure curve consists of two parts with different slopes related to different fluid phase distribution (in Fig. 8c) among the grid cells. Again, there is a slight difference in the temperature curve in comparison with AD-GPRS, but the difference is minor and within the acceptable range.

OBL convergence of 1D vertical model
Similar to the 1D horizontal case, convergence analysis of the 1D vertical model at different OBL resolutions is performed. Strong nonlinearity is present in the system due to the co-existence of multi-phase flow, heat transfer and buoyancy. The result comparisons for pressure, temperature and water saturation are shown in Fig. 9. The red line represents the solution with 512-points of OBL resolution and is considered as the reference solution again. The solution obtained with the lowest 8-points OBL resolution is the farthest away from the reference one, which indicates this OBL resolution is too coarse to accurately delineate the highly nonlinear equilibrium process. With the increase of OBL resolution, the solution error starts to decrease and 256-points resolution provides an accurate solution. We noticed that the OBL resolution required for accurate solution for this case is the same with the 1D horizontal case, which reflects that these two cases represent a similar level of nonlinearity.

Two-dimensional case
Realistic geothermal reservoirs are usually heterogeneous. A large permeability contrast requires a robust numerical scheme. Besides, the initial condition of a geothermal reservoir can vary from low-enthalpy to high-enthalpy conditions, depending on the thermal gradient and depth of the reservoir, which may lead to significant variations with respect to thermodynamic properties of the in situ fluids. All of these uncertainties in the subsurface cause difficulties for reservoir simulation. Therefore, the capability of DARTS in dealing with realistic models under different initial conditions should be verified.

Case 1
In this part, a one-layer model extracted from a synthetic geological   model from the West Netherlands Basin -WNB [2] is chosen for the two-dimensional comparison. Fig. 10a displays the porosity distribution of the model, which ranges from 0.1 to 0.37. Since the geological model represents a fluvial system, we can see the channelized distribution of porosity. The dimension of the model is × × 60 40 1 with grid size of × × 30 m 30 m 2.5m as is shown in Table 1. A closed boundary condition is used in the 2D comparison. Both low-enthalpy and high-enthalpy initial conditions are selected for the comparison with TOUGH2 and AD-GPRS. In addition, different well controls for injection and production wells are employed to make the comparison more representative.
• Comparison of DARTS and TOUGH2 Table 2 shows the reservoir initial conditions and well controls used in validation with TOUGH2. The results are shown in Fig. 11a and b for low-enthalpy and high-enthalpy conditions respectively. The left column shows the TOUGH2 solution, which is taken as the reference one in the comparison. The right column displays the relative difference between DARTS and TOUGH2 solutions in pressure and temperature. A good match is observed in both pressure and temperature maps for both low-enthalpy and high-enthalpy conditions. The maximum relative temperature difference is around 1.6% for low-enthalpy conditions, while for the high-enthalpy scenario, the maximum temperature difference is around 3.5% in very few grid cells around the displacement front.
Since a no-flow condition is assigned at the boundary, the pressure gradient building up between injector and producer guides the direction of fluid flow. Besides, fluid tends to flow within the high permeable channels, due to the channelized distribution of reservoir properties. In the high-enthalpy case, because of the higher mobility of steam, the water-swept area is larger than that of the low-enthalpy case, even with  the same production scheme.
• Comparison of DARTS and AD-GPRS Table 3 shows the reservoir initial conditions and well controls used in the validation with AD-GPRS. Fig. 12a and b show the solution and the difference between DARTS and AD-GPRS under low-enthalpy and   high-enthalpy conditions respectively, where AD-GPRS solutions are considered as the reference. In Fig. 12, the left column shows the AD-GPRS solutions in pressure and temperature and the right column displays the relative difference between DARTS and AD-GPRS solutions. As it is shown, a good match is observed in both pressure and temperature map under both low-enthalpy and high-enthalpy conditions. The maximum relative temperature difference is around 3.0% for low-enthalpy conditions, while for the high-enthalpy case, the maximum difference is around 3.5% in a few grid cells.

Case 2
Next, a fracture network extracted from an outcrop imaging of the Whitby Mudstone Formation [28] is taken to run and compare simulation results. Discrete Fracture Model (DFM) [29] is selected to characterize the fracture network with unstructured grid discretization, see more details on the DFM discretization in [30]. The model dimension is of × × 1200 m 1700 m 1 m. The geometry and discretized grids of the fracture network are depicted in Fig. 13. The model consists of 6998 matrix cells and 1382 vertical fracture cells. The basic parameter settings used in this model are listed in Table 4.
The model is run in both AD-GPRS and DARTS with fixed injection rate and production BHP under low-enthalpy condition. Fig. 14 shows the results of AD-GPRS, which is taken as the reference, and the relative difference between AD-GPRS and DARTS. A reasonable match is achieved between the two simulators with a maximum relative pressure difference of 6% and maximum relative temperature difference of 1.6%. Because of the complexity of heterogeneous data pre-processing in TOUGH2 and some convergence issues in AD-GPRS at high-enthalpy conditions, this model is only used to compare with AD-GPRS under low-enthalpy conditions. Table 5 shows the parameter settings used in the 3D comparison.    Fig. 16 shows the pressure and temperature comparison of the selected 20 th layer between DARTS and AD-GPRS after 100 years of simulation. As is displayed, the thermal breakthrough has already been reached for the specified simulation time. A good match (<2.0%) is achieved in the pressure solution and the maximum relative temperature difference is about 2.0%, which can be seen as a close match as well.
We noticed the distribution of temperature error corresponds to permeability distribution. Higher permeability provides faster fluid flow and sharper temperature fronts causing larger differences.
To show the solution difference of each layer between DARTS and ADGPRS, l 2 norm is taken to calculate the relative difference layer by layer. The normalized difference of k th layer can be evaluated as follows, The pressure and temperature differences are plotted in Fig. 17. As is shown, the solution difference of each layer is pretty small (below 1.0%), which indicates a good match is achieved between solutions (see Fig. 18). Table 6 shows the performance of different simulators on the desktop Intel(R) Xeon(R) CPU 3.50 GHz. All runs have been performed in a single thread regime. On average, DARTS provides a better performance than TOUGH2 and AD-GPRS run at default parameters. A small timestep of 20 days is required in the high-enthalpy case for robust convergence. Since the timestep strategy in DARTS is different from TOUGH2, there is a slight difference in the total number of timesteps. The fast simulation in DARTS can be attributed to the OBL approach, which significantly simplifies the calculation of state-dependent properties and Jacobian assembly. A slightly higher number of nonlinear iterations in DARTS in comparison to AD-GPRS for low enthalpy cases is related to differences in convergence criteria.

Conclusions
Numerical simulations have been widely used for the evaluation and optimization of energy production from the subsurface including geothermal applications. In this paper, we show that the Delft Advanced Research Terra Simulator (DARTS) can be used for the prediction and optimization of heat production in geothermal projects. A set of benchmark tests were devised and utilized to compare the solutions of DARTS with TOUGH2 and AD-GPRS. Comparison in the 1D horizontal model verifies the capability of DARTS to capture sharp temperature and saturation shocks resulting from the large mobility ratio between saturated steam and liquid water. The convergence analysis of Operator-Based Linearization (OBL) resolution validates the approach.   Besides, it suggests a reasonable resolution in physical space for highenthalpy simulation with strong non-linearity in physics. Another 1D vertical model with buoyancy validates DARTS capability to model a buoyancy-dominated flow in high-enthalpy systems. The phase-potential-upwinding (PPU) strategy was adopted for the OBL approach, and a close match of the simulation results indicates the reliability of DARTS handling buoyancy-induced flow coupled with phase equilibrium. For the 2D model, the capability of DARTS to simulate planar fluid and heat transport in a heterogeneous fluvial system with different boundary and initial conditions is verified by the close match with both TOUGH2 and AD-GPRS. In addition, the multiple options of well controls (e.g., constant bottom hole pressure, constant rate with constant temperature or enthalpy) integrated in DARTS were checked in these comparisons. Finally, the 3D synthetic geological model comparison displays the ability of DARTS to simulate realistic geothermal fields. The performance comparison among the 3 simulators demonstrates that DARTS allows simulation with a noticeable reduction in CPU time owing to the OBL approach and advanced programming concepts.
In this study, we focus on the general capability of non-isothermal simulation in DARTS framework with single-component (water) and multi-phase physics. Some essential aspects (e.g., chemical reaction between fluid and rock, multi-component multi-phase non-isothermal flow) will be taken into account in the future research.

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.