The numerical analysis of the physics of high temperature geothermal systems - which encompasses the physics of lower temperature systems - relies on the solution of a system of coupled nonlinear partial differential equations (PDE) and complementary equations of state (EOS) describing physical and thermodynamic properties over a very wide range of pressures and temperature. This system has many different characteristic timescales leading to severe stability constraints: permeability contrast between fracture networks and the 3D matrix, coupling of wellbore and reservoir models… This is the reason why flow simulations suffer limitations such as instabilities, repeated convergence failures, collapses of time step sizes…

Moreover, traditional numerical schemes fail to discretize properly this PDE system on the 3D unstructured meshes necessary to discretize accurately complex geological models. Much of current operational simulations are performed on Cartesian grids. Yet, simulating complex structures is intractable with such grids which do not allow the inclusion of the accurate position of surface discontinuities such as faults or geological interfaces and may lead to so-called grid effects. Beyond their inherent simplicity the use of orthogonal grids is linked to the fact that most of the commercial reservoir simulation codes use the two points finite volume scheme that require orthogonality constraints. This orthogonality is natural for grids. It is also an intrinsic property of Voronoï tessellations which may be used as an alternative.

Current software also suffers several limitations in terms of boundary conditions. Most of the time, there are only two possible boundary conditions types: fixed value/Dirichlet type for all primary variables or fixed fluxes/Neumann type for all conserved quantities. Mixed-type transient boundary conditions are not supported which impedes the convenient modeling of natural processes such as recharge or seepage or water table fluctuations. Workarounds may exist but are relatively tedious to implement and are not formulated in a generic way.

Transient complex upper or lower boundary conditions are mandatory to take into account some crucial processes. In volcanic island settings, the inland water table may be excessively deep and the interactions between the vadose zone and the fresh water recharge may hide geothermal resources. Though many groundwater simulation software can deal with the vadose zone, they are rarely designed to study multiphasic hydrothermal processes. Conversely, some geothermal reservoir simulators propose to take into account a gas/air component (e.g. EWASG module of TOUGH2, Battistelli, 2012) but they suffer from aforementioned numerical limitations. Regarding the bottom part of geothermal systems, when geothermal reservoir engineers set up a high temperature field model for production purposes, the intrusive heat sources are believed to be below the depth range of the model and are accounted for by choosing “appropriate” boundary conditions in the bottom layer of the model. This approach may seem incomplete from a scientific point of view if one wants to understand the whole hydrothermal system even if it is often considered sufficient for operational purposes.

Wells are obviously central features of geothermal exploitation as well as the coupling of wellbore heat and mass flow with the reservoir.