High technical and economic risks stemming from the lack of detailed knowledge of the subsurface hold back large-scale investments in geothermal energy. In a survey conducted on nine use cases from diverse geological settings across Europe and with different purposes (electricity/heating and cooling) and project objectives (scientific/commercial), we identify the “common practice” and the aspiration for the “state of the art” in geothermal exploration. For each use case, the survey investigates what workflows have been adopted and what data acquired by which methods at different stages of exploration. This provided a benchmark for exploration in a range geothermal play types. The survey shows that this industry-standard base-case can be adapted to improve exploration success and efficiency by (1) applying numerical modelling in early stages of exploration to guide strategic data collection, (2) novel application of innovative technologies and (3) closer integration of software tools for static geological interpretation and dynamic heat flow simulation.

To enable reliable exploration strategies for geothermal energy that have inherently lower economic and technical risks and hence increase public support, the multi-national, multi-disciplinary, and publicly funded FindHeat project is developing a novel, conceptual model-based geothermal exploration workflow. This workflow specifically focuses on faster turnaround times for exploration and appraisal of geothermal resources, making better use of legacy data and non-invasive geophysical techniques, and constraining uncertainties with respect to the size of the heat source and the range of possible heat production rates. Comprehensive social science research complements the technical work to set the foundation for new communication strategies that allow geothermal operators to earn the public trust that improved geothermal exploration and appraisal will lead to a more efficient and sustainable exploitation of geothermal energy. The workflow is being tested and validated at eight geologically diverse geothermal plays situated in Iceland, France, UK, Spain, and Netherlands, which allows us to demonstrate its economic and technical benefits as well as its societal impact.

Earthquake-triggered volcanic activity promoted by dynamic and static stresses are considered rare and difficult-to-capture geological processes. Calderas are ideal natural laboratories to investigate earthquake–volcano interactions due to their sensitivity to incoming seismic energy. The Campi Flegrei caldera, Italy, is one of the most monitored volcanic systems worldwide. We compare ground elevation time series at Campi Flegrei with earthquake catalogues showing that uplift events at Campi Flegrei are associated with large regional earthquakes. Such association is supported by (yet non-definitive) binomial tests. Over a 70-year time window we identify 14 uplift events, 12 of them were preceded by an earthquake, and for 8 of them the earthquake-to-uplift timespan ranges from immediate responses to 1.2 yr. Such variability in the response delay may be due to the preparedness of the system with faster responses probably occurring in periods during which the Campi Flegrei system was already in a critical state. To investigate the process that may be responsible for the proposed association we simulate the propagation of elastic waves and show that passing body waves impose high dynamic strains at the roof of the magmatic reservoir of the Campi Flegrei at about 7 km depth. This may promote a short-lived embrittlement of the magma reservoir's carapace otherwise marked by a ductile behaviour. Such failure allows magma and exsolved volatiles to be released from the magmatic reservoir. The fluids, namely exsolved volatiles and/or melts, ascend through a nominally plastic zone above the magmatic reservoir. This mechanism and the associated inherent uncertainties require further investigations but the new concept already implies that geological processes triggered by passing seismic waves may become apparent several months after passage of the seismic waves.

Thermohaline convection of subsurface fluids strongly influences heat and mass fluxes within the Earth's crust. The most effective hydrothermal systems develop in the vicinity of magmatic activity and can be important for geothermal energy production and ore formation. As most parts of these systems are inaccessible to direct observations, numerical simulations are necessary to understand and characterize fluid flow. Here, we present a new numerical scheme for thermohaline convection based on the control volume finite element method (CVFEM), allowing for unstructured meshes, the representation of sharp thermal and solute fronts in advection-dominated systems and phase separation of variably miscible, compressible fluids. The model is an implementation of the Complex Systems Modelling Platform CSMP++ and includes an accurate thermodynamic representation of strongly nonlinear fluid properties of salt water for magmatic-hydrothermal conditions (up to 1000°C, 500 MPa and 100 wt% NaCl). The method ensures that all fluid properties are taken as calculated on the respective node using a fully upstream-weighted approach, which greatly increases the stability of the numerical scheme. We compare results from our model with two well-established codes, HYDROTHERM and TOUGH2, by conducting benchmarks of different complexity and find good to excellent agreement in the temporal and spatial evolution of the hydrothermal systems. In a simulation with high-temperature, high-salinity conditions currently outside of the range of both HYDROTHERM and TOUGH2, we show the significance of the formation of a solid halite phase, which introduces heterogeneity. Results suggest that salt added by magmatic degassing is not easily vented or accommodated within the crust and can result in dynamic, complex hydrologies. We present a new numerical scheme for multiphase convection of salt water at magmatic-hydrothermal conditions based on the control volume finite element method. In a series of benchmarks with HYDROTHERM and TOUGH2, we find very good agreement of the simulated hydrothermal systems. Simulations at high-temperature, high-salinity conditions outside of the range of these models show the influence of solid halite on dynamic flow behavior and suggest that salt from magmatic degassing is not easily vented or accommodated within the crust.

