# A projected dynamical system approach to mineral precipitation-dissolution reactions in geochemistry

Jocelyne Erhel, Bastien Hamlat, Tangi Migot

<figure><img src="image_1.png" /></figure>

## To cite this version:

Jocelyne Erhel, Bastien Hamlat, Tangi Migot. A projected dynamical system approach to mineral precipitation-dissolution reactions in geochemistry. 2024. hal-04631094

HAL Id: hal-04631094

https://inria.hal.science/hal-04631094v1

Preprint submitted on 1 Jul 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
A projected dynamical system approach to mineral
precipitation-dissolution reactions in geochemistry

Jocelyne Erhel¹, Bastien Hamlat¹,², and Tangi Migot¹

¹Inria and IRMAR, Rennes, France

²IFPEN, Rueil Malmaison, France

July 1, 2024

Abstract

In this paper, we study kinetically-controlled precipitation-dissolution geochemical reactions. Ordinary differential equations are very often used to model aqueous reactions. However, non smoothness is induced by full dissolution. Moreover, one reaction can include several minerals and one mineral can participate in several reactions. We develop a new model, which is a specific projected dynamical system, with constraints on time derivatives. To prove existence of a solution, we discretize the system with an explicit scheme and prove convergence of the approximate solution. Moreover, we solve the approximate differential system thanks to an alternating projection algorithm. Numerical computations on some examples illustrate this approach.

Introduction

Many environmental studies rely on modeling geochemical reactions such as CO₂ sequestra-
tion, water injection, and geothermal energy. In several of the aforementioned applications,
the chemical reactions are coupled with hydrodynamic processes leading for instance to a
reactive transport problem [23, 4, 16, 10, 7], studied also in the thesis [28, 6, 15]. It is of
great interest to get an in-depth knowledge of the geochemical properties involved as it is
very often the source of non-linearity or discontinuity in the system. In particular, this is
true when considering heterogeneous reactions such as precipitation-dissolution reactions
[27]. The difficulty here comes from the change of phase, which implies discontinuities in
the reaction rate [11]. We also refer the reader to the survey on computational methods
for geochemistry and reactive transport modeling in [22] for more insights.

In some of these applications, when the reaction rate is relatively slow compared to the
environment, the reaction does not necessarily reach an equilibrium. We refer to [16, 10,
7, 2, 9] for more details on equilibrium-controlled reactions. In this paper, we focus on kinetically-controlled reactions involving aqueous and mineral species. The concentration evolution of each species is governed by an ideal rate, depending on the law of mass action, while satisfying non-negativity and conservation constraints. As mentioned in [11], the resulting differential system might be discontinuous. This constitutes a real challenge for studying the properties of the system and designing numerical methods.

Despite the difficulty of the problem, the past decade has seen some new developments in the study of such geochemical models [16, 3, 1, 17, 21]. These new developments are following up on some early discussions in [11, 19, 29, 30, 31]. In the literature, the discontinuity induced by the phase change has been handled in several ways: an Heaviside function in [1, 21, 19, 29, 30, 31], a complementarity condition in [17, 20], and an ordinary differential equation with discontinuous reaction rate in [11, 3]. Equivalence between these three approaches is discussed in [17]. Moreover, in [16, 15], using a complementarity model, the authors study a semi-smooth Newton method to solve reactive transport problems with kinetic reactions. To the best of our knowledge, the papers mentioned above assume that there is at most one mineral per reaction, except for some examples in [11]. Recently, a model allows to consider several minerals in precipitation or dissolution in the same reaction for kinetically-controlled systems [14, 13]. It is a differential inclusion problem with Heaviside set-valued functions in the right-hand side. Numerical approximations are obtained thanks to a regularization and an implicit scheme.

In this paper, we propose another model, also for kinetic geochemical systems, with several minerals in one reaction and one mineral in several reactions. The right-hand side reflects directly physical constraints, limitations of the reaction rate and mass conservation. Our new model is a projected dynamical system, with a differential equation dealing simultaneously with the dynamics and the set of constraints.

The paper is organized as follows. Our projected dynamical model is defined in section 2, along with some assumptions. In section 3, we compare this approach to a complementarity model [17] and a differential inclusion model [14, 13], using two reference examples. The first one is one reaction with one mineral, whereas the second system has two reactions with three minerals [14, 13]. To deal with discontinuities, we develop in section 4 an explicit scheme by discretizing the differential equation and by using a projection-type method. We prove convergence of the numerical sequence, deducing the existence of a solution constructively. We define an alternating projection algorithm to compute an approximate solution and run some experiments to show numerical results in section 5. Finally, we discuss some avenues for future work.

# 1 A Projected Dynamical System Model

In this section, we first recall some characteristics of kinetic chemical systems. Then we define some constraints and introduce our new model.
## 1.1 Kinetic system

We denote by $n$ the number of species participating into $m$ reactions, where we assume that $n > m$. The vector $x \in \mathbb{R}^n$ contains the species concentrations. We consider chemical reactions composed of aqueous components and minerals. There is no restriction on the number of minerals per reaction on either side of the reaction. We denote by $W \subseteq \{1, \dots, n\}$ the set of indices of the aqueous species.

$S \in \mathbb{Z}^{n \times m}$ denotes the stoichiometry matrix, where, by convention, the stoichiometric coefficients of the reactants are taken as negative integers. The sign convention is required here to define the rate of each reaction.

We assume that the activity of an aqueous species is equal to its concentration, whereas the activity of a mineral species is equal to 1. In these conditions, we use the kinetic mass action law to define the reaction rates $\tau(x) \in \mathbb{R}^m$. For each $j = 1, \dots, m$, the reaction rate is written

$$ \tau_j(x) = k_{r,j} \prod_{\substack{i \in W \\ s_{i,j} < 0}} x_i^{-s_{i,j}} - k_{p,j} \prod_{\substack{i \in W \\ s_{i,j} > 0}} x_i^{s_{i,j}}, \quad (1) $$

where we denote by $k_r > 0$, resp. $k_p > 0$, the constant vector of the left side, resp. right side, of reactions corresponding to the reactants, resp. products.

The table below summarizes notations used throughout the paper.

<table><thead><tr><th>Name</th><th>Definition</th></tr></thead><tbody><tr><td>$n \in \mathbb{N}$</td><td>Number of species</td></tr><tr><td>$m \in \mathbb{N}$</td><td>Number of reactions</td></tr><tr><td>$S \in \mathbb{Z}^{n \times m}$</td><td>Stoichiometry matrix</td></tr><tr><td>$x \in \mathbb{R}^n$</td><td>Concentrations</td></tr></tbody></table>

