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.@en

We have been able to solve a reservoir simulation problem which was previously thought of as intractable: We simulated multiphase displacement, including viscous, capillary, and gravitational forces, for highly resolved and geologically realistic models of naturally fractured reservoirs (NFR) at the sector, i.e. kilometre, scale with very reasonable runtime. This has been possible because we used massive parallelisation and hierarchical solvers in conjunction with a new discrete fracture and matrix modelling (DFM) technique that is based on mixed-dimensional unstructured hybrid-element discretisations. High-resolution DFM simulations are important to resolve the non-linear coupling of small scale capillary - viscous and large scale gravitational - viscous processes adequately for sector scale NFR. Cross-scale process coupling in NFR controls oil recovery and NFR often exhibit power-law fracture length distributions, i.e. they do not possess an REV, and highly permeable fractures can extend over the full hydrocarbon column height. As a consequence, emergent displacement patterns have been observed which are difficult to quantify using traditional means of upscaling. However, such patterns could now be used as benchmarks to reach a better consensus on the correctness of promising new upscaling techniques. The parallel DFM technologies presented here allow us to to obtain these results much more efficiently and hence explore the parameter space in greater detail. We observed a linear scaling behaviour for up to 64 processes and a significant decrease in runtime when applying our parallel DFM approach to three highly refined NFR simulations. These contain thousands of fractures, up to 5 million elements, and have local grid-refinements below 1 m for model dimensions between I and 10 kilometres. We achieved this excellent speedup because we reduced inter-processor communication by minimising the overlap between individual domains and decreased idle time of individual processors by distributing the number of unknowns equally among the processors.@en

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.@en

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.@en

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.

@en

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.@en