The dynamic behavior of magmatic hydrothermal systems entails coupled and nonlinear multiphase flow, heat and solute transport, and deformation in highly heterogeneous media. Thus, quantitative analysis of these systems depends mainly on numerical solution of coupled partial differential equations and complementary equations of state (EOS). The past 2 decades have seen steady growth of computational power and the development of numerical models that have eliminated or minimized the need for various simplifying assumptions. Considerable heuristic insight has been gained from process-oriented numerical modeling. Recent modeling efforts employing relatively complete EOS and accurate transport calculations have revealed dynamic behavior that was damped by linearized, less accurate models, including fluid property control of hydrothermal plume temperatures and three-dimensional geometries. Other recent modeling results have further elucidated the controlling role of permeability structure and revealed the potential for significant hydrothermally driven deformation. Key areas for future reSearch include incorporation of accurate EOS for the complete H₂O-NaCl-CO₂ system, more realistic treatment of material heterogeneity in space and time, realistic description of large-scale relative permeability behavior, and intercode benchmarking comparisons.

High-resolution numerical simulations give clear insights into the three-dimensional structure of thermal convection associated with black-smoker hydrothermal systems. We present a series of simulations that show that, at heat fluxes expected at mid-ocean ridge spreading axes, upflow is focused in circular, pipe-like regions, with the bulk of the recharge taking place in the near-axial region. Recharging fluids have relatively warm temperatures. In this configuration, the system maximizes its heat output, which can be shown to be linked to nonlinearity in the fluid properties. Furthermore, we present a series of simulations with different permeability scenarios. These show that when permeability contrasts are moderate, convection maintains this pipe-like fluid flow structure. The permeability contrast has a dominant effect on flow patters only at early, immature, stages of convection, focussing upflow in high-permeability regions and downflow in low-permeability regions. In such early stages of convection, diffusive vent styles can emerge, which look remarkably similar to diffuse vent fields in natural systems. Finally, simulations in which permeability is defined as a function of temperature indicate that the brittle-ductile condition is likely to occur at temperatures not lower than 650°C. At lower brittle-ductile transition temperatures, the system cannot remove the heat delivered from the magma chamber and vent temperatures are substantially lower than 4000C. This result is in agreement with estimates of the brittle-ductile transition temperature from rock-mechanical studies and the occurrence of earthquakes in the oceanic lithosphère.

Field data-based simulations of geologic systems require much computational time because of their mathematical complexity and the often desired large scales in space and time. To conduct accurate simulations in an acceptable time period, methods to reduce runtime are required. A parallelization approach is attractive because fast multi-processor clusters are nowadays readily available. Here we report on our recent efforts to parallelize our multiphysics code CSMP ++ (Complex System Modelling Platform). In particular, we describe a parallel finite element-finite volume method for multi-phase fluid flow in heterogeneous porous media. We take a domain partitioning approach where the finite element mesh is partitioned into sub-domains, assigning each of them to a single processor. For each sub-domain a local finite volume mesh is constructed. We can now solve advection-dispersion type equations taking an operator splitting approach: Pressure diffusion is calculated with an implicit finite element method and advection with an implicit or explicit finite volume scheme. We have tested the accuracy, robustness and computational speedup of our new parallel scheme on a Linux cluster by means of three geologic applications. All tests give excellent computational speedup with increasing number of up to 32 processors. These results broaden the range of possible simulations in terms of spatial and temporal scale and resolution as well as numerical accuracy up to two orders of magnitude.

