Georg Kosakowski

OpenGeoSys: an open-source initiative for numerical simulation of thermo-hydro-mechanical/chemical (THM/C) processes in porous media

In this paper we describe the OpenGeoSys (OGS) project, which is a scientific open-source initiative for numerical simulation of thermo-hydro-mechanical-chemical processes in porous media. The basic concept is to provide a flexible numerical framework (using primarily the Finite Element Method (FEM)) for solving multifield problems in porous and fractured media for applications in geoscience and hydrology. To this purpose OGS is based on an object-oriented FEM concept including a broad spectrum of interfaces for pre- and postprocessing. The OGS idea has been in development since the mid-eighties. We provide a short historical note about the continuous process of concept and software development having evolved through Fortran, C, and C++ implementations. The idea behind OGS is to provide an open platform to the community, outfitted with professional software-engineering tools such as platform-independent compiling and automated benchmarking. A comprehensive benchmarking book has been prepared for publication. Benchmarking has been proven to be a valuable tool for cooperation between different developer teams, for example, for code comparison and validation purposes (DEVOVALEX and CO2 BENCH projects). On one hand, object-orientation (OO) provides a suitable framework for distributed code development; however, the parallelization of OO codes still lacks efficiency. High-performance-computing efficiency of OO codes is subject to future research.

GEM-Selektor geochemical modeling package: revised algorithm and GEMS3K numerical kernel for coupled simulation codes

