SolveSAPHE-r 2 : 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 efficient approach at the basis of SOLVESAPHE (Munhoven, 2013), which determines the carbonate system speciation by calculating pH from total alkalinity (AlkT) and dissolved inorganic carbon (CT), and which converges from any physically sensible pair of such data, has been adapted and further developed for work with AlkT & CO2, AlkT & HCO3 and AlkT & CO 2− 3 . The mathematical properties of the three modified alkalinity-pH equations are explored. It is shown that the AlkT & CO2 and AlkT & HCO3 problems have one and only one positive root for any physically sensible pair 5 of data (i.e, such that, resp., [CO2]> 0 and [HCO3 ]> 0). The space of AlkT & 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 modified alkalinity-pH equations is far more demanding than that for the original AlkT-CT pair as they exhibit strong gradients and are not always monotonous. The two main algorithms used from SOLVESAPHE v. 1 had to be revised in depth to reliably process the three additional data input pairs. The AlkT & CO2 pair is numerically the most challenging. With the Newton-Raphson based solver, it takes about 10 five times as long to solve as the companion AlkT & CT pair, while AlkT & CO2 requires about four times as much time. All in all, it is nevertheless the secant based solver that offers the best performances. It outperforms the Newton-Raphson based one by up to a factor of four, to reach equation residuals that are up to seven orders of magnitude lower. Just like the pH solvers from routines from the v. 1 series, SOLVESAPHE v. 2 includes automatic root bracketing and efficient initialisation schemes for the iterative solvers. For AlkT & CO2− 3 pairs of data, it also determines the number of roots and calculates non-overlapping 15 bracketing intervals. An open source reference implementation in Fortran 90 of the new algorithms is made publicly available for usage under the GNU Lesser General Public Licence v. 3 or later.

i. e., eq. (21) from Munhoven (2013), where is that part of the total alkalinity not related to the water self-ionization, with i denumbering the acid systems resulting from the dissolution of acids A [i] whoe dissolution products contribute to total alkalinity. [H + ] is the proton concentration expressed 50 on one of the commonly used pH scales (total, seawater) and s is a factor to convert from that scale to the free scale. total concentrations of all the acid-base systems considered. We denote these two by Alk nWCinf and Alk nWCsup , respectively.
Eq. (1) is thus formally rewritten as The carbonate alkalinity term writes, as a function of C T 60 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 70 which will be used hereafter.

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 sulphate, phosphate, silicate, etc.) are known or can be derived from 75 adequate parametric relationships.
Eleven out of the fifteen different possible pairs of independent parameters of the carbonate system do not require any complex iterative schemes, but can be directly solved or require at most the resolution of a quadratic equation. For the sake of completeness -and with minimal details only -"recipes" for solving these straightforward cases are provided in the appendix.
Alternative approaches can be found in the literature, such as in the Guide to Best Practices for Ocean CO 2 Measurements 80 (Dickson et al., 2007). Dickson et al. (2007) also provide pathways for usage with triplets or quartets of input data. These 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.
(1) establish the analytical properties of the modified pH-alkalinity 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. 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.