<table><thead><tr><th>Name</th><th>Definition</th></tr></thead><tbody><tr><td>$Q \in \mathbb{R}^{n \times (n-m)}$</td><td>Conserved quantities matrix</td></tr><tr><td>$k_p \in \mathbb{R}^m$</td><td>Products species constant</td></tr><tr><td>$k_r \in \mathbb{R}^m$</td><td>Reactants species constant</td></tr><tr><td>$\tau \in \mathbb{R}^m$</td><td>Reaction rates</td></tr></tbody></table>

## 1.2 Constraints

Clearly, the concentrations $x(t)$ must be nonnegative. As long as $x(t) > 0$, reaction rates drive the dynamics of the system, so that $\dot{x}(t) = S\tau(x(t))$. Reactions are perturbed, but not necessarily stopped, when a mineral species is vanishing. We introduce constraints to model reactions in that case.

The first constraint concerns the concentrations, which must satisfy conservation rules. Let $Ker(S^T)$ be the kernel of $S^T$ and $Q$ a basis of $Ker(S^T)$. Let $x_0 \ge 0$ and $C_0$ be defined by linear inequality and equality constraints:

$$ C_0 = \{x \in \mathbb{R}_+^n : Q^T(x - x_0) = 0\}. $$

The set $C_0$ is non-empty and convex. The following lemma links the so-called conservative variables $Q$ to the derivative of $x$.
**Lemma 1.1.** Let $x: [0, +\infty[ \to \mathbb{R}^n$ be an absolutely continuous function such that $x(0) = x_0$ with $x_0 \ge 0$. For almost every $t \in [0, +\infty[$, it holds true that

$$Q^T(x(t) - x_0) = 0 \iff Q^T \dot{x}(t) = 0.$$

*Proof.* The function $y(t) := Q^T(x(t) - x_0)$ is also absolutely continuous with $\dot{y}(t) = Q^T \dot{x}(t)$ almost everywhere and $y(0) = 0$.

Then, using derivation properties, we get, almost everywhere, $y(t) = 0 \iff \dot{y}(t) = 0$ and the desired result. □

As noted previously, if $x(t) > 0$ then $\dot{x}(t) = S\tau(x(t))$ so that $Q^T \dot{x}(t) = 0$ and, thanks to Lemma 1.1, $x(t) \in C_0$. In our model, the concentration vector must always satisfy the constraint $x(t) \in C_0$.

The second constraint concerns the time derivatives when a mineral vanishes. We introduce the tangent cone to $C_0$ at $x \in C_0$:

$$\mathcal{T}(C_0; x) = \{d : Q^T d = 0, d_i \ge 0 \text{ for all } i \text{ such that } x_i = 0\}. \quad (2)$$

It is non-empty, closed and convex. We observe that, if $x(t) > 0$, then $\dot{x}(t) \in \mathcal{T}(C_0; x(t))$. In our model, the time derivative $\dot{x}(t)$ must satisfy $\dot{x}(t) \in \mathcal{T}(C_0; x(t))$ at any time.

Finally, we introduce a third constraint, connecting $\dot{x}(t)$ and $S\tau(x)$. Let $D(x)$ be defined by

$$D(x) = \{\text{diag}(S\tau(x))\lambda \text{ with } \lambda \in [0, 1]^n\}^1.$$

The set $D(x)$ is not empty, closed and convex. Clearly, if $x(t) > 0$ then $\dot{x}(t) \in D(x(t))$. In our model, the time derivative must satisfy $\dot{x}(t) \in D(x(t))$ at any time.

## 1.3 Projected Dynamical Model

We recall that the projection of a vector $v$ onto $C$ is defined by

$$\mathrm{Proj}(C; v) := \{w \in C : \|v - w\|_2 = \inf_{z \in C} \|v - z\|_2\}.$$

Combining the conditions described above, we define a specific projected dynamical system:

$$\left\{ \begin{array}{l} x(t) \in C_0, \\ \dot{x}(t) = \mathrm{Proj}(\mathcal{T}(C_0, x(t)) \cap D(x(t)); S\tau(x(t))), \\ x(0) = x_0. \end{array} \right. \qquad (3)$$

This model is similar to a classical projected dynamical system, with a projection onto the tangent cone:

$$\left\{ \begin{array}{l} x(t) \in C_0, \\ \dot{x}(t) = \mathrm{Proj}(\mathcal{T}(C_0, x(t)); S\tau(x(t))), \\ x(0) = x_0. \end{array} \right. \qquad (4)$$

<footnote n="1">¹For every vector $a \in \mathbb{R}^n$ and $b \in \mathbb{R}^n$, $\mathrm{diag}(a)b = \mathrm{diag}(b)a = (a_i b_i)_{i \in \{1, \dots, n\}}$.</footnote>
This model (4) has a unique solution [25, 5]. Since our system (3) is different with a projection onto the intersection of $\mathcal{T}(C_0, x(t))$ and $D(x(t))$, we must prove the existence of a solution, as defined below.

**Definition 1.1.** Let $x_0 \in \mathbb{R}_+^n$. A function $x : [0, t_f] \to C_0$ is a solution of (3) if $x$ is absolutely continuous and satisfies (3) for almost every $t \in [0, t_f]$.

Before proving existence, we compare our kinetic model with others from the literature.

# 2 Link with other models using reference examples

In this section, we use two reference examples, where minerals can react with aqueous or mineral species [13, 14]. Using the first example $W - M$, we prove that our projected dynamical model is equivalent in this case to complementary and discontinuous models. With the second example, we compare our model to a differential inclusion approach and illustrate the effect of constraints on time derivatives.

## 2.1 Example $W - M$

The first example is a single reaction with an aqueous species $W$ and a mineral species $M$ [14]:

$$W \Rightarrow M,$$

so that $S = (-1 \ 1)^T$ and $Q = (1 \ 1)^T$. The reaction rate is given by $\tau(x(t)) = k_r x_1(t) - k_p$, with two constants $k_r$ and $k_p$.

In [16, 17], the authors propose a complementarity formulation of kinetic systems, assuming one mineral per reaction and one reaction per mineral:

$$ \begin{cases} \dot{x}_1(t) + \dot{x}_2(t) = 0, \\ 0 \le x_2(t) \perp \dot{x}_2(t) - \tau(x(t)) \ge 0, \\ x(0) = x_0. \end{cases} \qquad (5) $$

We prove the following result.

**Proposition 2.1.** For the reference example $W - M$, the projected dynamical system (3) is pointwise equivalent to the complementarity system (5). Moreover, it is also equivalent to discontinuous models [16].

*Proof.* Let $x : [0, t_f] \to \mathbb{R}_+^2$ be a Lipschitz continuous function. We can split $[0, t_f]$ into a finite number of subintervals such that, in each open interval $I$, $x$ is differentiable and satisfies one of the three following conditions: first $\forall t \in I, x_1(t) = 0$, second $\forall t \in I, x_2(t) = 0$, third $\forall t \in I, x_1(t) > 0$ and $x_2(t) > 0$.
In the first case, since $x_1(t) = 0$ over $I$, it follows that $\forall t \in I, \dot{x}_1(t) = 0$. Hence, using the constraint $\dot{x}_1(t) + \dot{x}_2(t) = 0$, we get $\forall t \in I, \dot{x}_1(t) = \dot{x}_2(t) = 0$. Therefore, both models are equivalent in this case. The second case is similar using $x_2$ instead of $x_1$.

In the third case, the complementarity problem (5) can be rewritten

$$ \forall t \in I, \dot{x}(t) = (-\tau(x(t)), \tau(x(t)))^T = S\tau(x(t)). $$

Moreover, in this case,

$$ \mathrm{Proj}(D(x(t)) \cap \mathcal{T}(C_0; x(t)); S\tau(x(t))) = S\tau(x(t)), $$

so that the projected dynamical system is equivalently rewritten

$$ \forall t \in I, \dot{x}(t) = S\tau(x(t)). $$

Finally, both models are equivalent in each open interval $I$, thus almost everywhere in $[0, t_f]$.

Discontinuous or differential inclusion models [3, 14] are equivalent to the complementarity system [16] for this example, so that our model is also equivalent to them. □

## 2.2 Example $M_1 - M_2 - M_3$

The second example is a kinetic system with two reactions involving three minerals, where the second mineral $M_2$ participates in the two reactions [14]:

$$ M_1 \Leftrightarrow M_2 \text{ and } M_2 \Leftrightarrow M_3. $$

The stoichiometry matrix is $S = \begin{pmatrix} -1 & 1 & 0 \\ 0 & -1 & 1 \end{pmatrix}^T$ and a basis of $\mathrm{Ker}(S^T)$ is $Q = (1 \ 1 \ 1)^T$.

The reaction rate is constant and given by $\tau(x) = k_r - k_p$.

In this example, the differential inclusion model proposed in [14] is written

$$ \left\{ \begin{array}{l} x(t) \in C_0, \\ \dot{x}(t) = Sr(x(t)), \\ x(0) = x_0. \end{array} \right. \qquad (6) $$

We consider the numerical case where the initial condition is $x_0 = (1 \ 1 \ 1)^T$ and the reaction rate is $\tau = (2 \ 6)^T$ so that $S\tau = (-2 \ -4 \ 6)^T$. In this case, discontinuous reaction rates $r$ are defined in [14] by

$$ \left\{ \begin{array}{l} r_1(x(t)) = H(x_1(t))\tau_1, \\ r_2(x(t)) = H(x_2(t))\tau_2 + \delta(x_2(t))H(x_1(t))\tau_1, \end{array} \right. \qquad (7) $$
where $H$ is the heaviside function defined for $x \ge 0$ by $H(x) = 1$ if $x > 0$ and $H(0) = 0$, $\delta$ is the Dirac function defined for $x \ge 0$ by $\delta(0) = 1$ and $\delta(x) = 0$ if $x > 0$. The piecewise affine solution is plotted in Figure 1 and defined below, with $t_f = 1$:

$$x_1(t) = \begin{cases} 1 - 2t, & 0 \le t \le \frac{1}{2} \\ 0, & \frac{1}{2} \le t \le t_f, \end{cases}$$

$$x_2(t) = \begin{cases} 1 - 4t, & 0 \le t \le \frac{1}{4}, \\ 0, & \frac{1}{4} \le t \le t_f, \end{cases}$$

$$x_3(t) = \begin{cases} 1 + 6t, & 0 \le t \le \frac{1}{4}, \\ 2 + 2t, & \frac{1}{4} \le t \le \frac{1}{2}, \\ 3, & \frac{1}{2} \le t \le t_f. \end{cases}$$

<figure><img src="image_3.png" /><figcaption>Figure 1: Solution of the differential inclusion model (6) (7) on the reference example $M_1 - M_2 - M_3$ with $\tau = (2 \ 6)^T$ and $x_0 = (1 \ 1 \ 1)^T$.</figcaption></figure>

At time $t = \frac{1}{4}$, the second mineral $M_2$ is fully dissolved and at time $t = \frac{1}{2}$, the first mineral $M_1$ is also fully dissolved. Afterwards, only the third mineral $M_3$ remains and is constant.

The solution of the classical projected dynamical model (4), without the third constraint of our model (3), is now plotted in Figure 2.
<figure><img src="image_4.png" /><figcaption>Figure 2: Solution of the classical projected dynamical system (4) on the reference example $M_1 - M_2 - M_3$ with $\tau = (2\ 6)^T$ and $x_0 = (1\ 1\ 1)^T$.</figcaption></figure>

This solution is different from the previous one. At time $t = \frac{1}{4}$, the second mineral $M_2$ is fully dissolved, as in the differential inclusion model. But the first mineral $M_1$ is fully dissolved before time $t = 0.4$. Indeed, the solution of (4) satisfies $\dot{x}_1(t) = S\tau_1\lambda_1$, with $\lambda_1 > 1$, whereas the solution of (6) (7) satisfies $\dot{x}_1(t) = S\tau_1$. We observe that, for the solution of (4), the third constraint is not satisfied since $\dot{x}(t) \notin D(x(t))$. This example illustrates why we introduce constraints on time derivatives. Indeed, we prove the following result.

**Proposition 2.2.** The solution of the differential inclusion model (6) (7) applied to example $M1 - M2 - M3$ with $x_0 = (1\ 1\ 1)^T$ and $\tau = (2\ 6)^T$ is also solution of our specific projected dynamical model (3).

*Proof*. The solution $x$ of the differential inclusion model is absolutely continuous and differentiable for all $t \in [0, +\infty[ \setminus (\{\frac{1}{4}\} \cup \{\frac{1}{2}\})$.

The first conditions are easy to check: $\forall t, x(t) \in C_0$ and $\dot{x}(t) \in D(x(t)) \cap \mathcal{T}(C_0; x(t))$.

It remains to verify that $\dot{x}(t)$ is indeed the projection $y(t)$ of $S\tau$ onto the set $D(x(t)) \cap \mathcal{T}(C_0; x(t))$ for all $t$. It is true for $t \in [0, \frac{1}{4}[$ since $\dot{x}(t) = S\tau$. For $t \in ]\frac{1}{2}, 1]$), the solution $x$ satisfies $x_1 = x_2 = 0$ and $(S\tau)_1 < 0, (S\tau)_2 < 0$, thus $y(t) = (0, 0, 0)^T = \dot{x}(t)$. For $t \in ]\frac{1}{4}, \frac{1}{2}[$, $x_2 = 0$ and $(S\tau)_2 < 0$, thus $y_2(t) = 0$. Moreover, $Q^T y(t) = y_1(t) + y_3(t) = 0$ thus $y(t) = \text{diag}(S\tau)\mu = (-2\mu, 0, 2\mu)^T$ with $0 \le \mu \le 1$. Finally, $\mu = 1$ for the projection of $S\tau$ so that $y(t) = \dot{x}(t)$.

<figure><img src="image_5.png" /></figure>

# 3 Existence of a Solution

In this section, we prove the existence of a solution to (3). We give a constructive proof by introducing a numerical scheme. For a given time step $h > 0$, we consider a forward difference approximation of the derivative. Then, an approximate solution is defined by
a projection and an affine interpolation. We prove that this explicit scheme converges as $h \to 0^+$ to a solution of (3).

Throughout this section, we denote by $[a]$ the largest integer smaller than $a$.

## 3.1 Explicit scheme and interpolation

Consider a fixed step size $h > 0$ and an initial condition $x_0 \ge 0$. For any $x \in C_0$ we define a set $K$ by

$$K = C_0 \cap (x + hD(x)). \qquad (8)$$

The set $K$ is non empty, closed and convex. Therefore, the projection of a vector onto the set $K$ exists and is unique. For $k \in \mathbb{N}$ such that $kh \le t_f$, let $\{x_k\}_k$ be the sequence defined by

$$\begin{cases} x_0 \in C_0, \\ K_k = C_0 \cap (x_k + hD(x_k)), \\ x_{k+1} = \operatorname{Proj}(K_k; x_k + hS\tau(x_k)). \end{cases} \qquad (9)$$

This sequence is well-defined since $x_k \in C_0$. We now sum up in the following lemma some properties of the sequence $\{x_k\}_k$ defined by (9).

**Lemma 3.1.** Let $h > 0$, $k \in \mathbb{N}$ such that $hk \le t_f$, and $x_k \in C_0$.

Then

$$z \in K_k \iff \frac{z - x_k}{h} \in D(x_k) \cap \{d : Q^T d = 0, d \ge -\frac{x_k}{h}\}.$$

Moreover, (9) is equivalent to

$$\left\{ \begin{array}{l} x_{k+1} = x_k + hd_k, \\ \frac{x_{k+1} - x_k}{h} = \operatorname{Proj} (D(x_k) \cap \{d : Q^T d = 0, d \ge -\frac{x_k}{h}\} ; S\tau(x_k)). \end{array} \right. \qquad (10)$$

*Proof.* Let $z$ and $d$ be related by $z = x_k + hd$. Since $x_k \in C_0$, it holds that

$$z \in C_0 \iff Q^T d = 0, d \ge -\frac{x_k}{h}$$

and

$$z \in x_k + hD(x_k) \iff d \in D(x_k).$$

Thus $z \in K_k$ if and only if $d \in D(x_k) \cap \{d : Q^T d = 0, d \ge -\frac{x_k}{h}\}$. Also, the projection in (9) can be rewritten as

$$\|x_{k+1} - x_k - hS\tau(x_k)\| \le \|z - x_k - hS\tau(x_k)\|, \quad \forall z \in K_k.$$

Hence, dividing by $h > 0$ and defining $d_k$ by $x_{k+1} = x_k + hd_k$, previous inequality is equivalent to

$$\|d_k - S\tau(x_k)\| \le \|d - S\tau(x_k)\|, \quad \forall d \in D(x_k) \cap \{d : Q^T d = 0, d \ge -\frac{x_k}{h}\}.$$

The result follows again by definition of a projection. $\square$
Using the sequence {$x_k$}$_k$, for $h > 0$ and $x_0 \in C_0$, we define an interpolation function $x_h : [0, t_f] \to \mathbb{R}^n$ by

$$x_h(t) = x_k + (x_{k+1} - x_k) \frac{t - kh}{h}, \quad \text{with } k = \lfloor \frac{t}{h} \rfloor. \qquad (11)$$

Here, $k$ is chosen such that $t \in [kh, (k+1)h[$.

The function $x_h$ is piecewise affine and differentiable almost everywhere:

$$\dot{x}_h(t) = \frac{x_{k+1} - x_k}{h}, \quad \text{with } k = \lfloor \frac{t}{h} \rfloor, \text{ and } t \neq lh, l = 0, \dots, \lfloor \frac{t_f}{h} \rfloor. \qquad (12)$$

The function $x_h$ satisfies an approximate model:

**Lemma 3.2.** Let $h > 0$ and $x_0 \in C_0$. The function $x_h$ is solution of the following approximate projected dynamical system: for all $t \neq lh, l = 0, \dots, \lfloor \frac{t_f}{h} \rfloor$,

$$\left\{ \begin{array}{l} \dot{x}_h(t) = \operatorname{Proj}(D(x_k) \cap \{d: Q^T d = 0, d \geq -\frac{x_k}{h}\}; S\tau(x_k)), \\ x_h(t) \in C_0, \\ x_h(0) = x_0, \end{array} \right. \qquad (13)$$

where $k = \lfloor \frac{t}{h} \rfloor$.

*Proof.* For any $t \in [0, t_f]$, $x_h(t)$ is a convex combination of $x_k \in C_0$ and $x_{k+1} \in C_0$. Since $C_0$ is a convex set, $x_h(t) \in C_0$. By definition, $x_h(0) = x_0$.

Finally, by (12) and (10), the derivative satisfies the first equation of (13). $\square$

## 3.2 Bounded constraints

Kinetic systems have usually several invariants directly related to conservative variables $Q$. We assume the existence of a particular one [17, 14].

**Assumption 3.1.**

$$\{q \in \operatorname{Ker}(S^T) \mid q > 0\} \neq \emptyset. \qquad (14)$$

In particular, let $q$ be the vector of molar mass of species. It is natural to assume that $S^T q = 0$ and $q > 0$.

Using Motzkin alternative theorem [24] or Gordan's theorem [12], Assumption 3.1 is equivalent to

$$Q^T d = 0, d \geq 0 \Rightarrow d = 0.$$

Thanks to Assumption 3.1, the concentrations are necessarily finite.

**Lemma 3.3.** *Assuming 3.1, the set $C_0$ is bounded.*

The proof can be found in [14].

A direct consequence of this lemma is that the set $D(x)$ is uniformly bounded.
**Lemma 3.4.** *Assuming 3.1, there exists $\rho > 0$ such that, for any $x \in C_0$ and any $d \in D(x)$, $\|d\| \le \rho$.*

*Proof.* Let $x \in C_0$ and $d \in D(x)$. By definition, $d = \text{diag}(S\tau(x))\lambda$ for some $\lambda \in [0, 1]^n$. Therefore

$$ \|d\| \le \|\lambda\| \|S\tau(x)\| \le \|S\tau(x)\|. $$

Since $C_0$ is bounded by Lemma 3.3, and $\tau$ is continuous, then $\tau(x)$ is uniformly bounded. Hence, there exists $\rho > 0$ such that $\|d\| \le \rho$. $\square$

### 3.3 Convergence of the function $x_h$

Throughout this section, we make Assumption (3.1). We first prove convergence of $x_h$ towards a function $x$.

**Proposition 3.1.** *There exists a subsequence $\mathcal{H}$ and a Lipschitz continuous function $x : [0, t_f] \to C_0 \subset \mathbb{R}^n$ such that $x_h \to_{\mathcal{H}} x$ uniformly.*

*Proof.* By Lemma 3.3, the sequence $\{x_h\}_h$ is equi-bounded. Also, by Lemma 3.4, the sequence $\{x_h\}_h$ is equi-Lipschitz with constant $\rho$.

Then, all the assumptions of the Arzela-Ascoli theorem are satisfied. Therefore, there exists a subsequence $\mathcal{H}$, such that $h \to_{\mathcal{H}} 0$, and a Lipschitz mapping $x : [0, t_f] \to \mathbb{R}^n$ such that $x(0) = x_0$ and $\{x_h\}_{h \in \mathcal{H}}$ converges uniformly to $x$ in $[0, t_f]$. Besides, since $C_0$ is closed, then $x(t) \in C_0$ for all $t \in [0, t_f]$. $\square$

Since $k$ depends on $h$ and $t$, the sequence $x_k$ converges.

**Lemma 3.5.** Let $t \in [0, t_f]$. The sequence $x_k$ with $k = \lfloor \frac{t}{h} \rfloor$ converges towards $x(t): x_k \to_{\mathcal{H}} x(t)$.

*Proof.* By definition of $x_h$ in (11) and using Lemma 3.4, it follows

$$ \|x_h(t) - x_k\| \le \|(x_{k+1} - x_k)\frac{t - kh}{h}\| \le \rho(t - kh) \le \rho h. $$

Passing to the limit as $h \to 0^+$, and using convergence of $x_h(t)$ in Proposition 3.1, we get convergence of $x_k$ towards $x(t)$. $\square$

Another consequence of Proposition 3.1 is the weak convergence of the derivative.

**Lemma 3.6.** $\dot{x}_h \to_{\mathcal{H}} \dot{x}$ weakly in $L^1([0, t_f], \mathbb{R}^n)$.

*Proof.* The result is a straightforward application of [18, Lemma 2.2]. $\square$

Even though the convergence is weak, there exists a sequence of functions which converge strongly towards the derivative of $x$.
**Lemma 3.7.** There exists a sequence of functions {$y_h$}_h → H and satisfying

$$y_h = \sum_{l \in J_h} \alpha_l \dot{x}_{hl}, \text{ with } \sum_{l \in J_h} \alpha_l = 1, \text{ and } \alpha_l \ge 0, \forall l \in J_h,$$

where $J_h$ is a finite set of indices such that $(h_l)_{l \in J_h} \subset H$.

*Proof.* The sequence {$\dot{x}_h$}_{h∈H} converges weakly to $\dot{x}$, by Lemma 3.6. Then, applying Mazur's Lemma, there exists a sequence {$y_h$}_h, with values in the convex hull of the set {$\dot{x}_h$}_{h∈H}, which converges to $\dot{x}$. □

The following lemma is useful when passing to the limit in the constraint set. The set convergence has to be understood in the Painlevé-Kuratowski sense [26, Chapter 5.B].

**Lemma 3.8.** Let $z \in C_0$ and {$z_h$}_h ⊂ C_0 be a sequence such that $z_h \to_{h \to 0^+} z$. Then

$$\lim_{h \to 0^+} D(z_h) \cap \{d : Q^T d = 0, d \ge -\frac{z_h}{h}\} = D(z) \cap \{d : Q^T d = 0, d \ge -\lim_{h \to 0^+} \frac{z_h}{h}\}.$$

*Proof.* By Lemma 3.3 and Lemma 3.4, the sets $C_0$ and $D(z_h) \cap \{d : Q^T d = 0, d \ge -\frac{z_h}{h}\}$ are both non-empty and compact for all $h > 0$. Therefore, by [26, Theorem 5.7], it follows that $\limsup_{h \to 0^+} D(z_h) \cap \{d : Q^T d = 0, d \ge -\frac{z_h}{h}\} \subseteq D(z) \cap \{d : Q^T d = 0, d \ge -\lim_{h \to 0^+} \frac{z_h}{h}\}$.

On the other hand, let $\xi \in D(z) \cap \{d : Q^T d = 0, d \ge -\lim_{h \to 0^+} \frac{z_h}{h}\}$. Since $x \mapsto S\tau(x)$ is continuous, there exists a sequence $\xi_h \in D(z_h) \cap \{d : Q^T d = 0, d \ge -\frac{z_h}{h}\}$ which converges to $\xi$. Thus $D(z) \cap \{d : Q^T d = 0, d \ge -\lim_{h \to 0^+} \frac{z_h}{h}\} \subseteq \liminf_{h \to 0^+} D(z_h) \cap \{d : Q^T d = 0, d \ge -\frac{z_h}{h}\}$.

All in all, the limit exists and the result follows. □

## 3.4 Properties of the limit function

Again, we make Assumption 3.1 throughout this section.

In order to prove that the limit $x$ of $x_h$ is a solution of the model, we start by showing that it satisfies the constraints.

We define the set $I$ by $I = \{t \in [0, t_f] : \dot{x}(t)$ exists\}$. By Sard's theorem, the set $[0, t_f] \setminus I$ has a zero measure.

**Lemma 3.9.** For all $t \in I$, $\dot{x}(t) \in D(x(t)) \cap \{d : Q^T d = 0, d \ge -\lim_{h \to 0} \frac{x_{kl}}{h}\}$, with $l \in J_h$ and $k_l = \lfloor \frac{t}{h_l} \rfloor$.

*Proof.* Let $t \in I$ and $h > 0$. Considering $y_h$ defined in Lemma 3.7, and using (13), we get

$$y_h(t) \in \operatorname{conv} \left( \bigcup_{l \in J_h} D(x_{kl}) \cap \{d : Q^T d = 0, d \ge -\frac{x_{kl}}{h_l}\} \right),$$

where $\operatorname{conv}$ denotes the convex hull and $k_l = \lfloor \frac{t}{h_l} \rfloor$.

We note that $h_l \to_H 0$ and, using Lemma 3.5, $x_{kl} \to_H x(t)$. Thus the sequence $x_{kl}(t)$ satisfies assumptions of Lemma 3.8. Passing to the limit in $y_h$, we get the result. □
**Lemma 3.10.** For all $t \in I$, $\dot{x}(t) \in \mathcal{T}(C_0; x(t))$.

*Proof.* Recall that the tangent cone to $C_0$ is given by (2). By Lemma 3.9, $Q^T \dot{x}(t) = 0$ and $\dot{x}_i(t) = (S\tau(x(t)))_i \lambda_i$ for $i = 1, \dots, n$, with $0 \le \lambda_i \le 1$. Let $i$ such that $x_i(t) = 0$. We must show that $\dot{x}_i(t) \ge 0$.

If $(S\tau(x(t)))_i \ge 0$, then $\dot{x}_i(t) \ge 0$. On the other hand, if $(S\tau(x(t)))_i < 0$, by continuity of $\tau$ and $x$, there is an open interval $U$ centered in $t$ such that $(S\tau(x(s)))_i < 0$ so that $\dot{x}_i(s) \le 0$ for all $s \in U \cap I$. Also, $x_i(t) = 0$ and $x_i(s) \ge 0$. Since $x$ is absolutely continuous, we get $\dot{x}_i(t) = 0$. This concludes the proof. $\square$

**Lemma 3.11.** The limit $x$ satisfies the inclusion below, for $t \in I$ and $l \in J_h$:

$$D(x(t)) \cap \mathcal{T}(C_0; x(t)) \subseteq D(x(t)) \cap \{d : Q^T d = 0, d \ge -\lim_{h \in \mathcal{H}} \frac{x_{kl}}{h_l}\}.$$

*Proof.* Let $L = \lim_{h \in \mathcal{H}} \frac{x_{kl}}{h_l}$. Since $x_{kl} \ge 0$ the limit satisfies $L \ge 0$. Let $d \in \mathcal{T}(C_0; x(t))$. If $x_i = 0$, then $d_i \ge 0 \ge -L_i$. If $x_i > 0$, then $L = +\infty$ and $d_i \ge -L_i$. $\square$

Finally, we prove the main theorem of existence.

**Theorem 3.1.** The limit function $x$ is a solution of (3).

*Proof.* The function $x : [0, t_f] \to C_0$ is a Lipschitz continuous function such that $x(0) = x_0$. It satisfies the constraints of the model, so that it remains to verify the projection property.

Considering $y_h$ defined in Lemma 3.7, and using the triangle inequality, we get

$$\|S\tau(x_k) - y_h(t)\| \le \max_{l \in J_h} \|S\tau(x_k) - \dot{x}_{hl}(t)\|. \quad (15)$$

Let $j \in J_h$ which maximizes this inequality and $d \in D(x_{kj}) \cap \{d : Q^T d = 0, d \ge -\frac{x_{kj}}{h_j}\}$. Since $x_{h_j}$ is solution of (13), it satisfies

$$\|S\tau(x_{kj}) - \dot{x}_{h_j}(t)\| \le \|S\tau(x_{kj}) - d\|,$$

therefore

$$\|S\tau(x_k) - y_h(t)\| \le \|S\tau(x_k) - S\tau(x_{kj})\| + \|S\tau(x_{kj}) - d\|.$$

Thanks to Lemma 3.5, $\lim_{h \in \mathcal{H}} x_{kj} = \lim_{h \in \mathcal{H}} x_k = x(t)$. Using Lemma 3.7 and Lemma 3.8, we can pass to the limit:

$$\|S\tau(x(t)) - \dot{x}(t)\| \le \|S\tau(x(t)) - d\|, d \in D(x(t)) \cap \{d : Q^T d = 0, d \ge -\lim_{h \in \mathcal{H}} \frac{x_{kj}}{h_j}\}.$$

Finally, Lemma 3.11 provides the result. $\square$
# 4 Numerical Experiments

## 4.1 An Alternating Projection Algorithm

We propose an alternating projection method to compute the sequence $x_k$. This iterative algorithm consists in computing a sequence $\{y_j\}_j$, starting from $x_k + hS\tau(x_k)$, which converges towards $x_{k+1}$. The sequence is defined as follows:

$$ 
\begin{align} 
y_0 &= x_k + hS\tau(x_k), e_0 = e_{-1} = 0, \nonumber \\ 
y_{2j+1} &= \operatorname{Proj}(C_0; y_{2j} + e_{2j-1}), \nonumber \\ 
e_{2j+1} &= y_{2j} + e_{2j-1} - y_{2j+1}, \tag{16} \\ 
y_{2j+2} &= x_k + h\operatorname{Proj}\left(D(x_k); \frac{y_{2j+1} + e_{2j} - x_k}{h}\right), \nonumber \\ 
e_{2j+2} &= y_{2j+1} + e_{2j} - y_{2j+2}. \nonumber 
\end{align} 
$$

We prove the following result.

**Theorem 4.1.** The sequence $\{y_j\}_j$ computed by (16) converges to $x_{k+1}$ satisfying (9).

*Proof.* Using [8, Proposition 2.7 (ii) and (iv)], the fourth line of (16) can be rewritten

$$ y_{2j+2} = \operatorname{Proj}(x_k + hD(x_k); y_{2j+1} + e_{2j}). $$

The convergence of the sequence is a direct application of the Boyle-Dykstra Theorem [8, Theorem 9.24]. We already noted that $C_0$ and $x_k + hD(x_k)$ are closed convex set and that their intersection $K_k$ is non-empty. Therefore, the sequence $y_j$ converges and the limit is nothing else than $x_{k+1}$. □

We develop a software in Matlab in order to compute the sequence given by (16). For each $k = 0, \dots, \lfloor \frac{tf}{h} \rfloor$, and given $x_k$, we implement the alternating projection algorithm 1 to compute the new iterate $x_{k+1}$.

**Data:** starting vector $x_k$; iteration threshold $\epsilon$;

1. Begin;

2. Set $j := 0, y_0 := x_k + hS\tau(x_k), e_0 := e_{-1} := 0$;

3. repeat

4. $y_{2j+1} = \operatorname{Proj}(C_0; y_{2j} + e_{2j-1})$;

5. $e_{2j+1} = y_{2j} + e_{2j-1} - y_{2j+1}$;

6. $d_j := \operatorname{Proj}\left(D(x_k); \frac{y_{2j+1} + e_{2j} - x_k}{h}\right)$;

7. $y_{2j+2} := x_k + h d_j$;

8. $e_{2j+2} = y_{2j+1} + e_{2j} - y_{2j+2}$;

9. $j := j + 1$;

10. until $\|y_{2j} - y_{2j-1}\| \le \epsilon$;

11. return: $x_{k+1} := y_{2j}$.

**Algorithm 1:** Alternating projection algorithm.
Projections are computed approximately using Matlab lsqlin function with defaults parameters. In all experiments, the threshold $\epsilon$ is set to $10^{-7}$ and convergence is obtained when the step size $h$ tends towards 0. We show the results with a step size $h$ equal to $10^{-2}$.

We compute numerical solutions of the two reference examples defined in Section 2, and use also a third example, which is a variant of the second one.

## 4.2 Example $W - M$

<figure><img src="image_6.png" /><figcaption>Figure 3: Example $W - M$, with $k_r = 1$, $k_p = 2$ and $x_0 = (1 \ 0.5)^T$.</figcaption></figure>

The first numerical experiment deals with the reference example $W - M$. Figure 3 illustrates dissolution of the mineral $M$. As expected, the chemical reaction stops when $M$ is fully dissolved. Asymptotic equilibrium is given by $x_2 = 0$ and $Q^T x = Q^T x_0$.

## 4.3 Example $M_1 - M_2 - M_3$

The second numerical experiment concerns the reference example $M_1 - M_2 - M_3$. Figure 4 shows the evolution of the three minerals. The mineral $M_2$, which participates into both reactions, vanishes first. Then, the mineral $M_1$ decreases until full dissolution, until only the mineral $M_3$ remains and the reaction stops. Asymptotic equilibrium is given by $x_1 = x_2 = 0$ and $Q^T x = Q^T x_0$.

## 4.4 Example $W - M_2 - M_3$

This example is a variant where mineral $M_1$ is replaced by an aqueous species $W$:

$$W \rightleftharpoons M_2 \text{ and } M_2 \rightleftharpoons M_3.$$
<figure><img src="image_7.png" /><figcaption>Figure 4: Example $M_1 - M_2 - M_3$, with $k_r = (3\ 7)^T$, $k_p = (1\ 1)^T$ and $x_0 = (1.5\ 1\ 0)^T$.</figcaption></figure>

The stoichiometric matrix is the same, as well as the second reaction rate, but the first reaction rate is now

$$ \tau_1(x) = k_{r,1}x_1 - k_{p,1}. $$

<figure><img src="image_8.png" /><figcaption>Figure 5: Example $W - M_1 - M_2$, with $k_r = (3\ 7)^T$, $k_p = (1\ 1)^T$ and $x_0 = (1.5\ 1\ 0)^T$.</figcaption></figure>

Figure 5 shows the dynamics of the chemical system. The aqueous species W precipitates and produces the mineral $M_2$ which is dissolved into the other mineral $M_3$. Under the conditions chosen, $M_2$ is fully dissolved in a finite time. However, because the first reaction continues, the second one can also continue, at the reaction rate $\tau_1$ instead of $\tau_2$. In some sense, the system behaves as if W precipitates into mineral $M_3$, but at rate $\tau_1$. Asymptotic equilibrium, here at infinite time, is given by $x_2 = 0$, $Q^T x = Q^T x_0$ and $\tau_1(x_1) = 0$.
# 5 Concluding Discussion

In this paper, we consider systems of kinetically-controlled precipitation-dissolution geochemical reactions where minerals can react in several reactions and with other minerals. This extends the study in [17] where a complementarity model is applied to systems with one mineral in one reaction. It also provides a model and a numerical method different from the differential inclusion model with a regularization approach [14]. Our new model defines a projected dynamical system with additional constraints on time derivatives. In order to prove existence of a solution, we develop a numerical scheme, where the approximate system has a unique solution. Then we prove convergence towards a solution of the projected dynamical model. This numerical scheme is implemented using an alternating projection algorithm, yielding easy computations of approximate solutions. We also consider two reference examples, which serve to compare our projection model to complementarity and discontinuous models. Finally, these examples validate computations with our numerical algorithm.

Further research may focus on the numerical perspectives offered by this work. Our scheme is based on an explicit discretization. The perspective to use an implicit scheme is an interesting path but satisfying the constraints would be more difficult. Kinetic reactions are usually too slow to get equilibrium. However, when an asymptotic state is reached, it would be interesting to make a connection with models for equilibrium reactions [9]. Finally, an essential motivation in this work is to couple this model in reactive transport applications.

# References

[1] A. AGOSTI, L. FORMAGGIA, AND A. SCOTTI, Analysis of a model for precipitation and dissolution coupled with a Darcy flux, Journal of Mathematical Analysis and Applications, 431 (2015), pp. 752-781.

[2] L. AMIR AND M. KERN, A global method for coupling transport with chemistry in heterogeneous porous media, Computational Geosciences, 14 (2010), pp. 465-481.

[3] N. BOUILLARD, R. EYMARD, R. HERBIN, AND P. MONTARNAL, Diffusion with dissolution and precipitation in a porous medium: Mathematical analysis and numerical approximation of a simplified model, ESAIM: Mathematical Modelling and Numerical Analysis, 41 (2007), pp. 975-1000.

[4] J. CARRAYROU, J. HOFFMANN, P. KNABNER, S. KRÄUTLE, C. DE DIEULEVEULT, J. ERHEL, J. VAN DER LEE, V. LAGNEAU, K. U. MAYER, AND K. T. B. MACQUARRIE, Comparison of numerical methods for simulating strongly nonlinear and heterogeneous reactive transport problems—the MoMas benchmark case, Computational Geosciences, 14 (2010), pp. 483-502.
[5] M. G. COJOCARU, P. DANIELE, AND A. NAGURNEY, *Projected dynamical systems and evolutionary variational inequalities via Hilbert spaces with applications*, Journal of Optimization Theory and Applications, 127 (2005), pp. 549–563.

[6] C. DE DIEULEVEULT, *Un modèle numérique global et performant pour le couplage géochimie-transport*, PhD thesis, University of Rennes, 2008.

[7] C. DE DIEULEVEULT, J. ERHEL, AND M. KERN, *A global strategy for solving reactive transport equations*, Journal of Computational Physics, 228 (2009), pp. 6395–6410.

[8] F. R. DEUTSCH, *Best approximation in inner product spaces*, Springer Science & Business Media, 2012.

[9] J. ERHEL AND T. MIGOT, *Characterizations of solutions in geochemistry: existence, uniqueness, and precipitation diagram*, Computational Geosciences, 23 (2019), pp. 523–535.

[10] J. ERHEL AND S. SABIT, *Analysis of a global reactive transport model and results for the MoMas benchmark*, Mathematics and Computers in Simulation, 137 (2017), pp. 286–298.

[11] J. C. FRIEDLY AND J. RUBIN, *Solute transport with multiple equilibrium-controlled or kinetically-controlled chemical reactions*, Water resources research, 28 (1992), pp. 1935–1953.

[12] P. GORDAN, *Über die auflösung linearer gleichungen mit reellen coefficienten*, Mathematische Annalen, 6 (1873), pp. 23–28.

[13] B. HAMLAT, *Mathematical modeling of multiphase kinetic reactions in geochemistry*, PhD thesis, University of Rennes, 2020.

[14] B. HAMLAT, J. ERHEL, A. MICHEL, AND T. FANEY, *Analysis and numerical computation of geochemical systems with kinetic precipitation and dissolution reactions involving several minerals*, HAL, (2021).

[15] J. HOFFMANN, *Reactive transport and mineral dissolution/precipitation in porous media: efficient solution algorithms, benchmark computations and existence of global solutions*, PhD thesis, Friedrich-Alexander-Universität Erlangen-Nürnberg, 2010.

[16] J. HOFFMANN, S. KRÄUTLE, AND P. KNABNER, *A general reduction scheme for reactive transport in porous media*, Computational Geosciences, 16 (2012), pp. 1081–1099.

[17] J. HOFFMANN, S. KRÄUTLE, AND P. KNABNER, *Existence and uniqueness of a global solution for reactive transport with mineral precipitation-dissolution and aquatic reactions in porous media*, SIAM Journal on Mathematical Analysis, 49 (2017), pp. 4812–4837.
[18] A. JOURANI AND E. VILCHES, *Moreau-Yosida regularization of state-dependent sweeping processes with nonregular sets*, Journal of Optimization Theory and Applications, 173 (2017), pp. 91–116.

[19] P. KNABNER, C. J. VAN DUIJN, AND S. HENGST, *An analysis of crystal dissolution fronts in flows through porous media. Part 1: Compatible boundary conditions*, Advances in Water Resources, 18 (1995), pp. 171–185.

[20] S. KRÄUTLE, *General multi-species reactive transport problems in porous media: Efficient numerical approaches and existence of global solutions*, habilitation thesis, University of Erlangen-Nuremberg, Department of Mathematics, Erlangen, 2008.

[21] K. KUMAR, I. S. POP, AND F. A. RADU, *Convergence analysis for a conformal discretization of a model for precipitation and dissolution in porous media*, Numerische Mathematik, 127 (2014), pp. 715–749.

[22] A. M. LEAL, D. A. KULIK, W. R. SMITH, AND M. O. SAAR, *An overview of computational methods for chemical equilibrium and kinetic calculations for geochemical and reactive transport modeling*, Pure and Applied Chemistry, 89 (2017), pp. 597–643.

[23] P. C. LICHTNER, C. I. STEEFEL, AND E. H. OELKERS, *Reactive transport in porous media*, Walter de Gruyter GmbH & Co KG, 2019.

[24] O. L. MANGASARIAN, *Nonlinear programming*, Society for Industrial and Applied Mathematics, 1994.

[25] A. NAGURNEY AND D. ZHANG, *Projected dynamical systems and variational inequalities with applications*, Springer Science & Business Media, 1996.

[26] R. T. ROCKAFELLAR AND R. J.-B. WETS, *Variational analysis*, vol. 317, Springer Science & Business Media, 2009.

[27] J. RUBIN, *Transport of reacting solutes in porous media: Relation between mathematical nature of problem formulation and chemical nature of reactions*, Water Resources Research, 19 (1983), pp. 1231–1252.

[28] S. SABIT, *Numerical methods for reactive transport*, PhD thesis, University of Rennes, 2014.

[29] C. J. VAN DUIJN AND P. KNABNER, *Travelling wave behaviour of crystal dissolution in porous media flow*, European Journal of Applied Mathematics, 8 (1997), pp. 49–72.

[30] C. J. VAN DUIJN, P. KNABNER, AND R. J. SCHOTTING, *An analysis of crystal dissolution fronts in flows through porous media part 2: incompatible boundary conditions*, Advances in water resources, 22 (1998), pp. 1–16.
[31] C. J. VAN DUIJN AND I. S. POP, *Crystal dissolution and precipitation in porous media: Pore scale analysis*, Journal für die reine und angewandte Mathematik, 577 (2004), pp. 171–211.