Realistic simulation of structurally complex reservoirs (SCR) is challenging in at least three ways: (1) geological structures must be represented and discretized accurately on vastly different length scales; (2) extreme ranges and discontinuous variations of material properties have to be associated with the discretized structures and accounted for in the computations; and (3) episodic, highly transient and often localized events such as well shut-in have to be resolved adequately within the overall production history, necessitating a highly adaptive resolution of time. To facilitate numerical experiments that elucidate the emergent properties, typical states and state transitions of SCRs, an application programmer interface (API) called complex systems modelling platform (CSMP ++) has been engineered in ANSI/ISO C ++. It implements a geometry and process-based SCR decomposition in space and time, and uses an algebraic multigrid solver (SAMG) for the spatio-temporal integration of the governing partial differential equations. This paper describes a new SCR simulation workflow including a two-phase fluid flow model that is compared with ECLIPSE in a single-fracture flow simulation. Geologically realistic application examples are presented for incompressible 2-phase flow, compressible 3-phase flow, and pressure-diffusion in a sector-scale model of a structurally complex reservoir.

We present a new finite element - finite volume (FEFV) method combined with a realistic equation of state for NaCl-H₂O to model fluid convection driven by temperature and salinity gradients. This method can deal with the nonlinear variations in fluid properties, separation of a saline fluid into a high-density, high-salinity brine phase and low-density, low-salinity vapor phase well above the critical point of pure H₂O, and geometrically complex geological structures. Similar to the well-known implicit pressure explicit saturation formulation, this approach decouples the governing equations. We formulate a fluid pressure equation that is solved using an implicit finite element method. We derive the fluid velocities from the updated pressure field and employ them in a higher-order, mass conserving finite volume formulation to solve hyperbolic parts of the conservation laws. The parabolic parts are solved by finite element methods. This FEFV method provides for geometric flexibility and numerical efficiency. The equation of state for NaCl-H₂O is valid from 0 to 750°C, 0 to 4000bar, and 0-100 wt.% NaCl. This allows the simulation of thermohaline convection in high-temperature and high-pressure environments, such as continental or oceanic hydrothermal systems where phase separation is common.

We present the benchmarking of a new finite element - finite volume (FEFV) solution technique capable of modeling transient multiphase thermohaline convection for geological realistic p-T-X conditions. The algorithm embeds a new and accurate equation of state for the NaCl-H₂O system. Benchmarks are carried out to compare the numerical results for the various component-processes of multiphase thermohaline convection. They include simulations of (i) convection driven by temperature and/or concentration gradients in a single-phase fluid (i.e., the Elder problem, thermal convection at different Rayleigh numbers, and a free thermohaline convection example), (ii) multiphase flow (i.e., the Buckley-Leverett problem), and (iii) energy transport in a pure H₂O fluid at liquid, vapor, supercritical, and two-phase conditions (i.e., comparison to the U.S. Geological Survey Code HYDROTHERM). The results produced with the new FEFV technique are in good agreement with the reference solutions. We further present the application of the FEFV technique to the simulation of thermohaline convection of a 400°C hot and 10 wt.% saline fluid rising from 4 km depth. During the buoyant rise, the fluid boils and separates into a high-density, high-salinity liquid phase and a low-density, low-salinity vapor phase.

We present new, accurate numerical simulations of 2D models resembling hydrothermal systems active in the high-permeability axial plane of mid-ocean ridges and show that fluid flow patterns are much more irregular and convection much more unstable than reported in previous simulation studies. First, we observe the splitting of hot, rising plumes. This phenomenon is caused by the viscous instability at the interface between hot, low-viscosity fluid and cold, high-viscosity fluid. This process, known as Taylor-Saffman fingering could potentially explain the sudden extinguishing of black smokers. Second, our simulations show that for relatively moderate permeabilities, convection is unsteady resulting in transiently varying vent temperatures. The amplitude of these fluctuations typically is 40 °C with a period of decades or less, depending on the permeability. Although externally imposed events such as dike injections are possible mechanisms, they are not required to explain temperature variations observed in natural systems. Our results also offer a simple explanation of how seismic events cause fluctuating temperatures: Earthquake-induced permeability-increase shifts the hydrothermal system to the unsteady regime with accompanying fluctuating vent temperatures. We demonstrate that realistic modelling of these high-Rayleigh number convection systems does not only require the use of real fluid properties, but also the use of higher order numerical methods capable of handling high-resolution meshes. Less accurate numerical solutions smear out sharp advection fronts and thereby artificially stabilize the system.