Robust simulation of mineral precipitation–dissolution problems with variable mineral surface area

We consider a macroscale model of transport and reaction of chemical species in a porous medium with a special focus on mineral precipitation–dissolution processes. In the literature, it is frequently proposed that the reaction rate should depend on the reactive mineral surface area, and so on the amount of mineral. We point out that a frequently used model is ill posed in the sense that it admits non-unique solutions. We investigate what consequences this non-uniqueness has on the numerical solution of the model. The main novelty in this article is our proposal of a certain substitution which removes the ill-posedness from the system and which leads to better numerical results than some “ad hoc methods.” We think that the proposed substitution is a rather elegant way to get rid of the non-uniqueness and the numerical difficulties and is much less technical than other ideas. As a proof of concept, we present some numerical tests and simulations for the new model.

minerals and aqueous species. For many precipitation-dissolution processes, the reaction rate (per volume of the porous medium) should depend on the size of the reactive surface of the mineral within this volume (e.g., [2,5,[16][17][18][20][21][22][23][24]), where c is the concentration of a mineral (moles per total volume) and A is the reactive mineral surface per total volume. The factor r depends on the vector of aqueous species concentrations c aq ; we might consider the law of mass action, for example. For the reactive surface A, there are several models used and discussed in the literature. Some papers assume that the amount of precipitation-dissolution is independent of the mineral surface size, which can be modeled by setting A(c) = const in (1); see for example, [1,3,4,[9][10][11][12]14,15,29]. On the other hand, when assuming that the surface size determines the reaction velocity, then a frequently used classical model for a variable surface is with a constant A 0 > 0, for c ≥ 0. This concept is used by many authors, see [2,5,17,18,20,[22][23][24]. In particular, it is implemented in the simulation tools PFLOTRAN [19] and Crunch Flow [26] that are well known in the computational geoscientific community. The model (2) is justified by the following consideration on the microscale. If we assume that within a representative elementary volume (REV), the mineral consists of small equisized hemispheres, or spheres, or cubes of "radius" r , then the surface size of every mineral grain within the REV behaves like A grain ∼ r 2 and its volume or mass m grain like r 3 . Hence, A grain ∼ (m grain ) 2/3 , and therefore, reactive mineral surface and mineral amount should meet the relation (2) in every REV, and the mineral source term should contain a factor (2). Note that also the so-called shrinking particle/core models are based on the concept of these spherical grains.
There are other models possible, the assumption of platy mineral structures, growing only in two space dimensions (reactive surface A grain ∼r , mass m grain ∼r 2 ), may lead to (1) with A(c) = A 0 c 1/2 , instead of (2). There are also models focusing on large mineral amounts ("clogging" of the pore space) and describe a monotonically decreasing surface area-these scenarios are not considered in this article. There are models derived by upscaling relating the mineral surface to the porosity (i.e., to the mineral volume) [25,30], and the pore geometry is parametrized by the radius of the grains as a parameter (Sect. 4 in [30], Sect. 2.3 in [25]) to simplify the model-this radius corresponds to the variable ξ we are introducing in Sect. 6. There are also models assuming a population of crystals of different sizes within every REV or discretization point, instead of one radius r [7,16,27] and which have a closer look at the formation of mineral crystals (nucleation) and how the distribution of crystal sizes evolves in time (Ostwald ripening), see [27]. The price to pay for such a much more detailed model is its much larger computational amount of work. However, in this article, we are interested in models containing the term (2) and effects that arise for small mineral concentrations, caused by this term.
The ODE (1)-(2) will be coupled to a system of partial differential equations (PDEs) for the unknowns c aq (see Sect. 2), but for a moment-and for the sake of simplicity-let us just assume that c aq is given, so that (1)- (2) is an ODE c (t) = f (c(t)). One immediately realizes that (2), and therefore, the right-hand side of the ODE (1), is not locally Lipschitz continuous at c = 0, and therefore, the uniqueness of solutions is only clear as long as the solution stays away from the value c = 0. We may point out that a one-sided Lipschitz condition, i.e., if there is an L ∈ R such that would also be sufficient for the uniqueness of the solution of an initial value problem (see, for example, (12.2) and Theorem 12.1 in Chapter IV, in the textbook [8]). However, for our model of precipitation-dissolution reactions, even a one-sided Lipschitz condition does not hold. In fact, it would be fulfilled if r (c aq ) ≤ 0, i.e., if we would only consider scenarios in which the fluid is undersaturated with respect to the mineral reactions and so dissolution always exceeds precipitation; in an oversaturated setting with initial mineral concentration equal to zero, there are indeed infinitely many solutions (cf. Appendix A). Hence, first of all, we have to clarify which solution is the physically correct one. See Sect. 4 for a discussion. Then, an important issue is the question what impact the non-uniqueness of solutions has for numerical solvers. When there are multiple solutions, a numerical solver may follow the one or the other solution. The numerical scheme has to take care that the physical solution is computed. Possibilities are discussed in Sects. 5 and 6. Reactive transport is usually discretized using an implicit timestepping scheme, and the resulting nonlinear discrete problem has to be solved iteratively. The first method of choice is Newton's method, since it is (locally) quadratic convergent. We found that the term (2) causes difficulties for Newton's method to converge for mineral concentrations close to zero. The authors could not find any proposals in the literature how to deal with this specific difficulty, except the publication [24], where it is suggested to use a regularization of (2) or even a completely different type of law, such as a Monod-type law. Instead of modifying the law (2), another idea would be to hinder the mineral concentration to reach zero, i.e., to enforce c ≥ for an >0.
Both techniques, to replace or to regularize (2), or to force c ≥ , are approximations/modifications of the original model. Questions like how to choose or the regularization parameter in numerical simulations arise. However, we will show-and this is the main purpose of this publication-that no approximation/modification of the law (2) or of standard nonlinear solvers are necessary at all. Instead, the method we propose is to apply a simple substitution After the substitution is done, the model has a unique solution, and also the numerical difficulties vanish. The nicest thing about this substitution is that the idea is very simple, and also the implementation is easy. If at all, the implementation is even easier than for the original model. The paper is structured as follows. In Sect. 2, we present our multi-species multi-reaction model with mineral precipitation-dissolution reactions and with aqueous reactions. In Sect. 3, we present simplified versions of the model by omitting transport (hence, reducing the PDE-ODE model to an ODE model) and by limiting the complexity of the reactive network. The reason for this is the fact that the simpler model can be easier analyzed. In Sect. 4, we show that the model of Sect. 3 admits infinitely many solutions (the proof is given in Appendix A). In Sect. 5, we turn to the discretized version of the model and point out why we should expect convergence problems with Newton's method, caused by the term (2), for small mineral concentrations. We discuss "ad hoc methods," in particular, what we call the "epsilon method," i.e., the idea mentioned above that the solution should be forced to be larger or equal to an >0, to circumvent the difficulties, and we discuss how the numerical parameter should be chosen. Section 6 contains the most important novelty of this article. We present the substitution which lets the non-Lipschitz term vanish. In Appendix B, we prove that after the substitution, the model admits only one solution. In Sect. 7, we perform numerical tests for our "substitution method." We compare it in terms of accuracy, robustness, and number of Newton steps with the "epsilon method." We perform our numerical tests both for the ODE model from Sect. 3 and for the PDE-ODE reactive transport model of Sect. 2. It turns out that the new substitution method tends to smaller numerical errors and fewer Newton steps compared to the epsilon method. Finally, we discuss whether the substitution idea can be extended to "arbitrary" non-Lipschitz terms (Sect. 8), and we give a conclusion (Sect. 9).

The multi-species multi-reaction transport-reaction model
Let us consider a multi-species multi-reaction model with I ∈ N mobile (dissolved) species X i and J ∈ N kinetic reactions among the mobile species, as well as M mineral species Y j and mineral reactions. The vector of the mobile species concentrations is denoted by c aq = (c aq,1 , . . . , c aq,I ) T and the vector of the mineral species concentrations by c min = (c min,1 , . . . , c min,M ) T . The aqueous reactions and the precipitation-dissolution reactions are We assume that all the stoichiometric coefficients σ i j , σ i j , τ i j , τ i j are nonnegative integers, and we set We denote the aqueous reaction rates and the mineral reaction rates by r aq, j , j = 1, . . . , I , and r min, j , j = 1, . . . , M. The mass balance equations for the advection-diffusion-reaction problem then read i = 1, . . . , I , with porosity , and where (1) was implemented. If we want to take into account the evolution of the porous medium, we have to do this in the specific surface (e.g., by (2)), but equally one should take into account that the porosity is evolving with where ρ j are the molar densities of the pure minerals. This volume balance can be deduced from deriving the equations from upscaling by periodic homogenization [25,30] and is also used in [1], for one precipitation reaction only. Note that clogging may be possible but is not to be considered here; we focus on situations where mineral amounts change between zero and small positive amounts. Note that the assumption (6) leads to a back coupling of the reaction regime to the Darcy flow field q.
Our model with reaction rate vectors r aq for the aqueous reactions and r min for the mineral precipitationdissolution reactions then reads on a domain ⊂ R n , with T = (0, T )× , c, c min ≥ 0 (componentwise) and with some boundary conditions for c aq not specified here. In the following (cf. (12)-(18)), we will explain all the terms in our system (7)- (11). For an explanation of the entity, w see (18) and the subsequent text. In (7), L = (L 1 , . . . , L I ) T stands for an operator modeling diffusion/dispersion and advection, such as The product of vectors is meant componentwise, i.e., ( A r min ) j = A j r min, j , j = 1, . . . , M. As motivated in Sect. 1, we want to consider the surface model For what follows in the next chapters, it turns out that we may generalize this assumption to with α i ≥ 0 and A 0,i > 0. The case α i = 0 corresponds to the constant surface and for α i ≥ 1, the law (15) is locally Lipschitz continuous. We are mainly interested in the case (in particular, the case α i = 2 3 ) where the term (15) is not locally Lipschitz continuous at c min,i = 0. The aqueous reaction reaction rates in (7)-(11) are assumed to follow the law of mass action (LMA) with-for the sake of simplicity-ideal activities: The precipitation-dissolution reactions also follow the LMA, but with constant activity for the mineral species: j = 1, . . . , M, with precipitation and dissolution constants k p,i , k d,i > 0 and with w j ∈ H (c min, j ) with the Heaviside "function" or Heaviside graph By this Heaviside term in (8)- (9), the dissolution of minerals is switched off whenever a mineral concentration impends to become negative. This is an established method to implement the nonnegativity of constant-activity species (e.g., [1,[9][10][11]14,15,29]). In order to motivate this term, let us assume for a moment that we would replace r min, j by in (8), or, what is the same, to set w j = 1, j = 1, . . . , M, instead of (9). Then, if the fluid is undersaturated with respect to some mineral, i.e., ρ min, j (c aq ) < 0, and if the mineral concentration c min, j is zero, then negative mineral concentrations would be imminent-at least if the surface area term is non-zero, cf. (14). However, thanks to the Heaviside term, any solution of (7)-(11) is nonnegative. A proof of nonnegativity of solutions for the case (14) is standard and can be found for example in [10,29]; this proof easily carries over to non-constant surface laws like (15). There are other established formulations to implement nonnegativity of mineral concentrations. One can replace (8)-(9) by a so-called complementarity formulation [11] c min, j > 0 and ∂ t c min, j − A j (c min, j ) ρ min, j (c aq ) = 0 or c min, j = 0 and ∂ t c min, j = 1, . . . , M. It is well known that a complementarity condition can be transformed into an equation (see, e.g., [4,12] for the handling of mineral reactions at equilibrium), And another formulation can be found in [3]. An investigation in [10] showed that all these formulations are equivalent in some weak sense. For the model with a constant surface, (14), i.e., (7)-(11), (14), (16)- (18), under the simplification of a constant porosity instead of (6), existence of a global solution in a certain function space was proven in [9,10], provided that some mild assumption on the stoichiometry holds, uniqueness was proven in [10]. For a constant surface model restricted to one mineral reaction and without aqueous reactions, analysis (existence, uniqueness, convergence) can be found in [14,15,29], and such a model coupled to porosity and permeability changes is considered in [1]. For the model, we are mainly interested in, the model with the variable surface (15) in the specific form (7)-(11), (15)- (18); however, we could not find any mathematical analysis in the literature.
Let us discuss how a discretized version of the model could be solved numerically. While the Heaviside formulation (8)-(9) is well suited for analysis [1,9,10,15,29], we prefer the minimum formulation (21) for numerical solution. A discretization of (21) with the implicit Euler method reads where the subscript "old" for c min,old, j ≥ 0 denotes values from the previous time step. While f j is not very smooth (for the surface model (14), it has a kink where −c min,old, j − t A 0, j ρ min, j (c aq ) = 0), it is well known that such equations can be solved with the so-called semismooth Newton method [4,12]. From a practical viewpoint, the semismooth Newton method is in fact the same as the classical Newton method but applied to a "semismooth" function. 1 As for the classical Newton iteration, the semismooth Newton method is locally quadratic convergent for strongly semismooth functions-definitions [27] and proofs can be found in [6,Chapter 7]. Let us mention that other authors use the Heaviside formulation for discretization and apply a regularization to the discontinuous term [15] or event-driven algorithms detecting the moment when a trajectory meets a discontinuity [1]. We prefer to use the minimum formulation for discretization (21), (22) since it does not contain discontinuities.
In order to make (22) meaningful also for negative arguments c min, j for exponents α j < 1 (this might be reasonable in view of Newton's method, where iterates might be negative), we should discuss how to interpret a power with negative argument and fractional exponent. Let us define for the whole paper that This specific setting has another benefit. We have already pointed out that in general the Heaviside term, or, equivalently, the first of the two branches in (21) and in (22) is essential to avoid negative physically meaningless solutions. However, if we focus on the surface law (15) in combination with the setting (23), an adaptation of the nonnegativity proof in [10] shows that (8)-(9) might be simplified by setting w j ≡ 1 in our model, i.e., we may replace (8)-(9) by 2 -but only for the variable surface model (15), and not for the constant surface model (14). Analogously, for the surface law (15) and with the setting (23) and using c min,old, j ≥ 0, one can easily check that the first branch in (22) is dispensable and (22) may be replaced by

Batch models (ODE models)
In this section, we want to fix two simple ODE models as simplifications of the PDE-ODE model of the previous section. They will help us to understand the behavior of the full model from Sect. 2 and its discretization.
We assume that there is no transport, or that the initial values are space independent and the boundary conditions of (7)- (11) are such that solutions are constant in space. We assume that (6) can be approximated by = const and set this constant to one for the sake of simplicity. We consider the scenario that there are just two aqueous species and one mineral and one reaction The model from Sect. 2 then reads where c i is short for c aq,i . In view of (23), we may simplify this to We will focus on situations where i.e., there is no mineral at the initial state, and where the fluid is oversaturated, i.e., k p c 1 (0)c 2 (0) − k d > 0 (precipitation exceeds dissolution), e.g., Our system reads If we consider the even simpler chemistry instead of (25), we accordingly obtain the ODE model Since c aq + c min = const = c aq,0 (cf. (26)), we may eliminate c aq and get the ODE with

Non-uniqueness of solutions
Let us consider the ODE model (30)-which is equivalent to (29)-with initial condition (26). At the first glance, we get the impression that, since c min (0) = 0 and since κ 1 − κ 2 c min (t) ≈ const for small values of t > 0, the solution should behave like the solution of It is well known that this ODE admits infinitely many solutions t * ≥ 0 arbitrary. And indeed, one can show that also (30) with (26) admits infinitely many solutions. See Appendix A for a construction of a solution for (29)/(30) which is strictly positive for t > t * = 0. Besides, c min ≡ 0 is obviously another solution. It has to be expected that the solution is also non-unique for more complicated reactive networks, if the initial value for a mineral is zero and the fluid is oversaturated with respect to this mineral. Of course, the non-uniqueness issue is not only relevant for zero initial values for the mineral, but also if a mineral concentration drops to zero after some time and then some changes in the aqueous concentrations trigger precipitation. Concerning the PDE model from Sect. 2, we can expect that at least for constant-in-space initial values and certain boundary conditions (homogeneous Neumann), the non-uniqueness effect should carry over. Beyond that, we also can expect non-uniqueness for the PDE model (7)-(13), (16)-(18) for non-constant initial values: Suppose one mineral has zero concentration all over the computational domain at time t 0 = 0, and its concentration rises to a non-zero function c min, j (t) after some time t > 0 (for example, due to infiltrating aqueous species). Then also a second solution exists with mineral concentration c min, j remaining at zero. This second solution is the solution of the problem (7)-(13), (16)- (18), modified by omitting the j-th mineral reaction from the system. The following questions arise.
How to interpret this non-uniqueness. Which is the physically relevant solution? Second, what are the consequences of this non-uniqueness for numerical schemes? Which of the solutions is followed by the numerical scheme, are there any numerical difficulties coming along with the issue of nonuniqueness?
Third, how to overcome these numerical difficulties in a simple way? Concerning the first question, we postulate that in the scenario under consideration, i.e., a mineral concentration at zero in an oversaturated fluid, the solution which rises at once, without any delay time, is the relevant one. A justification could be given by considering models which have a closer look at the nucleation of minerals, such as [27]. In an oversaturated fluid, first nuclei should form with a certain non-zero probability per time, i.e., the expected value for the precipitated amount of mineral should immediately be-if tiny, but-positive, without delay. One might also mention here that the solution which rises without delay is distinguished from the other solutions by being the limit of solutions with initial value c min,0 > 0 for c min,0 → 0, i.e., it is stable in this sense. In several computer codes, the condition c min,0 ≥ 0 is replaced by c min,0 ≥ in order to evade any numerical difficulties arising from the non-uniqueness (cf. Sect. 5). Hence, these codes approximate that solution which rises without delay. In any case, the solution which rises without delay is the one we want to compute in this paper.
In Sect. 5, we tackle the second question. We analyze the simple ODE model (cf. Sect. 3) in order to demonstrate that the non-uniqueness issue causes numerical difficulties for the discretized problem.
The main purpose of this paper is to answer question 3. We make a proposal how the numerical solution, which rises immediately, can be computed in a robust and simple way, overcoming the numerical difficulties. See Sect. 6.

Numerical challenges in the discretized model
Let us consider a discretization in time of our model (26)- (28). We focus on Euler's implicit method, since this is probably the most used timestepping for reactive transport problems. See a short discussion of explicit timestepping methods at the end of this section.
Denoting the unknowns at the new time step by c 1 , c 2 , c 3 and the given values at the previous time step (or, the initial values in case of the computation of the first time-step) by c 1,old , c 2,old , c 3,old ≥ 0, we obtain the root-finding problem with Hence, a first difficulty is that the Jacobian is unbounded for c 3 close to zero. (Remember that our initial value, and so probably also the starting point for Newton's iteration, should be at zero.) So Newton's method or the Jacobian would have to be modified for iterates close to zero. Let us assume that we have overcome this first technical difficulty in one way or another. Let us turn to the first time-step, i.e., c 3,old = 0. Then one solution of (31) is obviously This is the solution which remains at zero, i.e., the solution which we do not want to find (see discussion in Sect. 4).
Another solution with c 1 < c 1,old , c 2 < c 2,old , c 3 > c 3,old = 0 could be computed explicitly by substituting the first two equations of (31) into the third. This is the solution that our numerical scheme should produce. For the moment, let us estimate this solution by exploiting that the rate ρ(c 1 , c 2 ) should not differ very much from ρ(c 1,old , c 2,old ) (if the time step is small, the supersaturation assumed for the initial condition should last over some time steps). Then we find (remember c 3,old = 0) That means that the two numerical solutions (32) and (33) of the first time-step are very close, and so it is not so easy to make sure that Newton's method converges to the wanted solution. See Fig. 1, left part. The bold line shows the root-finding problem for the first time-step (oversaturated scenario with initially no mineral present, c 3,old = 0), the physical solution (33) is displayed in green and the unphysical (32) in red. In this sketch only, the component f 3 as a function of (only) c 3 is displayed. The red line segment indicates: if the initial guess or any Newton iterate is here, we cannot expect that Newton's method converges to the physical solution.
Supposing that we have solved somehow the first time-step, let us consider now later time steps, i.e., we consider c 3,old small but strictly positive now. The graphs of the functions f 1 , f 2 , f 3 are just lowered by some constant c 3,old compared to the first time-step (dashed line in Fig. 1, left part). The unwanted solution at zero, which appeared in the first time-step, is not any more present. However, Fig. 1  Let us now discuss some straightforward techniques to evade the non-uniqueness issue and the numerical difficulties discussed above. Literature about the non-uniqueness issue for surface-controlled mineral precipitation reactions seems to be almost nonexistent. We found a report by [24] where this issue is discussed. One proposal in [24] is to just use another surface-to-mass relation which is Lipschitz continuous, or, to use a regularization of (15), which behaves like A 0 c α for large c and likeÃ 0 c withÃ 0 = A 0 c 1−α * for small nonnegative c. The regularization (34) is indeed Lipschitz continuous, so the issue of non-uniqueness vanishes. However, for c = 0, it hasÃ(0) = 0, i.e., c ≡ 0 constant in time would be the (unique) solution. So this regularization chooses a solution from the bunch of solutions, but the "wrong" one. The user would have to hope that numerical errors (e.g., in a PDE setting) will trigger the rise of the mineral concentration in the oversaturated case. We think that if a modification or regularizatioñ Instead of a modification of the relation A = A(c), another straightforward idea is to prevent that c is exactly zero, i.e., to enforce that c ≥ for some given >0. In the setting of Sects. 2-3, this could be easily implemented: in the complementarity formulation (20) we just replace "c min, j > 0," "c min, j = 0" by "c min, j > ," "c min, j = ." Correspondingly, in the Heaviside formulation, replace (9) by and in the minimum formulation, replace (21) by which is discretized by Euler's implicit method replacing the original equation The right part of Fig. 1 sketches the root-finding problem for the modified problem (36). The three lines represent three different choices of , here denoted by , , . The wanted root of the unmodified problem (37) is depicted in green. The value = is too small and can lead to nonconvergence of Newton's method (red part of the line).
For the choice = , the root is perturbed (black circle), i.e., we get only an approximate solution of (37). Note that we might also choose = j individually for every mineral reaction. If one chooses this "epsilon approach," the question arises how to set or j . Some rather crude analysis of (36) and right part of Fig. 1 indicates that j should meet the condition: to obtain/expect convergence for Newton's method in the case α j = 2 3 ; see Appendix B in [13] for details. For slightly larger j , in fact for j begin larger than the zero of (37), i.e., for j larger than one has to expect that the root of the function (36) is not equal to the root of the original equation (36), but only an approximation to it (see the shift from green to black circle in Fig. 1), right part. This range for j , between (38) and (39),where no perturbation is caused, is very small ( t 3 is much smaller than t) and this range depends not only on t and given rate constants, but also on dynamically changing variables (the aqueous concentrations-they vary not only from time step to time step by, let's say, O( t), they also vary within the Newton iteration. Hence, it is unlikely or impractical to hit this range, and so we have to accept that the introduction of an > 0 (Eqs. (36) instead of (37)) will influence numerical results. We present numerical tests on the choice of j and its effect on the accuracy in Sect. 7. Hence, the enforcement of c min ≥ is a modification, an approximation of the original model. The same, of course, holds true for any regularization of the non-Lipschitz surface term such as (34). The main issue of this publication is to point out that no modification is necessary. See Sect. 6 for the proposal of a substitution.
In this section and in the rest of this article, we focus on implicit timestepping methods, since it is well known that explicit schemes for diffusion problems would suffer from a severe restriction on the time step size. However, some reactive transport codes use a splitting of a time step into a transport substep and a reaction substep-note that such a splitting either introduces an additional error ("splitting error") or requires an iteration between transport and reaction ("sequential iterative approach"). However, if the reaction step is discretized with an explicit scheme such as Euler's explicit method-this may be an option if the mineral reactions are rather slow-then we do not have to solve a problem of the type (31), instead, we have the explicit computation For an initial value of the mineral c 3,old = 0, we would obtain c 3 = 0 even if ρ(c 1,old , c 2,old ) > 0, i.e., our scheme would pick the unphysical solution. One would have to hope that numerical errors trigger the rise of the solution somehow, but we would not call this a robust scheme. Similar as for the implicit scheme, an ad hoc way to resolve this issue would be the introduction of an >0 and to set c 3 = max{c 3,old + t A 0 c 2/3 3,old ρ(c 1,old , c 2,old ) , } and to postulate the initial value being equal to or larger than , i.e., a modification (perturbation) of the original model, or the method proposed in the next section.

A substitution
In this section, we apply a substitution to derive new model equations. The new model will not suffer from any non-uniqueness issues (cf. Sect. 4 and Appendices A, B) and the numerical challenges explained in Sect. 5 will disappear, without any perturbation of the model by an epsilon.
Let us present the idea for the case α = 2 3 , which seems to be the most relevant one. We apply the substitution ξ = c to the ODE model (28) (and componentwise to the PDE model (7)- (11)) and obtain, for c min > 0 where we have used the ODE (8)/(28) with (13) in the last step. The ODE model (28) now reads for the unknowns c 1 , c 2 , ξ. We state that the non-Lipschitz term in the ODE (28) For c min,i > 0, we obtain where we have used the ODE (8) with (15) in the second but last step. Replacing the ODE (8) by (44) and the unknown c min,i by ξ i , the non-Lipschitz term has vanished from our model. The resulting system reads What we can easily see is that in a setting of oversaturation (i.e., positive r min,i ), a solution with c min,i (t) = const = 0, which is ξ i (t) = const = 0, on an interval [0, t * ) is no longer possible. It is in fact possible to prove uniqueness of the solution. Since for the full PDE model-Eqs. (45)-(49) together with boundary conditions-the proof is rather technical and also very similar to the uniqueness proof in [10], we only present a uniqueness proof for the simpler ODE model (i.e., we omit the term L c aq in (45)) in Appendix B.
The substitution we made is an equivalence transformation as long as c min,i > 0. Hence, different from the ad hoc proposals from Sect. 5 (cf. (34) and (35)), our transformed system has precisely the same behavior as the original system. And for initial value c min,i (0) = 0 and oversaturation r min,i , ρ min,i > 0, we obviously obtain a solution which rises at once. So the transformation picks the "wanted" solution (see end of Sect. 4) from the bunch of infinitely many solutions. The other solutions disappear.
For numerical simulations, we replace (46)-(47) by min{ξ j , ∂ t ξ j −Ã 0, j ρ min, j (c aq )} = 0 and discrete in time (21), (22). Note that the mineral source term in the aqueous species equations (which was ∂ t c min, j in the old model) now reads or can be expressed as follows: However, in the discretized scheme, it is favorable to use the first version, ∂ t (ξ , for the sake of mass conservation, i.e., the term is discretized as 1 ]. We emphasize that-apart from the mineral source term in (45)-the variable surface model after the transformation, (45)-(49), looks like a constant surface model. And therefore, it is not admitted to replace (47) by w = 1 or to remove the first branch from the minimum function, except if there is some a priori knowledge that c min,i (t) >0 (or that the fluid is oversaturated) for all i on the whole time interval under consideration (0, T ].  By applying the substitution, the model with variable surface obtains a structure very similar to models with fixed surface, for which analysis is available. The only difference, if (45)-(49) is compared with the constant surface model in [9,10], lies in the mineral source term in the PDE. Apparently, the substitution establishes possibilities to extend analysis from [9,10] (existence of global solution) and possibly also from [14,15] (convergence of discrete schemes, for fixed-surface models of one aqueous and one immobile species) to the variable surface model.

Numerical simulations
We have implemented the new substitution method and, for comparison, the -method, in matlab. The code can be used both for the ODE setting (batch problem) as well as for the PDE-ODE setting (reactive transport). We use linear finite elements in space and implicit Euler time stepping. The nonlinear system is solved with Newton's method. We implemented Armijo's rule. For the sake of simplicity, we assume that the porosity changes are moderate so that we can approximate (6) by = const so that we can use a constant flow field in Sec. 7.2.

Batch model
At first, we consider the ODE model with three species (26)-(28), i.e., we have a binary precipitation-dissolution reaction (25), we have an oversaturated fluid, we start from a zero mineral concentration, exponent α = 2 3 , and A 0 = 1. We want to test our new substitution method (cf. Sect. 6) and compare its robustness and its accuracy with the -method (cf. Sect. 5). At first, we test how large has to be chosen for the -method so that Newton's method converges. The results are depicted in Table 1. We consider a time interval [0, 1], and we stop the Newton iteration when the residual is at most tol Newton = 10 −12 in the L ∞ -norm. As the considerations of Sect. 5 predicted, the tests confirm that the smallest for which Newton's method converges behaves as O( t 3 ). Note that for the -method, we have to start at c min,0 = instead of c min,0 = 0. On the other hand, for the new substitution method, we do not need any and we can start at c min,0 = 0.
Next, we measure the error of c min at time t = 1 for the -method and for the new substitution method for the same scenario as in Table 1. We consider values of t ranging from 10 −1 to 10 −3 . As a reference solution, we take the numerical solution for t = 2 · 10 −5 . The results are depicted in Fig. 2.
The red solid line is for the -method when we choose = ( t) 3 . Note that a computation with = 0.25 · ( t) 3 fails (while = 0.3 · ( t) 3 still works) in this test. This fits very well to the findings of Sect. 5 and [13, Appendix B] that should be chosen as ∼ max(A 0 tρ) 3 , more precisely as Fig. 2 Comparison of the substitution method (green) with the -method (red) in an ODE setting. For the red-bold line, we chose = ( t) 3 . The substitution method leads to more accurate results. Initial value for the mineral c 3 : for the substitution method, we are allowed to take zero; for the substitution method, we must take the inaccurate value . The black line uses the substitution method with the inaccurate initial value . Dotted lines: fixed where ρ is the value of the precipitation-dissolution rate, and the maximum is taken over those points in space and time where the mineral concentration rises from zero to positive values (it is ρ ≈ 1 in the scenario under consideration). The green line is for the new substitution method which works without the choice of any (and with = 0). Since we use implicit Euler timestepping, we observe first-order convergence for both methods. However, the new substitution method produces errors which are two powers of ten smaller than the -method. If we choose constant , i.e., independent of t, then we obtain no convergence to the exact solution at all (see the dotted red lines) for t → 0. This is not surprising, since we have perturbed the original problem by a fixed term. Let us discuss the question why the new substitution method is more accurate than the -method. The substitution method needs no perturbation, we can start our simulation exactly with the initial value c min (0) = 0, while for the -method the smallest possible value is c min (0) = , and this initial error of course spoils the whole solution. To investigate this effect in more detail, we also perturbed the substitution method by using the initial value c min (0) = , see the black line. It turns out that the solution of this variant is still a bit more accurate than the -method (red line). Hence, it seems that the disturbance of the initial value is not the only reason for the larger error. It seems that the term c 2/3 min in the right-hand side of our ODE (which is present using the -method, but which is not present when using the substitution) causes larger numerical errors in our timestepping scheme.

Reactive transport model (PDE-ODE model)
For the next test, we want to apply the new substitution method in a PDE setting. We consider a spatial domain = (0, 1)×(0, 1). We use linear finite elements on a Friedrichs-Keller mesh with mass lumping for the source terms and the accumulation terms, i.e., for the time derivatives. We have again two aqueous species and one mineral, and we consider the asymmetric binary reaction 2 X 1 + X 2 ←→ Y with rate coefficients k p = 2 · 10 4 , k d = 4 · 10 1 , α = 2 3 , and A 0 = 1. We use zero flux boundary conditions (−d∇c aq,i + qc aq,i ) · n = 0 at the inflow (=left) boundary and homogeneous Neumann conditions n · ∇c aq,i = 0 at the other boundaries. Our initial condition is where χ M is the characteristic function of the set M. Hence, the domains where the two aqueous species are present in the beginning do not overlap. A time series for the numerical solution is depicted in Fig. 4. It is computed using the new substitution method. We use t = 5·10 −4 and 128×128 mesh points. Due to diffusion (d = 0.15), the species mix, and where k p c 2 aq,1 c aq,2 > k d , the precipitation of the mineral sets in. This is the case at about t = 0.0085. The red line indicates where precipitation r prec = k p c 2 aq,1 c aq,2 and dissolution rate r diss = k d are equal, i.e., inside the red line the mineral concentration increases and outside it decreases. Because of advection (Darcy velocity q = (2, 0) T ) the concentrations of the aqueous species move to the right, and so does  Fig. 4 we see the precipitation of the mineral and its effect as a sink term on c aq,1 , c aq,2 (inside the red line). In the fourth and fifth the shape of the contour lines is influenced by the mineral dissolution serving as a source term for c aq,1 , c aq,2 . The test shows that the substitution method works also in a PDE setting, in a scenario that covers both the rise of the mineral concentration from zero and the decrease of the mineral down to zero.
We confirmed that the new substitution method converges of first order in t, see Fig. 3 (green line). For this test, we used 32×32 mesh points and computed the L ∞ -error at t = 0.32 for t ∈{16, 8, 4, 2, 1, 1 2 , 1 4 }·10 −3 . As a reference solution, we used the numerical solution with t = 2 −6 · 10 −3 on the same spatial mesh. (By extrapolation, we can estimate that the reference solution carries a temporal error of about 5.5·10 −5 .) We also tested the -method. For , we tried values 10 −i for i = 1, 2, . . . , 8, see the black lines in Fig. 3 from top to bottom. For a given > 0, for t too large, Newton's method does not converge; therefore, the black lines stop at a certain t-this is similar as in the ODE test setting of Sect. 7.1. Regarding the accuracy, Fig. 3 shows that there are some couples of and t for which the numerical error is similar as for the new substitution method (e.g., t = 2·10 −3 with = 10 −3 ). However, if for a given t-and the scenario under consideration-the value of is chosen too large then the error is larger than the error of the new substitution method. Hence, in order to obtain convergence and small numerical errors, the choice of the numerical parameter (or j , j = 1, . . . , M, if several mineral reactions are considered) is delicate. We consider it to be a major advantage of the substitution method that it does not require such a numerical parameter. Figure 5 shows the average number of Newton steps per time step. For the -method, we took some effort to take an optimal value of (at least up to a factor of two) for every t, i.e., we determined the smallest of the shape = 2 i · 10 −5 , i ∈ N, for which Newton's method converges. t by a factor of two, we roughly may reduce by a factor of about 2 6 (or sometimes 2 5 ). Is this a deviation or contradiction to the observations in the ODE test case (Sect. 7.1)? In fact, a major difference of the PDE-ODE test case to the test case of Sect. 7.1 is that we now start from undersaturation (at t = 0). In that moment when-at a certain point of the domain-saturation is reached and the mineral concentrations starts to rise from zero, then the rate ρ is zero at this point, and so the formula (51) seems difficult to interpret. If we assume that during the time step when-at a spatial point under consideration-the value of ρ rises from undersaturation (ρ < 0) to oversaturation  Symbol "x" for 32 × 32 grid points, symbol "•" for 64 × 64 grid points, symbol "+" for 128 × 128 points. In most cases, the substitution method requires fewer Newton steps than the epsilon method (ρ >0), the value of ρ is about ρ ∼ O( t), then we obtain O( t 6 ) from (51), which fits well to the numerical findings mentioned above.
The Newton iteration was carried out until the residual reached tol Newton = 10 −12 . Figure 5 shows that the number of Newton steps, computed with the optimal for the -method (in red), and with = 0 for the new substitution method (in green), is (in almost all simulations) smaller for the substitution method. This carries over to the overall numerical effort, since the problems to solve are very similar for the two methods (the sparsity of the matrix is slightly better for the new substitution method). We think that the smaller number of Newton steps for the substitution method compared to the epsilon method is caused by the elimination of one of the nonlinearities, c 2/3 min , from the system. Let us mention here that Armijo's rule was sometimes invoked with the epsilon method (in 1% of the Newton steps) but never for the -method. In the tests, we described in this chapter we never obtained a "nonconvergence" for the substitution method, but it occurred frequently for the -method when >0 was chosen too small. Finally, let us mention that for very large time steps, such as t = 16·10 −3 and, e.g., 64×64 grid points, we have 2 15 ·10 −5 ≈ 0.3 that means that the numerical results obtained with the epsilon method are completely wrong, since we have to raise the intended initial value c min (0) = 0 to c min (0) = ≈ 0.3; for comparison: the maximum value that the true solution attains after some time is about the same magnitude. With the substitution method, we can use c min (0) = = 0. Now let us consider another scenario. All the data are exactly the same as in the previous scenario, except that we change our initial value (52) to i.e., there is a domain where both aqueous species are present at t = 0. As a consequence, precipitation sets in at t = 0 at a rather high rate; in the first time-step, we can already expect a precipitation with rate ρ ≈ k p c 2 aq,1,0 c aq,2,0 −k d ≈ 2 · 10 4 in the overlap domain (instead of ρ ∼ O( t) in the previous setting). Hence, if we apply our estimator (which we originally developed for the ODE setting) to this scenario, we expect that the optimal for the -method should be about ≈ 8 3 , and not ∼ O( t) 6 as in the previous scenario. In fact, the numerical tests confirm such a behavior for small t, see Table 2.
So, for example, for t = 10 −5 , we must use ≥ 0.003 for the epsilon method, i.e., we start with an initial error of (at least) 0.003 in the mineral concentration. (For comparison, the maximal value that the mineral concentration reaches after some time is about 1.3, so the initial error is less than one percent.) However, a simulation with this timestep size requires 75,000 time steps to reach time T = 0.75, the time when the aqueous species are advected out of the computational domain. Hence, one might want to use larger time steps. However, with t = 8 · 10 −5 (9375 time steps), we must use at least = 0.13, and so the initial error is already about ten percent of the expected mineral amount. Hence, for larger time steps, the -method cannot be used, i.e., if it is used, it requires such a large that the numerical results are meaningless. On the other hand, the new substitution method has no problem with Table 2 The -method for a scenario where the precipitation happens at a large rate (initial condition (53) instead of (52)) t z = A 0 t ρ pred fail ok 65 × 65 grid points. The value z = A 0 t (k p c 2 aq,0,1 c aq,0,2 − k d ) ≈ 2 · 10 4 · t indicates how well the initial precipitation pulse is resolved by the time step size. The optimal value for , as predicted by pred = 8 27 · z 3 . For = fail , the method did not converge, for = ok , it converged. Since ∼ O( t), for larger t, the -method requires such a large that the numerical results become meaningless such a large time step size, or even much larger ones. For the much larger timestep size t = 2.5 · 10 −3 (which corresponds to 300 time steps) and = 0, it runs with an average of 3.05 Newton steps per time step, and with t = 10 −2 (i.e., only 75 time steps), it runs without difficulties with 4.12 Newton steps per time step. Of course, the initial fast precipitation in the overlap domain is not any longer resolved with these larger t. However, all the other processes (advection, diffusion, precipitation-dissolution in other places or for larger time t) are well resolved. The numerical solution obtained by the substitution method is depicted in Fig. 6.
To summarize, if strong precipitation processes occur, the necessity of choosing a sufficiently large for convergence of Newton's method acts as a kind of constraint for the timestep size. For large time steps, the method requires such large that the numerical results are polluted by a large error, making the numerical results meaningless. The new substitution method, however, runs stable even for large t (and = 0).

Extension of the substitution method
The question arises if the issue on non-uniqueness/the presence of a non-Lipschitz term can be avoided if "arbitrary" surface dependences of the rate beyond (15) are used in the model. For example, one might be interested in a model where the surface behaves like (15) for small mineral concentrations c min,i , but where the surface area is bounded for c min,i → ∞. This might be achieved by (a) a surface model of "product type," with the non-Lipschitz factor c α i min,i and a smooth factor f (c min,i ), or (b) by a concatenation with a Monod term, For model extension (a), we try the same substitution which worked for the standard model (15), i.e., (43). We obtain (compare (44)) Hence, the removement of non-Lipschitz terms was successful, for arbitrary α i ∈ (0, 1). Let us turn to model extension (b). Since for small mineral amount (54) is close to (15), we might try the same substitution which worked for (15) and for case (a), i.e., (43). We obtain Hence, the removal of non-Lipschitz terms was successful for α i ∈ [ 1 2 , 1) but surprisingly failed for α i ∈ (0, 1 2 ). However, the most important scenario of α i = 2 3 is covered. The question arises whether there is any other substitution which removes the non-Lipschitz term also for α i ∈ (0, 1 2 ), or, even more general, for an "arbitrary" law for a function F i which is not specified and which is not Lipschitz at zero. In order to find the matching substitution, we use the ansatz . Without loss of generality, we assume s i (0) = 0. Following (55), we obtain Hence, our task is to find s such that the mapping This ODE is autonomous and might be solved if F i is sufficiently simple. For the model (b) discussed above, i.e., for F i (c min,i ) = c α i min,i /(K i +c α i min,i ), we find the solution of (56) : R + 0 → R + 0 are monotonically increasing functions. We have found a substitution that removes the non-Lipschitz term for (54) not only for α i ∈[ 1 2 , 1), but also for α i ∈ (0, 1): Indeed, the mineral ODEs then read while the PDEs contain a source term We have to realize that for the surface model (54) with α i ∈ (0, 1 2 ), the situation is not as nice as for the original surface model (15): We do not have the inverse transformation c min,i = s i (ξ i ) in an explicit form 4 , and so we would have to compute c min,i iteratively whenever it is needed-that is for the PDE source term (58)-or to use look-up tables and interpolation.
Let us mention that a surface model which behaves like (15) for small mineral concentrations c min,i and where the surface area is bounded for c min,i → ∞ is a cut-off version of (15), for some c i, * > 0. This surface model is of type (a), and so for this surface model, the original substitution (43) works, and there are no difficulties with the inverse transformation as for (57). In fact, we have implemented this version in our code. Of course, the applicability of the substitution (43) remains true if we replace the kink in (59) by some smoother transition. Let us summarize. The substitution technique can be generalized to many surface area models beyond (15), to remove a non-Lipschitz term. However, depending on the chosen model, it might happen that a suitable substitution for which also the inverse function is explicitly given, is not easily at hand.

Summary
We have pointed out that the typical reactive transport model with mineral surface A ∼ c 2/3 min is ill posed in the sense that it admits non-unique solutions whenever a mineral concentration rises from zero to positive values. We demonstrated by theoretical investigations and by numerical tests that convergence problems of Newton's method have to be expected for a numerical scheme based on the ill-posed model. We compared two strategies to remove the ill-posedness and the nonconvergence issue: the-what we call-epsilon method, which forces any mineral concentration to be larger than or equal to > 0, and the substitution method. At least the latter is new to our best knowledge. To summarize in short, we propose to use the new substitution method. Disadvantages of the epsilon method are -A parameter > 0 has to be chosen. If it is chosen too small, then nonconvergence of Newton's method is imminent. The larger it is chosen, the more error is introduced in the numerical solution. Let us mention that the necessity to carefully choose a parameter also is true for other potential regularizations of the model. In practice, the parameter is probably chosen by trial and error. If a simulation is started with another timestep size, then also should be changed. Questions occur for numerical codes with adaptive timestepping. -While, in some scenarios, tests show that the additional numerical error of the epsilon method may be tolerable, we found scenarios (transition from zero to nonzero mineral concentration by fast precipitation, precipitation that is poorly resolved by the timestep size) where the accuracy of the numerical solution becomes extremely bad or very short time steps have to be used. -In most of our tests, the number of Newton steps was larger for the epsilon method than for the substitution method. This means a larger computational effort for the same spatial and temporal resolution.
We think that the main advantages of the new substitution method are -no parameter has to be chosen, -it does not introduce additional numerical error, -it might influence the number of Newton steps in a positive way, -the method is elegant and easy to implement, and no technical modifications of Newton's method are required, -and it allows to carry over proofs of existence and uniqueness of solutions for the constant surface model to the A ∼ c α -surface model.
An interpretation of the substitution and its performance might be that it seems to be advantageous to use grain size (radii) instead of mineral volume or mass in a macroscale model. We hope that the substitution is a contribution for the improvement of robustness and efficiency of reactive transport codes. It should not be too difficult to implement the substitution into existing reactive transport codes. For the future, it would be desirable to obtain experiences and comparisons also for more complex scenarios (parameters and concentrations varying over many magnitudes, or aqueous reactions modeled by equilibrium conditions, or adaptive timestepping).
Appendix B: Proof of uniqueness of solutions for the ODE version of the problem formulation from Sect. 6 Let us consider the reactive transport formulation after the substitution of Sect. 6 was performed, in an ODE version. That means that we consider the equations (45)-(49), but we omit the term L c aq in (45). In Sect. 6, we have claimed that the substitution leads to a problem that has a unique solution. We are going to prove that there is at most one solution in the following sense: Theorem For the proof, we need Gronwall's lemma in the following version (e.g., [28]): Let ψ, c 1 , c 2 be continuous and c 2 ≥ 0 and c 1 monotonically increasing. If ψ(t) ≤ c 1 (t) + t 0 c 2 (τ )ψ(τ ) dτ for all t, then ψ(t) ≤ c 1 (t) exp( t 0 c 2 (τ ) dτ ). Proof of the theorem The proof consists of two parts. In the first part, we estimate the difference of the transformed mineral concentrations ξ (t),ξ (t) in terms of the difference of the aqueous concentrations c aq (t),c aq (t). For this we use (46). In the second part, we estimate the difference of the aqueous concentrations and show that the difference grows at most exponentially in time. For this we use (45) and the result of part 1. Exploiting the initial conditions, it follows that the differences are zero on [0, T ]. Before, let us introduce the notation (compare (19)) For later use, we state that and we state that there are constants L aq , L prec , L diss , L pow (depending on the functions (60)) such that S r aq (c aq (t)) − S r aq (c aq (t)) ≤ L aq c aq (t) −c aq (t) = L aq u(t) , ρ prec (ξ (t)) − ρ prec (ξ (t)) ≤ L prec ξ (t) −ξ (t) = L prec v(t) , for all t ∈ [0, T ] where · is the Euclidean norm and where the power ξ β is meant componentwise. For (62) to hold, we have used the fact that r aq , ρ prec , ρ diss and the mapping x → x β i = x 1/(1−α i ) for α i ∈ (0, 1) are locally Lipschitz continuous, together with the fact that the values of c aq (t),c aq (t), ξ (t),ξ (t) are bounded due to our assumption (see definition of V ). Part 1 of the proof. We consider the i-th component of Eq. (46) for ξ i and forξ i and take the difference, and then we multiply this equation by v i . We obtain Since the mapping ξ i −→ w i ∈ H (ξ i ) is monotonically increasing, we know that the estimate −(w i −w i ) v i ≤ 0 holds and so the last of the three terms in (63) can be estimated by zero. By using (62) and |w i | ≤ 1, we obtain 1 2 ∂ t |v i | 2 ≤Ã 0,i (L prec + L diss ) u |v i | . By summing over i = 1, . . . , M, we obtain where we set C 1 = Ã 0 (L prec + L diss ). We integrate and use (61) to get v(t) 2  From this exponential growth limitation follows (Gronwall's lemma again) that u ≡ 0, and using (64) again, also v ≡ 0 holds, which completes the proof. Some final remarks on the proof: The proof also shows that the solution depends continuously on the initial value. The substitution does not remove all non-Lipschitz terms from the model, since the Heaviside function, which is non-Lipschitz and even discontinuous, is still contained. However, the mineral source term in (46)-(47) is one-sided Lipschitz (cf. (3)), and this is enough for a proof of uniqueness of solutions, as we have seen.
Note that the mineral source term in the PDE (cf. (50), second and third formulations) contains the factor ξ α i /(1−α i ) i which is not Lipschitz if α i < 1 2 . Hence, in some sense, we have removed a non-Lipschitz term but obtained another non-Lipschitz term instead. However, the proof of uniqueness above shows that this is irrelevant: for the proof-as well as for the numerical computation-we used the first of the representations in (50) with ξ i to the power 1/(1−α i ), which is Lipschitz for all α i ∈ (0, 1) under consideration.
Reference [10] shows how to extend such a uniqueness proof to a PDE setting.