Reactive mass transport (RMT) simulation is a powerful numerical tool to advance our understanding of complex geochemical processes and their feedbacks in relevant subsurface systems. Thermodynamic equilibrium defines the baseline for solubility, chemical kinetics, and RMT in general. Efficient RMT simulations can be based on the operator-splitting approach, where the solver of chemical equilibria is called by the mass transport part for each control volume whose composition, temperature, or pressure has changed. Modeling of complex natural systems requires consideration of multiphase–multicomponent geochemical models that include nonideal solutions (aqueous electrolytes, fluids, gases, solid solutions, and melts). Direct Gibbs energy minimization (GEM) methods have numerous advantages for the realistic geochemical modeling of such fluid–rock systems. Substantial improvements and extensions to the revised GEM interior point method algorithm based on Karpov’s convex programming approach are described, as implemented in the GEMS3K C/C+ + code, which is also the numerical kernel of GEM-Selektor v.3 package (http://gems.web.psi.ch). GEMS3K is presented in the context of the essential criteria of chemical plausibility, robustness of results, mass balance accuracy, numerical stability, speed, and portability to high-performance computing systems. The stand-alone GEMS3K code can treat very complex chemical systems with many nonideal solution phases accurately. It is fast, delivering chemically plausible and accurate results with the same or better mass balance precision as that of conventional speciation codes. GEMS3K is already used in several coupled RMT codes (e.g., OpenGeoSys-GEMS) capable of high-performance computing.

A parallel finite element scheme for thermo-hydro-mechanical (THM) coupled problems in porous media

Many applied problems in geoscience require knowledge about complex interactions between multiple physical and chemical processes in the sub-surface. As a direct experimental investigation is often not possible, numerical simulation is a common approach. The numerical analysis of coupled thermo-hydro-mechanical (THM) problems is computationally very expensive, and therefore the applicability of existing codes is still limited to simplified problems. In this paper we present a novel implementation of a parallel finite element method (FEM) for the numerical analysis of coupled THM problems in porous media. The computational task of the FEM is partitioned into sub-tasks by a priori domain decomposition. The sub-tasks are assigned to the CPU nodes concurrently. Parallelization is achieved by simultaneously establishing the sub-domain mesh topology, synchronously assembling linear equation systems in sub-domains and obtaining the overall solution with a sub-domain linear solver (parallel BiCGStab method with Jacobi pre-conditioner). The present parallelization method is implemented in an object-oriented way using MPI for inter-processor communication. The parallel code was successfully tested with a 2-D example from the international DECOVALEX benchmarking project. The achieved speed-up for a 3-D extension of the test example on different computers demonstrates the advantage of the present parallel scheme.

Analysis of field observations of tracer transport in a fractured till.

We analyze a set of observations from a recently published, field-scale tracer test in a fractured till. These observations demonstrate a dominant, underlying non-Fickian behavior, which cannot be quantified using traditional modeling approaches. We use a continuous time random walk (CTRW) approach which thoroughly accounts for the measurements, and which is based on a physical picture of contaminant motion that is consistent with the geometric and hydraulic characterization of the fractured formation. We also incorporate convolution techniques in the CTRW theory, to consider transport between different regions containing distinct heterogeneity patterns. These results enhance the possibility that limitations in predicting non-Fickian modes of contaminant migration can be overcome.

DIFFUSION OF Na AND Cs IN MONTMORILLONITE

The state and dynamics of water and cations in pure and mixed Na-Cs-montmorillonite as a function of the interlayer water content were investigated in the present study, using Monte Carlo and classical, molecular-dynamics methods. While highly idealized, the simulations showed that the swelling behavior of hetero-ionic Na-Cs-montmorillonite is comparable to the swelling of a homo-ionic Na- or Cs-montmorillonite. The mixed Na-Cs-montmorillonite is characterized by intermediate interlayer distances compared to homo-ionic Na- and Cs-montmorillonites. Dry, hetero-ionic Na-Cs-montmorillonite is characterized by a symmetric sheet configuration, as is homo-ionic Cs-montmorillonite.We found that at low degrees of hydration the absolute diffusion coefficient of Cs+ is less than for Na+, whereas at greater hydration states the diffusion coefficient of Cs+ is greaterthan for Na+. An analysis of the relative diffusion coefficients (the ratio between the diffusion coefficient of an ion in the interlayer and its diffusion coefficient in bulk water) revealed that water and Na+ are always less retarded than Cs+. With large interlayer water contents, tetralayer or more, Na+ ions preferentially form outer-sphere complexes. The mobility perpendicular to the clay surface is limited and the diffusion is equivalent to two-dimensional diffusion in bulk water. In contrast, Cs+ ions preferentially form ‘inner-sphere complexes’ at all hydration states and their two-dimensional diffusion coefficient is less than in bulk water.The question remains unanswered as to why experimentally derived relative diffusion coefficients of Cs+ in the interlayer of clays are about 20 times less than those we obtained by classical molecular dynamics studies.

Flow pattern variability in natural fracture intersections

We use numerical simulations to examine the variability of flow patterns in representative fracture intersection geometries. In contrast to existing studies of perfectly orthogonal intersections, we demonstrate that more realistic geometries lead to a rich spectrum of flow patterns. Moreover, numerical solutions of the Navier-Stokes equations in these fracture intersections indicate that non-linear inertial effects become important for Reynolds numbers as low as 1–100. Such Reynolds numbers often exist in naturally fractured formations, particularly in karst systems and in the vicinity of wells during pump tests.

An ab initio molecular dynamics study of hydronium complexation in Na-montmorillonite

The Car–Parrinello molecular dynamics simulation technique was used to predict the structure and dynamics of hydronium solvation in mono-, bi- and trihydrated Na-montmorillonite. In monohydrated montmorillonite, hydronium ions are located within the hexagonal rings of the basal clay plane. Oxygen sites of hydronium ions point towards the clay surface and hydrogen atoms towards the water layer. In bi- and trihydrated montmorillonite, hydronium ions form water-solvated, outer-sphere complexes. Similar to the solvation mechanism in bulk water, hydronium ions donate three hydrogen bonds to interlayer water molecules. In all studied hydration states, hydronium ions do not form hydrogen bonds with the basal oxygen sites. Similar to bulk water, the free energy barrier for a classical proton transfer between interlayer water molecules is of the order of kT and therefore not the limiting factor for the proton diffusion. The diffusivity of hydrogen in the interlayer is controlled by the structural rearrangements of the solvating water molecules.

Dissolution–precipitation processes in tank experiments for testing numerical models for reactive transport calculations: Experiments and modelling

In the context of testing reactive transport codes and their underlying conceptual models, a simple 2D reactive transport experiment was developed. The aim was to use simple chemistry and design a reproducible and fast to conduct experiment, which is flexible enough to include several process couplings: advective-diffusive transport of solutes, effect of liquid phase density on advective transport, and kinetically controlled dissolution/precipitation reactions causing porosity changes. A small tank was filled with a reactive layer of strontium sulfate (SrSO4) of two different grain sizes, sandwiched between two layers of essentially non-reacting quartz sand (SiO2). A highly concentrated solution of barium chloride was injected to create an asymmetric flow field. Once the barium chloride reached the reactive layer, it forced the transformation of strontium sulfate into barium sulfate (BaSO4). Due to the higher molar volume of barium sulfate, its precipitation caused a decrease of porosity and lowered the permeability. Changes in the flow field were observed with help of dye tracer tests. The experiments were modelled using the reactive transport code OpenGeosys-GEM. Tests with non-reactive tracers performed prior to barium chloride injection, as well as the density-driven flow (due to the high concentration of barium chloride solution), could be well reproduced by the numerical model. To reproduce the mineral bulk transformation with time, two populations of strontium sulfate grains with different kinetic rates of dissolution were applied. However, a default porosity permeability relationship was unable to account for measured pressure changes. Post mortem analysis of the strontium sulfate reactive medium provided useful information on the chemical and structural changes occurring at the pore scale at the interface that were considered in our model to reproduce the pressure evolution with time.

Decalcification of cracked cement structures

The benchmark problem presented in this paper deals with the leaching of calcium from hardened cement paste. The leaching of calcium results in the dissolution of the cement minerals which affects physical, chemical and mechanical properties of porous cement matrix. The dissolution of cement minerals in this case progresses heterogeneously as a consequence of a small-scale geometrical feature (crack) within a domain. Complexity of transport through cracked porous media combined with complex cement chemistry can lead to considerable modelling uncertainties. One possible way to get an insight into the robustness of modelling results is to perform benchmark based on (i) different transport models and solution methods (finite volume, finite element, etc.), (ii) different geochemical solvers and (iii) different coupling algorithms (sequential iterative and non-iterative). This benchmark is designed to gradually increase the complexity of the problem and in this way recognize modelling elements that are the most sensitive in terms of modelling results, e.g. evolution of physical and chemical properties. Five international teams participated in this benchmark exercise. The reactive transport codes used (HYTEC, MIN3P, OGS-GEM, ORCHESTRA, COMSOL Multiphysics-iPHREEQC) give similar patterns in terms of predicted concentrations of elements and the mineralogy. The level of agreement depends on the problem complexity related mainly to the weighting and conservation properties of different numerical methods, to the coupling between transport and reactive solver and the agreement of thermodynamic database.

Differential alteration of basaltic lava flows and hyaloclastites in Icelandic hydrothermal systems

Geological field observations evidence that active and fossil Icelandic hydrothermal systems are typically embedded into an intercalation of almost completely altered and nearly unaltered volcanic rock layers. We investigated the reasons for this finding with help of geochemical reaction path calculations, by studying the mineralogical evolution of contrasting lithofacies–basalt flows and hyaloclastites at various temperatures and pressures, different recharge water composition, and gas content. From this study, we conclude that the initial porosity of protoliths and volume changes due to their transformation into secondary minerals are sufficient to explain the different extents of alteration as observed in field studies. In addition, we present a generalized kinetic model to estimate the alteration time of glassy fragments in hyaloclastite as a function of grain size, surface roughness, and temperature. This time was found to be rather short, ranging from a few hours to a few years.