Alk
The Alk T & CO 2 pair can be dealt with in a similar way to the Alk T & C T pair in the original SOLVESAPHE. The Alk C ([H + ]) term in eq. (2) is written as in eq. (7) and eq. (2) becomes 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 has also 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,  Any bracketing root-finding algorithm can then be used to solve the modified pH-alkalinity equation (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 pH-alkalinity equation and any bracketing root-finding algorithm can be used to solve it.  Fig. 1b and 1c. On one hand, there are two compatible C T , and equivalently two pH values for most Alk T -[CO − 3 ] pairs. This little-known fact was already documented in the 1960s (see, e. g., Deffeyes (1965)). On the other hand, there are also

Alk
3 ] pairs that do not allow for any solution, as they lead to negative carbonate alkalinity. To our best knowledge, none 150 of the currently used carbonate system speciation programs takes this possibility into account.
The solution of the Alk T -[CO − 3 ] problem thus requires a more in-depth preliminary mathematical analysis. To start, we write out eq. (2) with the Alk C expression for [CO − 3 ] (eq. (9)): Let us collect all the terms that are related to carbonate or water self-ionization alkalinity at the left-hand side, introduce the 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 as the usual pH-alkalinity equation (eq. (1). It has exactly one root which can be calculated using similar procedures as in the original SOLVESAPHE. Please notice though that this means that [CO 2− 3 ] < K2 s . Since K 2 is of the order of 10 −9 mol/kg-SW 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 and less.
2. If γ = 0 (i. e., if [CO 2− 3 ] = K2 s ), the equation has exactly one root if Alk T − 2[CO 2− 3 ] − Alk nWCinf > 0, no root otherwise. 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.

Mathematical analysis and root bracketing 170
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.  , for a given alkalinity A. The band moves up and down without being distorted as A is increased, resp., decreased. For a given pair of AlkT and CO 2− 3 concentrations, the actual equation Case γ < 0 The first case can be handled similarly to the Alk T & CO 2 and Alk T & HCO − 3 pairs. Eq. (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.

185
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, we 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 else. The root can be easily bracketed from below. It is sufficient to chose H inf such that The analogue equation for H sup , with Alk nWCinf replaced by Alk nWCsup (cf. eqs. (15) and (16) 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 is equivalent to L(H; γ) − R(H; Alk T ) = 0 thus has one single root between H inf and H sup .
Case γ > 0 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 There are two ranges of Alk T values where firm conclusions can be drawn right away.
1. If R(H min ; Alk T ) > L min , i. e., if Alk T > L min + Alk nWC (H min ), eq. (12)  if Alk T − Alk nWCsup > L min , the abscissae of the intersection points P LL and P LR (see Fig. 2), which are solutions of The limiting Alk T value for which the two curves are tangent and the corresponding H tan value can be calculated with a 220 common algorithm to characterize a bracketed local minimum, such as Brent's algorithm (Brent, 1973). To start, we reconsider for given parameter value γ (or, equivalently, [CO 2− 3 ]) and A, but rather as an implicit definition for A as a function of [H + ], for a given γ (here γ > 0). This implicit function definition can actually be solved explicitly here:   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 on Fig. 2 with the grey band shifted down to include the 245 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.

Initialisation: rationale
Since we have bracketing intervals for each of the root(s), we may always use the fall-back initial value H 0 = H inf H sup .

250
This value is, however, often far from optimal. The efficient initialisation strategy of Munhoven (2013) can be generalized and adapted to each of the three pairs. For each case, we chose 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 255 quadratic does not have any positive roots, the fall-back initial value is used. The roots for that quadratic are then determined.  which are able to process problems that have two roots. They return the number of roots of the problem, as well as the actual roots.
In 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 285 whereon the bisection step is performed to that delimited by the lower and the upper brackets of the root, which are always different. 1 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 v. 1 iterations are stopped as soon as the relative difference between successive iterates falls below a set tolerance ( = 10 −8 by default). However, iterations 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 290 only slowly changed due to the extreme gradients in the Alk C term of the equation function. 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 resp. the upper and lower brackets of the root, which are continuously updated as iterations progress.

Test cases
Results from the three test cases SW1, SW2 and SW3 from Munhoven (2013) were used to define ranges of CO 2 , HCO − 3 and 310 CO 2− 3 concentrations to drive the test case experiments (see Figs. S1-S4 in the Additional Results in the Supplement). For the sake of simplicity, we keep the corresponding denominations here. A supplementary 'SW4' was added for brackish water (with S = 3.5). For CO 2 and CO 2− 3 , which are most conveniently handled on a logarithmic scale, the representative ranges were adapted so that the range endpoints are integer powers of ten. The adopted ranges and scales are reported in Tab. 1. Each of the tests cases is complemented with temperature, salinity and pressure conditions for four typical environments (surface 315 cold, surface warm and deep cold seawater, as well as surface brackish water.   (c) and (d) shows the critical limit above which the AlkT & CO 2− 3 always has two roots. Below this limit further calculations are required to determine the number of solutions. More details are given in the text and in the Supplement. Please notice the different scales on the horizontal axes and for the pH colour coding in the four panels.

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 sea-water samples, thus encompassing SW1, and conditions expected to occur over the next 50,000 years as derived from simulation experiments carried out with MBM-MEDUSA (Munhoven, 2009) AlkT & CO 2− 3 ), for each of the four test cases, carried out with solve_at_general_sec (using a hybrid secant-regula falsi-bisection method). The absolute maximum numbers of iterations were respectively 22, 23, 29 and 27, for SW1 to SW4, resp. All these observations are also reflected in the execution times of the two solvers. The Newton-Raphson based solver takes 335 more than five 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 four times as much (for both roots though). 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 secant based one takes less than one fourth of the time taken by the Newton-Raphson based one.

340
Another key factor that influences the execution times is the initialisation scheme, although the comparisons are not as clear-cut as in Munhoven (2013) . Number of iterations required by Brent's algorithm in the SW2 test case to solve the auxiliary minimisation problem whose solution determines the number of roots of the AlkT & CO 2− 3 pair and also provides the separation between the two roots The white areas covers the region where the solution of the minimisation problems is not required as AlkT is sufficiently high so that there were two roots). The black line in the lower right corner traces the limit between regions with two roots and without roots -cf. Fig. 4c. initialisation. For Alk T & CO 2 and Alk T & CO 2− 3 , the differences are much smaller and range between a decrease or an 345 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 section 2.3.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 350 minimisation procedure required to determine H tan is computationally expensive as can be seen on Fig. 7. The most probable number of iterations is in all cases 25, the median number is 27-28, due to the skew-symmetric nature of the distribution of the number of iterates (not shown). The subsequent computation of the roots is much 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. 10,500 (0.54%) of these do 355 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 meq/kg. This is is due to the hydrogen sulphate acid system which strongly dominates the Alk nW minimum in seawater, because of the high total sulphate concentration in sea-water (S T 28 mmol/kg). For carbonate ion concentrations below 400 µmol/kg, i.e., for most of the naturally occurring waters, the Alk T & CO 2− 3 problem will always have two roots and 360 the solution of the auxiliary minimisation problem is not required to characterise them.
The approach adopted in SOLVESAPHE (Munhoven, 2013) to safely determine carbonate speciation in particular, and speci- (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 & 370 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 v. 1 have been adapted to allow for the solution of problems 375 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 v. 2, was prepared and used to evaluate the performances of the different methods for solving four benchmark problems. While the secant-bisection method was already 380 slightly superior to the Newton-Raphson-bisection method in SOLVESAPHE v. 1, that advantage has now become overwhelming: in SOLVESAPHE v. 2, it typically requires two to four times less 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