SolveSAPHE-r2 (v2.0.1): revisiting and extending the Solver Suite for Alkalinity-PH Equations for usage with CO 2 , HCO − 3 or CO 2 − 3 input data

. The successful and efﬁcient approach at the basis of the Solver Suite for Alkalinity-PH Equations (SolveSAPHE) (Munhoven, 2013), which determines the carbonate system speciation by calculating pH from total alkalinity (Alk T ) and dissolved inorganic carbon ( C T ), and which converges for any physically sensible pair of such data, has been adapted and further developed to work with Alk T –CO 2 , Alk T –HCO − 3 , and Alk T –CO 2 − 3 . The mathematical properties of the three modiﬁed alkalinity–pH equations are explored. It is shown that the Alk T –CO 2 , and Alk T – HCO − 3 problems have one and only one positive root for any physically sensible pair of data (i.e. such that [ CO 2 ] > 0 and [ HCO − 3 ] > 0). The space of Alk T –CO 2 − 3 pairs is partitioned into regions where there is either no solution, one solution or where there are two. The numerical solution of the modi-ﬁed alkalinity–pH equations is far more demanding than that for the original Alk T – C T pair as they exhibit strong gradients and are not always monotonous. The two main algorithms


Introduction
Among all the aspects of the ongoing global environmental changes (climate change, ocean acidification, etc.), the solution chemistry of carbon dioxide (CO 2 ) is one of the best known. The related chemistry of the carbonate system in the oceans and other aqueous environments is well understood and routinely monitored and modelled. The equilibria between the carbonate system species involves four variables: [CO 2 ] (or equivalently the partial pressure of CO 2 , pCO 2 , or its fugacity, f CO 2 ), , [CO 2− 3 ] and [H + ] (or equivalently pH). The speciation, i.e. the determination of the concentrations of the individual species, therefore also requires four constraints. Two of these are given by the equilibrium relationships that characterise the equilibria between dissolved CO 2 and HCO − 3 on one hand, and between HCO − 3 and CO 2− 3 on the other hand, assuming that the respective equilibrium constants are known or can be calculated. Two more independent constraints are thus required to completely characterise the system. These two additional constraints are generally chosen among the four traditional measurables of the system (see, e.g. Dickson et al., 2007): (1) the total concentration of dissolved inorganic carbon (DIC), equilibrium constants adopted, pressure corrections applied, etc. Here, we do not focus on these aspects but on the design of algorithms that can solve the underlying mathematical problem with as little user input as possible. The aim is to reduce user input to the bare essentials: besides the fundamental information about temperature, salinity, pressure and the thermodynamic data, this ideally had to be any physically meaningful data pair only; the algorithm should be able to derive any other auxiliary information, such as root brackets or starting values for iterations, on its own.
In the companion paper (Munhoven, 2013), such autonomous algorithms with robust convergence properties for a wide range of environmental conditions are presented for usage with the Alk T -C T pair. Here, we are revisiting that approach, extending and adapting it so that the Alk T -CO 2 , Alk T -HCO − 3 , and Alk T -CO 2− 3 pairs can be processed with the same ease and reliability. For the sake of completeness -and with minimal details only -"recipes" for solving the other 11 explicit cases are provided in the Appendix. Alternative approaches can be found in the literature, such as in Zeebe and Wolf-Gladrow (2001) or Dickson et al. (2007). Dickson et al. (2007) also provide pathways for using triplets or quartets of input data, which only require the knowledge of one of the two dissociation constants or of their ratio, or none of them. That kind of approach is, however, not considered in this study.

Theoretical considerations
In the following, it is assumed that the temperature T , salinity S and applied pressure P are given and that adequate values for all the required stoichiometric equilibrium constants are available. It is furthermore assumed that the total concentrations of all the other relevant acid systems (borate, hydrogen sulfate, phosphate, silicate, etc.) are known.

