* corresponding author: jerome.carrayrou@unistra.fr

Keywords

Geochemical modelling; instantaneous equilibrium chemistry; activity correction, outer fixed-point algorithm, Davis activity, Debye-Hückel activity, B-dot model.

Abstract

Reactive transport codes are very useful elements of environmental research. They now contain multiphysics with very complex algorithms, including flow, transport, chemical and sometimes heat transport, mechanical and/or biological algorithms. Because of this complexity, some parts of these algorithms still have not been sufficiently studied. Here, we present a comparison of 3 algorithms for activity correction, a specific subset of equilibrium chemistry algorithms. We show that the most used algorithm (the inner fixed-point algorithm) or the most rigorous algorithm (the full Newton) might not be the most efficient, and we propose the outer fixed-point algorithm, which is more robust and faster than other algorithms.

Introduction

The problem of groundwater management has received increasing attention, and many tools have been developed to address this issue. One of these tools, reactive transport modeling, was first limited to laboratory experiments [1] and then extended to the comprehension of problems in various fields [2]. Reactive transport modeling is actually a mature research field that has produced important results in many environmental domains, such as water management, sea water intrusion [3], long-term nuclear waste storage [4] and heavy metal contamination [5]. Numerous reactive transport codes are available, and some review articles [6–12] propose an overview of them. Examining these articles, it can be seen that all these codes include one or more activity correction models. Even though the different models of activity correction are usually well-detailed, the algorithmic method used to compute the activity coefficients and to incorporate these calculations into the entire chemistry algorithm is not given.
We show here that different algorithms can handle activity corrections and that they are not all equivalent. It seems that the most used algorithm, named the inner fixed-point algorithm, leads to numerical instabilities when handling highly-charged chemical species. We then recommend the new algorithm presented in this work: the outer fixed-point algorithm. It is more robust, faster and less sensitive to the initial conditions.

General concepts

A general formulation of a chemical reaction leading to the formation of one of the Nc species (Ci) from the number Nx of chosen components Xj is written as:
Instantaneous equilibrium chemistry is usually described using two fundamental concepts: mass conservation equations and mass action laws. According to the classical formulation stated by Morel and Morgan, mass conservation equations describe the conservation of the total concentrations of the components (Tj), and mass action laws describe the formation of each chemical species as a combination of the Nx chosen components.
Mass conservation equations are written using the species concentrations [Ci]:
On the other hand, mass action laws are written using the species {Ci} and components {Xj} of the activities:
To ensure the closure of the system, an activity model is used. The activity coefficient (i) is less than one and determines the species activity from its concentration.
Several activity models have been developed, all of which use the ionic strength of the solution (I):