Project Topics
Get Complete Project Material File(s) Now! »
Coupling techniques: operator splitting and global approach
Systems (1.1) and (1.2) are decoupled, therefore they can be solved subsequently. On the contrary, transport and reaction equations are coupled. The approach chosen to solve these two groups of equations represents another watershed in reactive transport modeling. On one hand, it is possible to substitute reaction equations in system (1.2) and to solve a wide, heavy system of differential equations through for example Newton Raphson Method (global approach); on the other hand, it s possible to choose to sequentially (and sometimes iteratively) solve transport and reaction equations (operator splitting approach). A comprehensible introduction to both techniques can be found in (Steefel and MacQuarrie 1996).
Global implicit approach (GA) results harder to implement into a code but avoids coupling errors; although the power computing has grown significantly, the solution of the global system of differential equations may be heavy to solve for complex 3D problems. On the other hand, operator splitting (OS) approach engenders some coupling errors, but results more flexible and allows the coupling of different simulation tools (i.e. one code for the transport and another for the chemical reactions) chosen according to their respective strengths.
The differences between the two approaches and the different coupling techniques (SNIA, SIA, Strang…) for the OS approach will be detailed in Chapter 6, but it is important to understand that almost philosophical differences exist between the two approaches. During a reactive transport modeling workshop it occurred to me to hear the sentence « operator splitting is dead » and although people with such strong opinions exist, the superiority of one method with respect to the other is still to be proven and is strictly problem dependent (Fahs et al. 2008) and references therein.
OS approach opens to other debates that are specific to bio-geochemical solvers, especially related to the technique applied for modeling thermodynamic equilibrium. Two are the possible approaches: minimization of Gibbs Energy or Mass Action Law. The two approaches present strengths and liabilities, which are often linked to the phenomenon that is intended to be simulated; nevertheless, models and codes for biogeochemical reactions based on Mass Action Law constitute the majority.
A universe of reactive transport codes
A detailed study of the state of the art of reactive transport codes can be found in (Steefel et al. 2015), where twelve world-known codes are compared according to their capabilities and numerical methods. The review produced by Steefel and co-authors provides an overview on PHREEQC, HP1, PHT3D, OpenGeoSys, HYTEC, ORCHESTRA, TOUGHREACT, eSTOMP, HYDRO-GEOCHEM, CrunchFlow, MIN3P, PFLOTRAN. Some of the codes analyzed in the review started as mere tools for reaction simulations that were extended to one-dimensional transport (PHREEQC), others are ambitious platforms engineered to connect software for different chemical-physical phenomena (OpenGeoSys). These codes represent a percentage of all the reactive transport codes and reaction modules that are available to an interested user. Other examples are: CHEPROO (Bea et al. 2009), CORE2D (Samper et al. 2009), GEM-Selektor (Kulik et al. 2012), GEOCHEM-EZ (Shaff et al. 2010), GWB Geochemist s WorkBench (Bethke 2008), RETRASO (Saaltink et al. 2004), WITCH (Goddéris et al. 2006).
Regarding flow, ten over twelve codes revised in Steefel et al. (2015) are capable of tri-dimensional simulations, nine embed Richards equation for unsaturated flow, eight are able to deal with variable density flow and six with non-isothermal flow. Regarding transport, the most remarkable characteristic is that only one code over twelve has a full dispersion tensor, while of the eleven remaining ten have a simple diagonal tensor while TOUGHREACT doesn t use one. Nine codes are declared able to treat transport in multiple continua. From a geochemical point of view, extended Debye-Hükel and Davies activity models, non-isothermal geochemistry, surface complexation, ion exchange, aqueous gas exchange, kinetic mineral precipitation, Monod kinetics are considered the basics and implemented in nearly every code. Less common are Pitzer activity model and mineral nucleation, which are respectively embedded into four and six codes. Nine codes are declared to be able to treat mineral solid solutions.
Some extra-attention is due to the implementation of isotopes. Equilibrium and kinetic isotope fractionation are implemented in seven and eight codes, and theses implementations are somehow recent with respect to first versions of the codes (Wanner and Sonnenthal 2013; Druhan et al. 2013); one of the reasons for the growing interest in isotopes is the recent improvement of analytic techniques that allow the detection of isotope fractionation.
From a numerical point of view, only three codes analyzed in Steefel et al. (2015) apply a global approach (MIN3P, CrunchFlow, PFLOTRAN) and two of them (CrunchFlow, PFLOTRAN) allow to choose between global approach and operator splitting. Although code parallelization is becoming more and more frequent, allowing sensible reductions of the computational time, only the half of the codes exploited parallelization at the time of the publication of the review.
LHyGeS numerical models: a garden worth gardening
The patrimony of numerical models and codes at the Laboratory of Hydrology and Geochemistry at the University of Strasbourg is well documented. Codes have been developed for the simulation of mineral reactivity (KINDIS (Madé et al. 1994)), mineral reactivity and 1D transport (KIRMAT, (Gérard et al. 1998)), geochemical reactions (SPECY, (Carrayrou et al. 2002)), conservative transport (TRACES) and a number of other codes for more specific hydrological problems (Pan et al. 2015; Weill et al. 2017, Jeannot et al. 2018). Nevertheless, the evolution of geochemical models, numerical methods, and computing techniques is continuous and a constant update of these models is mandatory to ensure their survival.
Despite being a very powerful instrument, with numerous practical applications that continue at present time, KINDIS has two substantial limitations: it is not well suited for coupling with a multidimensional transport code and it is has been specifically designed for mineral reactivity. SPECY, on the other hand, has been conceived from the beginning as a tool with broader fields of application and to be coupled with transport codes. Nevertheless, it has never been exploited to its full potential, being used more as a laboratory for numerical tests than for the systematic simulation of actual geochemical problems.
In TRACES, who embeds the powerful numerical scheme of mixed and discontinuous finite elements, and who is very efficient for simulations of conservative transport, only very basic reactions were taken into account at the moment of its development.
The aim of this thesis is to develop a tool for the simulation of bio-geochemical reactions gathering all the capabilities of the existing models (and adding some, in particular, the treatment of isotopes) in a modern numerical environment. The tool is conceived as a flexible module to be coupled with other codes, primarily for the simulation of transport (TRACES) to produce a reliable and up to date tool for the simulation of reactive transport.
With paragraph 1.4 well in mind, one might question the usefulness of another reactive transport code. It is important to highlight that the point is creating neither clone nor a competitor of the already existing, mainstream reactive transport codes, which would be a quite ambitious goal, but to create an alternative which is transparent, controllable, adaptable, reliable and efficient to the maximum extent. Carrayrou et al. (2010), with the application of the MoMaS benchmarks, showed that different codes may still provide different results, thus underlying the necessity of plurality and intercomparisons. It is also useful to remind that although many of the codes mentioned before are available to the scientific public, only a few of them are Open source. Moreover, as clear as their user-guides can be, the comprehension and usage of these tools is not immediate (courses to learn how to model with some mainstream codes are often put in place); this implies that scientist performing experimental or fundamental studies are far from just downloading powerful software and running a simulation.
In conclusion, the tool is intended to serve as a complement to experimental activity and a solid basis for further numerical improvements; this manuscript is mainly directed to those potentially continuing the work, to provide an honest picture of reactive transport modeling at LHyGeS, its potential and limitations.
Structure of the work
The work is ideally organized into three main sections. The first section is focused on the reaction solver and begins with a preliminary study on the implementation of isotopes in reactive transport models (Chapter 2). An extensive bibliographic study was conducted on the topic, concluding that isotopes are implemented as autonomous species in all models taken into consideration. Chapter 2 also gives an overview of the growing interest in isotopes and explains their role (mainly as tracers) in reactive transport studies. Following conclusions of Chapter 2, no substantial modifications of the algorithm are required to implement isotopes in the model; on the other hand, the increase in the number of species and the potentially extremely low concentrations treated stress once more the importance of reliability and efficiency of computations.
Chapters 3 and 4 explore challenges in the solution of systems of equations arising from mathematical formulations of geochemical problems. Chapter 3 is dedicated to Thermodynamic equilibrium formulation. The method of primary and secondary species for Mass Action Law formulation is introduced as well as solution techniques of resulting algebraic non-linear systems. Liabilities of Newton Raphson method for the solution of non-linear systems are examined and solutions are proposed. In Chapter 3 are also listed and briefly explained the capabilities of the code at thermodynamic equilibrium.
Table of contents :
Chapter 1 – General introduction
1.1 Reactive transport
1.2 Governing equations
1.3 Coupling techniques: operator splitting and global approach
1.4 A universe of reactive transport codes
1.5 LHyGeS numerical models: a garden worth gardening
1.6 Structure of the work
Chapter 2 Implementing isotopes
2.1 Theoretical background
2.1.1 What is an isotope?
2.1.2 Notation
2.1.3 Isotopic abundance and its variations
2.1.5 Why do we care about isotopes?
2.2 Modeling isotopes
2.2.1 Modeling stable isotopes equilibrium fractionation
2.2.2 Modeling stable isotopes kinetic fractionation
2.2.3 Conclusions about modeling isotopes
Chapter 3 Thermodynamic equilibrium
3.1 Thermodynamic equilibrium solutions through a modified Newton Raphson method
3.1.1 Abstract
3.1.2 Introduction
3.1.3 Thermodynamic equilibrium: governing equations
3.1.4 Newton Raphson algorithm
3.1.5 Condition of the linear system
3.1.6 Working on a logarithmic base
3.1.7 Preconditioning
3.1.8 Scaling procedures in this work
3.1.9 Positive continuous fraction method
3.1.10 Numerical experiments
3.1.11 Numerical simulations: discussion
3.1.12 Conclusions about the strategies to improve Newton Raphson method
3.2 Thermodynamic capabilities of the code
3.2.1 Modeling surface complexation
3.2.2 Modeling ion exchange
Chapter 4 Mixed equilibrium and kinetics
4.1 Theoretical background and generic formulation
4.1.1 Generic formulation I
4.1.2 Generic formulation II
4.1.3 Systems of equations
4.2 Solving the systems of equations
4.2.1 Implicit and explicit, one-step or multistep methods of integration
4.2.1.1 Implicit and explicit methods
4.2.1.2 One-step and multi-step methods
4.2.1.3 Variable stepsize
4.2.2 An implemented explicit method: Richardson extrapolation of QSSA method
4.2.3 An implemented implicit method: BDF in DASPK
4.2.4 Solving systems with DASPK
4.2.4.1 Residual computation for DASPK 1
4.2.4.3 Residual computation for DASPK 2
4.2.4.3 Residual computation for DASPK 3
4.3 Numerical simulations
4.3.1 TST model, verification of results with PHREEQC and KINDIS
4.3.1.1 Description of the problem
4.3.1.2 Numerical simulations: results
4.3.2 Chilakapati test case: verification with publication
4.3.2.1 Description of the problem
4.3.2.2 Numerical simulations: results
4.3.2.3 Numerical simulations: effect of convergence criteria
4.4 Conclusions about mixed equilibrium and kinetics
Chapter 5 Solid solutions
5.1 Introduction and theoretical background
5.1.1 The interest in solid solutions
5.1.2 Theoretical background: Thermodynamics of solid solutions
5.1.3 Modeling solid solutions and their interaction with the aqueous phase
5.1.3.1 Equilibrium models for solid solutions
5.1.3.2 Kinetic models for solid solutions
5.1.3.3 Exploiting solid solutions concept for stable kinetic isotope fractionation
5.2 Numerical simulations of solid solutions
5.2.1 Modeling solid solutions at thermodynamic equilibrium
5.2.1.1 Verification with PHREEQC
5.2.1.2 Fe-Cr redox reaction, a reactive transport example
5.3 Conclusions about solid solutions
Chapter 6 Building SpeCTr, a reactive transport code
6.1 Coupling flow, transport and reaction
6.1.1 Governing equations
6.1.2 Global approach
6.1.3 Operator splitting
6.1.4 War of the approaches
6.1.5 Multicomponent reactive transport
6.2 TRACES
6.2.1 Code capabilities
6.2.2 Numerical schemes
6.2.3 Coupling with reaction module
6.2.4 SpeCTr
6.3 Validation: coupling and implementation of isotopes
6.3.1 Interest of validation
6.3.2 Presentation of the problem
6.3.2.1 Spatial discretization, flow characteristics and ground properties
6.3.2.2 Boundary and initial conditions, transport parameters
6.3.2.3 Reaction network
6.3.2.4 Considerations about Courant number
6.3.3 Results of numerical simulations: SpeCTr
6.3.3.1 Results: t(CFL), L = 0 m, no Cr fractionation
6.3.3.2 Results: t(CFL=1), L = 0 m, Cr fractionation
6.3.3.3 Results: reduced time step, L = 0 m, Cr fractionation
6.3.3.2 Results: t (CFL=1), L = 1.0 m, Cr fractionation
6.3.3.3 Results: reduced time step, L = 0.54 m, Cr fractionation
6.4 Conclusions about SpeCTr validation
Chapter 7 Application of SpeCTr: modeling Calcite dissolution & precipitation
7.1 From mixed flow reactor to column experiments and modeling: upscaling of calcite dissolution rate
7.1.1 Abstract
7.1.2 Introduction
7.1.3 Materials and experimental methods
7.1.3.1 Sample preparations
7.1.3.2 Aqueous solution preparations
7.1.3.3 Mixed flow reactor experiments
7.1.3.4 Column experiment
7.1.3.5 Aqueous sample analyses and thermodynamic calculations
7.1.3.6 Determination of calcite dissolution rate
7.1.4 Mathematical modeling of flow and reactive transport for the column experiments
7.1.5 Results and discussion
7.1.5.1 Mixed-flow reactor experiments
7.1.5.1.a Etching and etch pits morphology
7.1.5.1.b R- and R-G relationships as determined from VSI measurements
7.1.5.2 Column experiment
7.1.5.2.a Dissolution rates determined from VSI measurements
7.1.5.2.b Etch pit morphology
7.1.5.2c Comparison of the mean dissolution rates retrieved with VSI to those inferred from pit morphology
7.1.5.3 Modeled dissolution rates using 2D reactive transport simulations of the column experiment
7.1.5.4 Modeled dissolution rates: 1D versus 2D simulations.
7.1.5.5 Simulation of the Calcium breakthrough curve.
7.1.5.6 1D and 2D-reactive transport simulations of the column experiment: overview and perspectives
7.1.5.6a Mineralogical considerations
7.1.5.6.b Calcium breakthrough
7.1.5 Conclusion
7.2 Preparing Calcite dissolution rate modeling
7.2.1 Computation of reactive surface area for TST and SWM models
7.2.2 Time and spatial discretization
7.3 Mixing induced CaCO3 precipitation
7.3.1 Presentation of the test case
7.3.1.1 Spatial discretization, flow characteristics and ground properties
7.3.1.2 Boundary and initial conditions, transport parameters
7.3.1.3 Reaction network
7.3.1.4 Algorithms for porosity changes
7.3.2 Results of numerical simulations: SPeCTr
7.3.2.1 Results: constant porosity equilibrium and kinetic CaCO3 precipitation
7.3.2.2 Results: variable porosity – equilibrium CaCO3 precipitation
7.3.2.3 Results: variable porosity – kinetic CaCO3 precipitation
7.3.2.4 Conclusion about Calcite precipitation and porosity changes
7.4 – 3D Calcite dissolution modeling
7.4.1 Presentation of the problem
7.4.2 Results of 3D simulation
7.4.2 Conclusions about 3D simulation