Revisiting the mathematics of the alkalinity-pH equation
Cornerstone to the speciation calculation is the resolution of the following equation, which I call the "alkalinity-pH equation", as it derives from the definition of total alkalinity: i.e. Eq. (21) from Munhoven (2013). In this equation, [H + ] is the proton concentration expressed on one of the commonly used pH scales (total, seawater) and s is a factor to convert from that scale to the free scale. s depends on temperature, pressure and salinity of the sample and its value is close to 1 (typically between 1.0 and 1.3). The first term on the right-hand side is that part of the total alkalinity that is not related to the water self-ionisation: Alk nW ([H + ]) = i Alk A [i] ([H + ]), where i denumbers the acid systems resulting from the dissolution of acids A [i] whose dissolution products contribute to total alkalinity. For the purpose of this study, Alk nW ([H + ]) is partitioned into carbonate alkalinity, Alk C ([H + ]), and non-carbonate alkalinity, Alk nWC ([H + ]), since the relevant carbonate system parameters (the concentrations of CO 2 , HCO − 3 and CO 2− 3 and their sum, C T ) are all directly related to Alk C : Alk nW ([H + ]) = Alk C ([H + ]) + Alk nWC ([H + ]). Similarly to Alk nW , Alk nWC admits an infimum and a supremum which can both be derived from the total concentrations of all the acid-base systems considered. We denote these two by Alk nWCinf and Alk nWCsup , respectively. Equation (1) is thus formally rewritten as The carbonate alkalinity term writes, as a function of C T where K 1 and K 2 are the first and second stoichiometric dissociation constants of carbonic acid. The individual carbonate species fractions of C T can be expressed as a function of [H + ]: Accordingly, Alk C ([H + ]) may be rewritten in one of the following forms: In order to get a first idea about the complications that we might encounter for the solution of the three new data pairs, we start with an exploratory analysis using the Solver Suite for Alkalinity-PH Equations (SolveSAPHE) version 1.0.3 (Munhoven, 2013. The three panels in the upper row of Fig. 1 show the pH and the HCO − 3 and CO 2− 3 concentration distributions for a reduced SW3 test case from Munhoven (2013), with the C T range extending from 0 to 4 mmol kg −1 and the Alk T range from −1 to 3 mmol kg −1 only. The pH distribution from Fig. 1a is then used to derive the corresponding CO 2 (not shown), HCO − 3 and CO 2− 3 concentration distributions ( Fig. 1b and c). Since we intend to solve the alkalinity-pH equation for given Alk T , and either of [CO 2 ],  or [CO 2− 3 ], we furthermore produce the concentration isolines for the three species on a pH-Alk T graph (Fig. 1d, e and f). For the three latter, we first calculated Alk C from Alk T by using Eq. (2). Positive Alk C values were then used with Eqs. (7), (8) and (9) to derive the corresponding [CO 2 ], [HCO − 3 ] and [CO 2− 3 ], respectively. Blank areas represent the pH-Alk T combinations that lead to negative Alk C .
The V-or U-shaped isolines for HCO − 3 on the C T -Alk T graph and for CO 2− 3 on the C T -Alk T and on the pH-Alk T graphs show that the C T -HCO − 3 and the Alk T -CO 2− 3 pairs will not always provide unambiguous results. This is illustrated by the blue and black stars in Fig. 1c  As will be shown below, the SolveSAPHE approach of Munhoven (2013), which is based upon the use of a hybrid iterative solver safeguarded by intrinsic brackets that can be calculated a priori, can be easily adapted for the Alk T -CO 2 and Alk T -HCO − 3 pairs. According to the outcome of our preliminary analysis above, the Alk T -CO 2− 3 pair requires a more in-depth analysis. We show that it is nevertheless possible to diagnose the different cases that can theoretically be encountered and, in case there are two solutions, to derive bracketing intervals for each of the two and to isolate them efficiently. For each pair, we (1) establish the analytical properties of the modified alkalinity-pH equation; (2) derive brackets for the root(s); (3) develop a reliable and safe algorithm to solve the problem; (4) design an efficient initialisation scheme. The C T -HCO − 3 pair, which requires only a quadratic equation to be solved, is straightforward to diagnose a priori (see the corresponding recipe in the Appendix).
We will now in turn analyse the mathematical properties of the alkalinity-pH equation that results from the substitution of C T by the concentration of one of its individual species. (2) is written as in Eq. (7). Equation (2)

Root bracketing
Intrinsic brackets for the solution of Eq. (10) can be derived similarly to what is done in Sect. 5.1 in Munhoven (2013). The lower bound H inf can be chosen such that i.e. as the positive root of the cubic equation: Let us denote this cubic by P (H ). It is important to notice that P (0) = −2K 1 K 2 [CO 2 ] < 0 and P (0) = −(K 1 [CO 2 ] + K W ) < 0. The equation P (H ) = 0 has therefore one and only one positive root. Similarly, the upper bound H sup can be chosen such that i.e. as the positive root of the cubic equation which also has one and only one positive root for the same reasons as above.
The positive roots of these cubic equations can be found by adopting a strategy similar to that used for the cubic initialisation of the iterative solution in SolveSAPHE (Munhoven, 2013, Sect. 3 3. solve Q(H ) = 0 which has two roots and choose the one that is greater than H min .
In this particular case, it is, however, not necessary to solve these equations exactly as we only need approximate bounds of the root for safeguarding the iterations while solving Eq.
(2). For H inf , we may actually choose the H min of the first cubic which is lower than the positive root and thus sufficient.
Regarding H sup , it should be noticed that P (H ) = Q(H ) + (H − H min ) 3 /s. Accordingly, P (H ) > Q(H ) for H > H min and therefore the greater root of Q(H ) for the second cubic is greater than the positive root of that cubic. The greater of the two roots of Q(H ) is therefore a sufficient upper bracket and may be used instead of the exact H sup . Any bracketing root-finding algorithm can then be used to solve the modified alkalinity pH equation (Eq. 10).

Root bracketing
The lower bound H inf can be chosen such that i.e. as the positive root of the quadratic equation Similarly, the upper bound H sup can be chosen such that i.e. as the positive root of the quadratic equation Both equations always have two roots: one positive and one negative; their product is negative as indicated by the constant term. With the respective positive roots, we have again bounds for the solution of the modified alkalinity-pH equation and any bracketing root-finding algorithm can be used to solve it.

The Alk T -CO 2− 3 problem
Whereas any physically meaningful Alk T -[CO 2 ] or Alk T -[HCO − 3 ] concentration pairs will always provide one and only one [H + ] (or equivalently pH) value as demonstrated above, this cannot be the case for every Alk T -[CO − 3 ] pair, as can be deduced from Fig. 1c and f. On one hand, there are two compatible C T values, and equivalently two pH values, for most Alk T -[CO 2− 3 ] pairs. This little-known fact was already documented in the 1960s (see, e.g. Deffeyes, 1965) 1 . On the other hand, there are also Alk T -[CO 2− 3 ] pairs that do not allow for any solution, as they lead to negative carbonate alkalinity. To our best knowledge, none of the currently available carbonate system speciation programmes take this possibility into account.

Mathematical analysis and root bracketing
The solution of the Alk T -CO 2− 3 problem thus requires a more in-depth mathematical analysis. To start, we write out Eq. (2) with the Alk C expression for [CO 2− 3 ] (Eq. 9): Let us collect all the terms that are related to carbonate or water self-ionisation alkalinity on the left-hand side, introduce K 2 − 1 s and rewrite the equation as The value of γ is one of the main controls on the number of roots that this equation has.
1. If γ < 0, the equation has similar mathematical characteristics to the usual alkalinity-pH equation (Eq. 1). It has exactly one root which can be calculated using similar procedures as in the original SolveSAPHE. Please note though that this means that [CO 2− 3 ] < K 2 s . Since K 2 is of the order of 10 −9 mol (kg SW) −1 and s is of the order of 1, this case is only relevant for CO 2− 3 concentrations of the order of 1 nmol (kg SW) −1 and less.
3. If γ > 0, the left-hand side is not monotonous: it decreases from +∞ in [H + ] = 0 + to a minimum (see below) and then increases back to +∞ as [H + ] → +∞. The right-hand side is bounded and strictly increasing over the same interval (Munhoven, 2013). As a result, the equation has no root if the right-hand side is too low, exactly one if the two curves become tangent and two roots if the right-hand side is great enough.
To alleviate notation let us define the two parametric functions: where [H + ] is the independent variable and γ and A (alkalinity) are parameters. With these two function definitions, Eq. (12)  Case γ < 0 The first case can be handled similarly to the Alk T -CO 2 and Alk T -HCO − 3 pairs. Equation (12) always has exactly one root with γ < 0 as the equation function is monotonous and strictly decreasing with [H + ]. Upper and lower bounds for that root can be derived by solving the (quadratic) equations: for H inf and for H sup , and retaining the respective positive roots of each.
Case γ = 0 The second case might be considered to be only mathematically of importance as it only applies for one exact (and thus improbable) CO 2− 3 concentration value. For the sake of completeness, I nevertheless solve it.
As mentioned above, if γ = 0, Eq. (12) has one solution if and only if Alk T − Alk nWCinf > 2[CO 2− 3 ], and no solution otherwise. The root can be easily bracketed from below. It is sufficient to choose H inf such that where i denumbers the acid systems considered, except for the carbonate system, [ A [i] ] is the total amount of the acid i dissolved and K 1, [i] is the first dissociation constant of the acid system i. This equation always has a solution and, taking into account that (12), which is equivalent to L(H ; γ ) − R(H ; Alk T ) = 0 thus has a single root between H inf and H sup .
The third case is the most commonly encountered and the most challenging. With γ > 0, L([H + ]; γ ) has a minimum and the location of that minimum is a critical parameter in the analysis of this case. Let us denote the location of that minimum by H min and the value that L takes there by L min : There are two ranges of Alk T values where firm conclusions can be drawn right away.  The limiting Alk T value for which the two curves are tangential and the corresponding H tan value can be calculated with a common algorithm to characterise a bracketed local minimum, such as Brent's algorithm (Brent, 1973  . The determination of H tan is costly, generally more costly than the subsequent resolution of the pH equation itself. As mentioned right at the beginning of this section, there are extended ranges of Alk T values for which the exact knowledge of H tan is not indispensable. In these situations H min may be a sub-optimal but nevertheless sufficient separation limit for the roots (or equal to the double root itself) and cheap to calculate. If available, H tan can be used as an upper bound for the lower and as a lower bound for the greater of the two roots. To start the minimisation algorithm to derive H tan , we can use the three characteristic [H + ] values from Fig. 3 as initial conditions. These are H min together with the abscissae H L and H R of the intersection points between L([H + ]; γ ) and the horizontal line at Alk min − Alk nWCinf , which are the roots of Once H tan is known, the root brackets can be completed by the intersection points between L([H + ]; γ ) and the horizontal line at Alk T − Alk nWCinf -corresponding to the P UL and P UR points in Fig. 2 with the grey band shifted down to include the minimum -i.e. by solving the same quadratic equation than for H L and H R , with Alk min replaced by Alk T . We have again Alk T − Alk nWCinf > L min , and the equation has two positive roots. With these brackets on the two roots, any safeguarded iterative procedure, such as those implemented in SolveSAPHE, can be used to find the two roots in a controlled way.

Two roots: which one to choose?
Since every physical aqueous sample has a pH, the case without roots is essentially a theoretical one: it can actually arise only if the adopted alkalinity composition is not appropriate or if measurement errors are large. The case where an Alk T -CO 2− 3 data pair is compatible with two different pH values is, on the contrary, the most common one. SolveSAPHE-r2 has been designed as a universal pH solver and as such returns all the roots that a data pair may offer, since there is no universal criterion to decide which one of the two roots is preferable over the other.
However, at the end of the calculations, one of the two has to be chosen. Additional information, qualitative or quantita-tive, is required to make that decision. This could be a third measurable, but often even qualitative information about, say, the expected pH or the C T range might be sufficient. For typical seawater samples, the greater of the two [H + ] solutions will typically be the adequate one, following the "use the larger one" advice of Zeebe and Wolf-Gladrow (2001).
In the analysis of the Alk T -CO 2− 3 problem above, we determined the Alk T ranges that would, respectively, lead to two, one or no roots, for a given γ , i.e. [CO 2− 3 ]. To better understand the reasons why and when there are two, one or no roots and what other implications the individual roots have, it is instructive to perform the analysis the other way around: figure out how [CO 2− 3 ] evolves as a function of pH for a given value of Alk T , and determine the [CO 2− 3 ] ranges that would, respectively, lead to two, one or no roots. Such an analysis is presented in Fig. 4. For that figure, we have reconsidered the sample composition previously used in Fig. 1c and f. We thus start with Alk T fixed at 2.3 mmol kg −1 and draw the evolution of [CO 2− 3 ] as a function of pH following obtained by first using Eq. (9) to express [CO 2− 3 ] as a function of Alk C and [H + ], and then Eq.
(2) to calculate Alk C as a function of Alk T and [H + ] (and the total concentrations of all the other contributing acid-base systems, which we assume to be known, as stated initially). CO 2 and HCO − 3 evolution curves can be derived similarly, by using Eqs. (7) and (8), respectively, instead of Eq. (9). The concentration evolution for the other Alk T contributors can be calculated from their species fraction equations (see, e.g. Munhoven, 2013). Figure 4a shows the concentration curves for all the species contributing to total alkalinity and dissolved inorganic carbon, for pH values ranging from 3 to 12, and for Alk T = 2.3 mmol kg −1 ; Fig. 4b shows the [CO 2− 3 ] evolution curves for different Alk T values, ranging from 0.5 to 2.5 mmol kg −1 . Solving the Alk T -CO 2− 3 problem for our showcase sample where [CO 2− 3 ] = 0.1 mmol kg −1 thus means drawing a horizontal line through the 0.1 mmol kg −1 concentration level and locating the intersection points with the CO 2− 3 curve, if any. There are actually two of them, located at pH = 8.03 and pH = 11.43. With increasing target values for [CO 2− 3 ], i.e. moving the horizontal line upward, the two pH roots will move closer and closer together, until the maximum of the CO 3 curve is touched, at CO 3 max = 841.7 µmol kg −1 . At this exact value, there will only be one root: pH = 10.1943. Positioning the line higher up does not allow any intersection with the CO 3 curve any more: there are no roots for [CO 2− 3 ] > CO 3 max . As illustrated in Fig. 4b, the value of CO 3 max grows as Alk T increases, thus extending the range of [CO 2− 3 ] that allows for roots. Another noteworthy fact in Fig. 4b is that all the displayed CO 2− 3 curves increase to a maximum before declining and reducing to zero at some finite pH. While not all possible curves have that shape (it is possible to show that curves for Alk T < Alk nWCinf + 2 K 2 s are monotonously decreasing), all of them nevertheless go to zero. The equation CO 3 ([H + ]; Alk T ) = 0 actually always has exactly one root, for any given Alk T , since this simply requires that the numerator on the right-hand side of Eq. 17) is 0, i.e. that [H + ] is the solution of a standard alkalinity-pH equation where C T = 0. Such an equation always has exactly one positive solution, for any physically meaningful set of total concentrations of the different acid-base systems at play and any Alk T value (Munhoven, 2013). The pH value at that zero-crossing point is furthermore the maximum possible pH for the given Alk T : beyond that value, the ever-growing [OH + ] would inevitably make Alk T increase above the fixed value, thus requiring Alk C to become negative, which is not possible.
As can be seen in Fig. 4a, the high-pH solution goes together with C T [CO 2− 3 ]: CO 2− 3 represents only 4.6 % of the C T at pH = 8.03, typical for seawater, but 99.2 % at pH = 11.43. Accordingly C T = 0.1008 mmol kg −1 at the high-pH root, which is unrealistically low for many natural samples. In the marine realm, this observation regarding the high-pH root is actually correct in general. CO 2− 3 represents more than 80 % of DIC for pH > 10, and more than 90 % for pH > 10.3, as can be calculated from Eq. (6). In Fig. 4b, one can see that the maxima of the CO 3 curves are greater than 0.47 mmol kg −1 for Alk T ≥ 1.5 mmol kg −1 and that they are located at pH > 10. Since the larger of the two solutions is always at a pH greater than or equal to the maximum of the curve that it must intersect, we may conclude that for [CO 2− 3 ] < 0.47 mmol kg −1 and Alk T ≥ 1.5 mmol kg −1 , the greater of the two pH roots always implies that CO 2− 3 represents more than 80% of C T . Accordingly, even a rough estimate of one of the other relevant parameters of the carbonate system might be sufficient to reject one of the two roots.

Initialisation: rationale
Since we have bracketing intervals for each diagnosed root, we may always use the fall-back initial value H 0 = H inf H sup . This value is, however, often far from optimal. The efficient initialisation strategy of Munhoven (2013) can be generalised and adapted to each of the three pairs. For each case, we choose the most complex Alk T approximation that leads to a cubic equation. If the cubic polynomial behind that equation does not have a local minimum and a local maximum, we use the fall-back value. If such a local minimum and maximum exist, we use the quadratic Taylor expansion around the relevant extremum -this will normally be the maximum if the coefficient of the cubic term is negative, and the minimum if that coefficient is positive. If that quadratic does not have any positive roots, the fall-back initial value is used. The roots for that quadratic are then determined. For problems that have only one positive [H + ] solution (Alk T -CO 2 , Alk T -HCO − 3 and Alk T -CO 2− 3 with γ < 0), we consider that root of the quadratic expansion that is greater than the greatest location of the two extrema: if that root is lower than H inf , we use H 0 = H inf ; if it is greater than H sup , we set H 0 = H inf . For problems that have two positive [H + ] solutions (Alk T -CO 2− 3 with γ > 0 and sufficiently great Alk T ), the initial value for determining the greater of the two [H + ] solutions can be chosen exactly the same way; the initial value required to calculate the lower of the two [H + ] solutions may be more tricky. If the location of the right-hand side extremum is too close to 0, the estimated root of the cubic may be negative. In this case, the quadratic fitted to lefthand extremum should be considered as well and the greater of its roots tested. Because of the symmetries of a cubic, that root can be calculated with a few extra additions only.
The developments for each of the three input pairs are presented in full detail in the "Mathematical and Technical Details" report in the Supplement.

Reference Fortran 90 implementation
The SolveSAPHE Fortran 90 library from Munhoven (2013) -hereafter SolveSAPHE v1 -has been revised, cleaned up and upgraded to allow the processing of the additional three pairs. For the purpose of this paper, only the two main solvers have been kept: these are solve_at_general, which uses a Newton-Raphson method, and solve_at_general_sec, which uses the secant method. Both can be still be used with the same application programming interface (API) as in v1. The instances in SolveSAPHE-r2 are, however, only wrappers to the newly added Newton-Raphson-based solve_at_general2 and secant (or more precisely regula falsi) based solve_at_general2_sec both of which are able to process problems that have two roots. They return the number of roots of the problem, as well as their actual values, if any.
In the course of the developments related to the Alk T -CO 2 pair, the Newton-Raphson-based algorithm showed a few weaknesses. With the Alk T -C T pair that SolveSAPHE v1 had been designed for, each non-water alkalinity term was bounded, just like its derivative. Once CO 2 takes the role of C T these favourable properties are lost: with [CO 2 ] fixed, the carbonate alkalinity term and its derivative with respect to [H + ] become unbounded. Newton iterates can then change by large amounts and floating point over-and underflow errors on the exponential correction became common. The rate of change for Newton-Raphson iterates during each step was therefore limited to a factor of 100. With high CO 2 concentration values prescribed, there was another loss of control on the iteration sequence that had not been encountered before. At some iterations, most often at the first one, it happened that one of the two root brackets, say the upper one, was reduced to the iteration value. In the next iteration, that same bound was exceeded by the trial Newton-Raphson iterate, which was then rejected and replaced by a bisection iterate on the interval delimited by the previous iterate and the upper bracket. Since both were identical, the bisection actually produced no variation and falsely led to convergence diagnosis. This has been fixed by changing the interval whereon the bisection step is performed to that delimited by the lower and the upper brackets of the root, which are always different. 2 The unbounded variations of the carbonate alkalinity term when one of the individual species was used instead of C T furthermore required to modify the stopping criterion for the iterations: in SolveSAPHE v1, iterations are stopped as soon as the relative difference between successive iterates falls below a set tolerance ( = 10 −8 by default). However, itera-2 Both corrections have been backported to the version 1 branch of SolveSAPHE and are included in v1.1 in the SolveSAPHE archive on Zenodo (Munhoven, 2013 tions for Alk T -CO 2 and for the greater root of Alk T -CO 2− 3 were prone to early termination with that stopping criterion, as iterates only slowly changed due to the extreme gradients in the Alk C term of the equation. The stopping criterion is therefore now based upon the width of the bracketing interval, and iterations are stopped as soon as (H max − H min ) < 1 2 (H max + H min ), where H max and H min are, respectively, the upper and lower brackets of the root, which are continuously updated as iterations progress. As a consequence of this change, the number of bisection steps considerably increased. In order to speed up convergence, most bisection steps were replaced by regula falsi steps on [H min , H max ]. Bisection steps are only used occasionally when either the minimum or maximum root bracket gets updated too often in a row (three times by default) which indicates that the equation values at H max and H min have strongly different magnitudes. Unfortunately, the number of iterations required for the original SolveSAPHE pair Alk T -C T increases with this stopping criterion, without any appreciable gain in precision (compare, e.g. the number of iterations from Fig. 3b and the residuals from Fig. 1d in Munhoven (2013), with the number of iterations required here as reported in Fig. 6 for SW3 and the synthetic overview of the equation residuals reported in Tables S4 and S5 in the "Additional Results" section in the Supplement). For modelling purposes, where Alk T -C T is generally the relevant pair of data, SolveSAPHE v1 remains the most efficient choice. Tests have shown that the two safeguarded algorithms from SolveSAPHE v1 typically require 40 %-45 % less computing time than their SolveSAPHE-r2 counterparts.
Finally, as explained above, some Alk T -CO 2− 3 combinations require the solution of an auxiliary minimisation problem. For this purpose, Brent's algorithm was implemented into SolveSAPHE (translated to Fortran 90 from the Algol 60 version in Brent (1973, Sect. 5.8), taking into account the author's errata reported on https://maths-people.anu.edu. au/~brent/pub/pub011.html (last access: 24 June 2021) and his modifications to the original algorithm as implemented in https://www.netlib.org/go/fmin.f (last access: 24 June 2021).

Test case definitions
Results from the three test cases (SW1, SW2 and SW3) from Munhoven (2013) were used as starting points to define sets of Alk T -CO 2 , Alk T -HCO − 3 and Alk T -CO 2− 3 concentration pairs to drive the test case experiments. Two supplementary cases were added here: BW4 for surface brackish water with S = 3.5 and ABW5 (based upon the data of Yao and Millero (1995) for the Framvaren Fjord, Norway) for anoxic brackish water with Alk T and C T values 7 to 9 times higher than in the open ocean, as well as comparatively high alkalinity contributions from phosphates, silicates, sulfides, phosphates and ammonium. For CO 2 and CO 2− 3 , which are most conve-  10 −5 10 −3 10 −6 10 −3 10 −14 10 −2 10 −9 10 −3 10 −5 10 −3 niently handled on a logarithmic concentration scale, the representative ranges were adapted so that the range endpoints are integer powers of 10. The adopted ranges and scales are reported in Table 1. Each of the SW1, SW2 and SW3 test cases is complemented with three sets of temperature, salinity and pressure conditions for typical environments (surface cold, surface warm and deep cold seawater); for BW4, only one such set for cold surface dilute or brackish water is used (T = 275.15 K, S = 3.5) and for ABW5 one set for subsurface brackish water (P = 13.5 bar, S = 22.82).
For the comparison of the computational requirements for the processing of each set of samples, the adopted  Table 1 have been defined on the basis of their respective distributions calculated from the Alk T -C T results and although the adopted grids have the same dimensions, they do not cover exactly the same "samples" in any given test case. To overcome that inconsistency, each test experiment for the intercomparison is first carried out with the Alk T -C T pair and the results stored. For the other three pairs, the pH distribution obtained with the Alk T -C T pair for the corresponding set of temperature, salinity and pressure is first read in and the corresponding [CO 2 ],  or [CO 2− 3 ] distribution calculated on the underlying C T -Alk T grid. The soobtained arrays of species concentrations are then used to define the set of Alk T -CO 2 , Alk T -HCO − 3 and Alk T -CO 2− 3 data pairs for the benchmark calculations. This way, the test case experiments for the four different characteristic carbonate system concentrations cover exactly the same set of samples in each test case. These sample sets cannot be represented on rectangular [CO 2 ]-Alk T , [HCO − 3 ]-Alk T or [CO 2− 3 ]-Alk T grids, respectively, which is nevertheless irrelevant for the histogram syntheses presented in Figs. 6 and 7. These variants of the test cases are denoted SW1 CT , SW2 CT , SW3 CT , BW4 CT and ABW5 CT .

Results
While all the test cases have their specific relevance, we are going to focus on SW2 for most of our discussion here. SW2 covers currently observed seawater samples, thus encom-passing SW1, and conditions expected to occur over the next 50 000 years as derived from simulation experiments carried out with MBM-MEDUSA (Munhoven, 2009). A wider selection of results also for the other cases is presented in the "Additional Results" in the Supplement. pH distributions for the SW2 test case are shown in Fig. 5.
The difficulties posed by Alk T -CO 2 that were at the origin of most of the amendments to the solver algorithms show up in the histograms for the number of iterations required to reach convergence shown in Fig. 6 for solve_at_general which uses the hybrid Newton-Raphson-regula falsi-bisection scheme and in Fig. 7 for solve_at_general_sec which uses the hybrid secantregula falsi-bisection scheme. With each one of the two solvers, Alk T -CO 2 problems require in general more iterations to conclude than the other three pairs. This is especially pronounced with solve_at_general (Fig. 6), where a considerable fraction of the Alk T -CO 2 samples requires 45 to 55 and more iterations. In comparison, Alk T -C T samples typically require about four to eight iterations for naturally occurring compositions, and only in some rare instances more than 20 for the extreme SW3. The other pairs range between these two, Alk T -HCO − 3 coming closest to Alk T -C T . ABW5 shows a few deviations from the other tests cases. Here, solving the Alk T -CO 2 problem with solve_at_general nearly always takes more than 50 iterations, with solve_at_general almost always nine. The solution of Alk T -C T for ABW5 with solve_at_general_sec takes considerably more iterations than Alk T -CO 2− 3 (the fastest), and Alk T -CO 2 (the second fastest).
Finally, Figs. 6 and 7 demonstrate the superiority of solve_at_general_sec over solve_at_general. All in all, the former requires only one-fourth to one-half of the number of iterations than the latter, and it produces root approximations characterised by equation residuals that are up to 7 orders of magnitude lower than those obtained with the former (see Tables S4 and S5 in the "Additional Results" section in the Supplement). ABW5 again presents an exception to this general pattern: solve_at_general_sec re- quires typically about twice as many iterations to solve the Alk T -C T problem than solve_at_general.
All these observations are also reflected in the execution times of the two solvers. The Newton-Raphson-based solver takes more than 5 times as much time for the SW2 test case with Alk T -CO 2 than with Alk T -C T ; for Alk T -CO 2− 3 , it takes 4 times as much (for both roots though, including the solution of the minimisation problem for part of the domain). For Alk T -HCO − 3 , the difference is only 20 %. With the secant-based method, the picture is completely different: Alk T -CO 2 takes only about 30 % more time than Alk T -C T , Alk T -CO 2− 3 twice as much, whereas Alk T -HCO − 3 executes even about 5 % faster. For the Alk T -CO 2 pair of input data, the difference between the two solvers is greatest: the secantbased one takes less than one-fourth of the time taken by the Newton-Raphson-based one.
Another key factor that influences the execution times is the initialisation scheme, although the comparisons are not as clear cut as in Munhoven (2013). Safe initialisation with the geometric mean of the root brackets (the fall-back initialisation value mentioned in Sect. 2.5) results in 40 %-60 % increases of the execution times for the Alk T -C T and the Alk T -HCO − 3 input pairs, compared to the standard cubic polyno-mial one. Similar increases are obtained with a constant uniform pH = 8 initialisation. For Alk T -CO 2 , and Alk T -CO 2− 3 , the differences are much smaller and range between a decrease or an increase of up to 5 %. With these two, the quality of the root brackets seems to be more critical than the initial value.
In the analysis in Sect. 2.4.1, two characteristic thresholds for Alk T have been made out for γ > 0: an upper one at L min + Alk nWC (H min ), above which the problem always has two [H + ] solutions, and a lower one at L min +Alk nWCinf , below which the problem does not have any solution at all. For intermediate values of Alk T , it is necessary to determine H tan and Alk tan to find out how many roots the problem has, and, in case there are two, where the separation between them lies. The minimisation procedure required to determine H tan is computationally expensive as can be seen in Fig. 8 (for SW2-sc). The most probable number of iterations is in all experiments between 21 and 25; the median number is each time 0.9 ± 0.5 higher than the most probable number, due to the skew-symmetric nature of the distribution of the number of iterates, as illustrated in the insert in Fig. 8 (see also Fig. S23 in the "Additional Results" section in the Supplement). The subsequent computation of the roots is much Figure 6. Number of iterations to convergence required by the various data pairs (separately for the lower and the greater [H + ] roots of the Alk T -CO 2− 3 problem), for each of the four test cases, carried out with solve_at_general (using a hybrid Newton-Raphson-regula falsi-bisection method). The absolute maximum numbers of iterations were 58, 67, 64 and 56, for SW2, SW3, BW4 and ABW5, respectively, and 58 for SW1 (not shown).   Brent's algorithm in the SW2 test case to solve the auxiliary minimisation problem whose solution determines the number of roots of the Alk T -CO 2− 3 pair and also provides the separation between the two roots. The white area covers the region where the solution of the minimisation problems is not required as Alk T is sufficiently high so that there were two roots. The insert shows the frequency distribution of the number of iterations required. The black line in the lower right corner traces the limit between regions with two roots and without roots (compare with Fig. 5c and d).
cheaper: for the lower root, the secant-based algorithm most probably takes five iterations, and only occasionally 15-16, and for the greater root, most probably four and only rarely more than nine. The total number of samples in the SW2 test case is 1.95 million. Overall, 10 500 (0.54 %) of these do not have any root for the Alk T -CO 2− 3 pair and the solution of the minimisation problem is required for 173 445 samples (8.89 %), because H tan is required to separate the two roots. The lower threshold essentially turns out as useless: it ranges at about −28 mmol kg −1 . This is due to the hydrogen sulfate acid system which strongly dominates the Alk nWC minimum in seawater, because of the high total sulfate concentration in seawater (S T 28 mmol kg −1 ). For carbonate ion concentrations below 400 µmol kg −1 , i.e. for most of the naturally occurring waters, the Alk T -CO 2− 3 problem will always have two roots and the solution of the auxiliary minimisation problem is not required to characterise them.

Conclusions
The approach adopted in SolveSAPHE (Munhoven, 2013) to safely determine carbonate speciation in particular, and speciation calculations of mixtures of acids in aqueous solution in general, knowing only the total concentrations of the different acid systems and the total alkalinity of the system was adapted and extended here to use [CO 2 ], [HCO − 3 ] and [CO 2− 3 ] instead of the total inorganic carbon concentration, C T . The rationale can be entirely transposed to these three pairs: (1) the amended alkalinity-pH equations for Alk T -CO 2 and for Alk T -HCO − 3 still have one and only one pos-itive solution, while Alk T -CO 2− 3 may have no solution, or one or two; (2) intrinsic brackets that only depend on a priori available information can be derived for the root of the Alk T -CO 2 and Alk T -HCO − 3 problems, as well as for the two roots of Alk T -CO 2− 3 problems that may have to be solved for naturally occurring sample compositions. More uncommon but physically realistic Alk T -CO 2− 3 problems may additionally require the solution of an auxiliary minimisation problem to determine the threshold Alk T value below which the problem does not have any roots and above which it has two of them. The solution of this problem also provides a separation value of the two roots. To our best knowledge, SolveSAPHE is the first package to offer a complete solution of the Alk T -CO 2− 3 problem, autonomous above all.
The two safeguarded numerical solvers from SolveSAPHE v1 have been adapted to allow for the solution of problems that may have up to two roots. The Newton-Raphson-bisection-based solver required extensive modifications for the reliable solution of the numerically far more challenging Alk T -CO 2 , Alk T -HCO − 3 , and Alk T -CO 2− 3 problems. Most bisection steps have been replaced by regula falsi steps for increased convergence speed. The secant-bisection solver only required minimal adaptations. A Fortran 90 reference implementation, SolveSAPHE-r2, was prepared and used to evaluate the performance of the different methods for solving four benchmark problems. While the secant-bisection method was already slightly superior to the Newton-Raphson-bisection method in SolveSAPHE v1, that advantage has now become overwhelming: in SolveSAPHE-r2, it typically requires 2 to 4 times fewer iterations, and for the newly handled pairs, the equation residuals are orders of magnitude lower than the Newton-Raphson-regula falsi-bisection-based solver (typically of the order of 10 −19 -10 −18 compared to 10 −13 -10 −12 ).
For carbonate speciation problems posed by Alk T and either one of [CO 2 ],  or [CO 2− 3 ], the secant-based routine from SolveSAPHE-r2, solve_at_general2_sec, is thus clearly the method of choice; for calculations on the basis of Alk T -C T , both solve_at_general and solve_at_general_sec from SolveSAPHE v1 will perform better, although the secant-based solver is marginally faster once again. for me to reconsider SolveSAPHE and to complete it in order to combine Alk T with the individual carbonate system species instead of C T -thank you. Thanks are also due to Jean-Pierre Gattuso for pointing out issues with the colour schemes adopted in the originally submitted version of the manuscript and for guidance about how to address them. The comments and recommendations by the two anonymous referees were most helpful to improve the framing of the paper and the scope of the results. The analysis based upon Fig. 4 was suggested by Anonymous Referee #1 and recommended for inclusion in the manuscript by Anonymous Referee #2. Review statement. This paper was edited by Sandra Arndt and reviewed by two anonymous referees.