SolveSAPHE-r2 (v2.0.1): revisiting and extending the Solver Suite for Alkalinity-PH Equations for usage with CO₂, HCO₃⁻ or CO₃²⁻ input data

The successful and efficient 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₂, Alk_T–${\mathrm{HCO}}_{\mathrm{3}}^{-}$, and Alk_T–${\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}$. The mathematical properties of the three modified alkalinity–pH equations are explored. It is shown that the Alk_T–CO₂, and Alk_T–${\mathrm{HCO}}_{\mathrm{3}}^{-}$ problems have one and only one positive root for any physically sensible pair of data (i.e. such that [CO₂]>0 and $\left[{\mathrm{HCO}}_{\mathrm{3}}^{-}\right]>\mathrm{0}$). The space of Alk_T–${\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}$ 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 Alk_T–C_T pair as they exhibit strong gradients and are not always monotonous. The two main algorithms used in SolveSAPHE v1 have been revised in depth to reliably process the three additional data input pairs. The Alk_T–CO₂ pair is numerically the most challenging. With the Newton–Raphson-based solver, it takes about 5 times as long to solve as the companion Alk_T–C_T pair; the Alk_T–${\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}$ pair requires on average about 4 times as much time as the Alk_T–C_T pair. All in all, the secant-based solver offers the best performance. It outperforms the Newton–Raphson-based one by up to a factor of 4 in terms of average numbers of iterations and execution time and yet reaches equation residuals that are up to 7 orders of magnitude lower. Just like the pH solvers from the v1 series, SolveSAPHE-r2 includes automatic root bracketing and efficient initialisation schemes for the iterative solvers. For Alk_T–${\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}$ data pairs, it also determines the number of roots and calculates non-overlapping bracketing intervals. An open-source reference implementation of the new algorithms in Fortran 90 is made publicly available for usage under the GNU Lesser General Public Licence version 3 (LGPLv3) or later.

1 Introduction

Among all the aspects of the ongoing global environmental changes (climate change, ocean acidification, etc.), the solution chemistry of carbon dioxide (CO₂) 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₂] (or equivalently the partial pressure of CO₂, pCO₂, or its fugacity, fCO₂), [${\mathrm{HCO}}_{\mathrm{3}}^{-}$], [${\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}$] 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₂ and ${\mathrm{HCO}}_{\mathrm{3}}^{-}$ on one hand, and between ${\mathrm{HCO}}_{\mathrm{3}}^{-}$ and ${\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}$ 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), $C_{\text{T}}=\left[{\mathrm{CO}}_{\mathrm{2}}\right]+\left[{\mathrm{HCO}}_{\mathrm{3}}^{-}\right]+\left[{\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}\right]$; (2) total alkalinity, Alk_T; (3) pH and (4) pCO₂ or fCO₂. Recently, a procedure to measure $\left[{\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}\right]$ has been developed as well, thus increasing the number of measurables to five (Byrne and Yao, 2008; Patsavas et al., 2015; Sharp and Byrne, 2019). This latter has, however, not yet been widely adopted. With these two additional constraints, the concentrations of all the individual species as well as C_T and Alk_T can then be calculated.

There are 10 different data pairs that can be composed from the set of five independent measurable variables of the carbonate system; there are 15 if we further include ${\mathrm{HCO}}_{\mathrm{3}}^{-}$ as a sixth independent variable, although currently not (yet) measurable. Most modellers will call upon C_T and Alk_T which, besides being measurable, are also conservative and thus convenient for a budgeting approach. Experimentalists will use the pair that best suits their analytical equipment and expertise. Depending on the study requirements, not all pairs are equally attractive though. The analysis of Sharp and Byrne (2019) reveals that it is always advisable to measure pCO₂ directly if that variable is required: the uncertainty of calculated pCO₂ is always 5 to 10 times as large as that of the directly measured one. [${\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}$] can be calculated with lower uncertainty from Alk_T–C_T data than it can be directly measured; calculating it from other data pairs always bears greater uncertainty than directly measuring it. The most peculiar combination of uncertainties affects the results derived from paired measurements of pCO₂ and ${\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}$: they allow to calculate pH with the same uncertainty as if directly measured, thus providing nearly optimal values for the three individually measurable species. The uncertainties of the Alk_T–C_T calculated from it are, however, about 20 times as large as those directly measured. The currently most attractive pairs are Alk_T–pCO₂ and C_T–pCO₂, both of which allow to calculate pH with better and [${\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}$] with only slightly larger but still acceptable uncertainty than the direct measurement would offer. Direct [${\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}$] measurements, which might be most advisable for tracing carbonate mineral saturation states, are best paired with Alk_T or C_T (Sharp and Byrne, 2019). It can nevertheless be expected that, once it becomes more widely used, the measurement uncertainty currently affecting that still young measurable can be reduced and eventually become better than that of [${\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}$] calculated from Alk_T–C_T, which is currently the best option (Sharp and Byrne, 2019).

Overall, 11 out of these 15 possible pairs of independent parameters of the carbonate system can be directly solved or require at most the resolution of a quadratic equation. The remaining four pairs require iterative procedures. Besides the Alk_T–C_T pair which was addressed in full detail by Munhoven (2013) these are (1) Alk_T–CO₂, (2) Alk_T–${\mathrm{HCO}}_{\mathrm{3}}^{-}$ and (3) Alk_T–${\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}$. Such calculations are performed to an advanced level of detail with dedicated and highly specialised packages. The review of Orr et al. (2015) offers a systematic analysis of subsisting uncertainties and inconsistencies between 10 such packages, focusing on the sets of 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₂, Alk_T–${\mathrm{HCO}}_{\mathrm{3}}^{-}$, and Alk_T–${\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}$ 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.

2 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.

2.1 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:

$$\begin{array}{}\text{(1)}& R_{\text{T}}\left(\right[{\mathrm{H}}^{+}\left]\right)\equiv {\text{Alk}}_{\text{nW}}\left(\right[{\mathrm{H}}^{+}\left]\right)+{\displaystyle \frac{K_{\text{W}}}{\left[{\mathrm{H}}^{+}\right]}}-{\displaystyle \frac{\left[{\mathrm{H}}^{+}\right]}{s}}-{\text{Alk}}_{\text{T}}=\mathrm{0},\end{array}$$

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: ${\text{Alk}}_{\text{nW}}\left(\right[{\mathrm{H}}^{+}\left]\right)={\sum }_i{\text{Alk}}_{{\text{A}}_{\left[i\right]}}\left(\right[{\mathrm{H}}^{+}\left]\right)$, 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₂, ${\mathrm{HCO}}_{\mathrm{3}}^{-}$ and ${\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}$ and their sum, C_T) are all directly related to Alk_C: ${\text{Alk}}_{\text{nW}}\left(\right[{\mathrm{H}}^{+}\left]\right)={\text{Alk}}_{\text{C}}\left(\right[{\mathrm{H}}^{+}\left]\right)+{\text{Alk}}_{\text{nWC}}\left(\right[{\mathrm{H}}^{+}\left]\right)$. 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

$$\begin{array}{}\text{(2)}& \begin{array}{rl}{R}_{\text{T}}\left(\right[{\mathrm{H}}^{+}\left]\right)& \equiv {\text{Alk}}_{\text{C}}\left(\right[{\mathrm{H}}^{+}\left]\right)+{\text{Alk}}_{\text{nWC}}\left(\right[{\mathrm{H}}^{+}\left]\right)\\ & +{\displaystyle \frac{K_{\text{W}}}{\left[{\mathrm{H}}^{+}\right]}}-{\displaystyle \frac{\left[{\mathrm{H}}^{+}\right]}{s}}-{\text{Alk}}_{\text{T}}=\mathrm{0}.\end{array}\end{array}$$

The carbonate alkalinity term writes, as a function of C_T

$$\begin{array}{}\text{(3)}& {\text{Alk}}_{\text{C}}\left(\right[{\mathrm{H}}^{+}\left]\right)={\displaystyle \frac{K_{\mathrm{1}}\left[{\mathrm{H}}^{+}\right]+\mathrm{2}{K}_{\mathrm{1}}{K}_{\mathrm{2}}}{[{\mathrm{H}}^{+}{]}^{\mathrm{2}}+K_{\mathrm{1}}[{\mathrm{H}}^{+}]+K_{\mathrm{1}}{K}_{\mathrm{2}}}}{C}_{\text{T}}\end{array}$$

where K₁ and K₂ 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⁺]:

$$\begin{array}{}\text{(4)}& {\displaystyle }\left[{\mathrm{CO}}_{\mathrm{2}}\right]& {\displaystyle }={\displaystyle \frac{[{\mathrm{H}}^{+}{]}^{\mathrm{2}}}{[{\mathrm{H}}^{+}{]}^{\mathrm{2}}+K_{\mathrm{1}}[{\mathrm{H}}^{+}]+K_{\mathrm{1}}{K}_{\mathrm{2}}}}{C}_{\text{T}}\text{(5)}& {\displaystyle }\left[{\mathrm{HCO}}_{\mathrm{3}}^{-}\right]& {\displaystyle }={\displaystyle \frac{K_{\mathrm{1}}\left[{\mathrm{H}}^{+}\right]}{[{\mathrm{H}}^{+}{]}^{\mathrm{2}}+K_{\mathrm{1}}[{\mathrm{H}}^{+}]+K_{\mathrm{1}}{K}_{\mathrm{2}}}}{C}_{\text{T}}\text{(6)}& {\displaystyle }\left[{\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}\right]& {\displaystyle }={\displaystyle \frac{K_{\mathrm{1}}{K}_{\mathrm{2}}}{[{\mathrm{H}}^{+}{]}^{\mathrm{2}}+K_{\mathrm{1}}[{\mathrm{H}}^{+}]+K_{\mathrm{1}}{K}_{\mathrm{2}}}}{C}_{\text{T}}.\end{array}$$

Accordingly, Alk_C([H⁺]) may be rewritten in one of the following forms:

$$\begin{array}{}\text{(7)}& {\displaystyle }{\text{Alk}}_{\text{C}}\left(\right[{\mathrm{H}}^{+}\left]\right)& {\displaystyle }={\displaystyle \frac{K_{\mathrm{1}}\left[{\mathrm{H}}^{+}\right]+\mathrm{2}{K}_{\mathrm{1}}{K}_{\mathrm{2}}}{[{\mathrm{H}}^{+}{]}^{\mathrm{2}}}}\left[{\mathrm{CO}}_{\mathrm{2}}\right]\text{(8)}& {\displaystyle }{\text{Alk}}_{\text{C}}\left(\right[{\mathrm{H}}^{+}\left]\right)& {\displaystyle }={\displaystyle \frac{K_{\mathrm{1}}\left[{\mathrm{H}}^{+}\right]+\mathrm{2}{K}_{\mathrm{1}}{K}_{\mathrm{2}}}{K_{\mathrm{1}}\left[{\mathrm{H}}^{+}\right]}}\left[{\mathrm{HCO}}_{\mathrm{3}}^{-}\right]\text{(9)}& {\displaystyle }{\text{Alk}}_{\text{C}}\left(\right[{\mathrm{H}}^{+}\left]\right)& {\displaystyle }={\displaystyle \frac{K_{\mathrm{1}}\left[{\mathrm{H}}^{+}\right]+\mathrm{2}{K}_{\mathrm{1}}{K}_{\mathrm{2}}}{K_{\mathrm{1}}{K}_{\mathrm{2}}}}\left[{\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}\right].\end{array}$$

Figure 1 (a) pH; (b) ${\mathrm{HCO}}_{\mathrm{3}}^{-}$ and (c) ${\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}$ concentration isolines in C_T–Alk_T space; (d) CO₂, (e) ${\mathrm{HCO}}_{\mathrm{3}}^{-}$ and (f) ${\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}$ concentration isolines in pH–Alk_T space. Blank areas in panels (d), (e) and (f) represent the pH–Alk_T combinations that lead to negative Alk_C. The blue and black stars in panels (c) and (f) locate the two possible C_T and pH roots for a sample with Alk_T=2.3 mmol kg⁻¹ and $\left[{\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}\right]=\mathrm{0.1}$mmol kg⁻¹ (one of the dashed isolines, as indicated by the open star symbol in the colour scale). Figure 3 in Deffeyes (1965) is similar to panel (c).

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–2021). The three panels in the upper row of Fig. 1 show the pH and the ${\mathrm{HCO}}_{\mathrm{3}}^{-}$ and ${\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}$ concentration distributions for a reduced SW3 test case from Munhoven (2013), with the C_T range extending from 0 to 4 mmol kg⁻¹ and the Alk_T range from −1 to 3 mmol kg⁻¹ only. The pH distribution from Fig. 1a is then used to derive the corresponding CO₂ (not shown), ${\mathrm{HCO}}_{\mathrm{3}}^{-}$ and ${\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}$ concentration distributions (Fig. 1b and c). Since we intend to solve the alkalinity–pH equation for given Alk_T, and either of [CO₂], $\left[{\mathrm{HCO}}_{\mathrm{3}}^{-}\right]$ or $\left[{\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}\right]$, 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₂], $\left[{\mathrm{HCO}}_{\mathrm{3}}^{-}\right]$ and $\left[{\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}\right]$, respectively. Blank areas represent the pH–Alk_T combinations that lead to negative Alk_C.

The V- or U-shaped isolines for ${\mathrm{HCO}}_{\mathrm{3}}^{-}$ on the C_T–Alk_T graph and for ${\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}$ on the C_T–Alk_T and on the pH–Alk_T graphs show that the C_T–${\mathrm{HCO}}_{\mathrm{3}}^{-}$ and the Alk_T–${\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}$ pairs will not always provide unambiguous results. This is illustrated by the blue and black stars in Fig. 1c and f: they both lie on the 100 µmol kg⁻¹ $\left[{\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}\right]$ isoline and on the horizontal line drawn through Alk_T=2.3 mmol kg⁻¹. For that pair of data values, there are thus two compatible C_T and correspondingly two possible pH values. On the other hand, with the same Alk_T, a $\left[{\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}\right]$ of 1 mmol kg⁻¹ would not provide any solution as the 1 mmol kg⁻¹ isoline has its minimum at Alk_T=2.63 mmol kg⁻¹. Similarly, there are pairs of C_T and $\left[{\mathrm{HCO}}_{\mathrm{3}}^{-}\right]$ values that are compatible with two Alk_T values and thus two pH values, and others with none: a vertical line drawn through C_T=2.2 mmol kg⁻¹ crosses the 2.0 mmol kg⁻¹ isoline for $\left[{\mathrm{HCO}}_{\mathrm{3}}^{-}\right]$ twice and so that pair of data values leads to two pH solutions; a vertical line drawn through C_T=2.05 mmol kg⁻¹ does not cross that 2.0 mmol kg⁻¹ isoline for $\left[{\mathrm{HCO}}_{\mathrm{3}}^{-}\right]$ at all, and that pair of data values does not have any pH solution.

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₂ and Alk_T–${\mathrm{HCO}}_{\mathrm{3}}^{-}$ pairs. According to the outcome of our preliminary analysis above, the Alk_T–${\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}$ 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–${\mathrm{HCO}}_{\mathrm{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.2 The Alk_T–CO₂ problem

2.2.1 Mathematical analysis

The Alk_T–CO₂ 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). Equation (2) then becomes

$$\begin{array}{}\text{(10)}& \begin{array}{rl}& \left({\displaystyle \frac{K_{\mathrm{1}}}{\left[{\mathrm{H}}^{+}\right]}}+{\displaystyle \frac{\mathrm{2}{K}_{\mathrm{1}}{K}_{\mathrm{2}}}{[{\mathrm{H}}^{+}{]}^{\mathrm{2}}}}\right)\left[{\mathrm{CO}}_{\mathrm{2}}\right]+{\text{Alk}}_{\text{nWC}}\left(\right[{\mathrm{H}}^{+}\left]\right)\\ & +{\displaystyle \frac{K_{\text{W}}}{\left[{\mathrm{H}}^{+}\right]}}-{\displaystyle \frac{\left[{\mathrm{H}}^{+}\right]}{s}}-{\text{Alk}}_{\text{T}}=\mathrm{0}.\end{array}\end{array}$$

Just like the Alk_C([H⁺]) expression from Eq. (3) is monotonously decreasing with [H⁺] for C_T fixed, that from Eq. (7) is monotonously decreasing with [H⁺] for [CO₂] fixed. The expression on the left-hand side of Eq. (10) decreases from +∞ to −∞ for [CO₂]>0 as [H⁺] varies from 0⁺ to +∞. Equation (10) thus always has exactly one positive solution.

2.2.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

$$\left({\displaystyle \frac{K_{\mathrm{1}}}{H_{\text{inf}}}}+{\displaystyle \frac{\mathrm{2}{K}_{\mathrm{1}}{K}_{\mathrm{2}}}{H_{\text{inf}}^{\mathrm{2}}}}\right)\left[{\mathrm{CO}}_{\mathrm{2}}\right]+{\displaystyle \frac{K_{\text{W}}}{H_{\text{inf}}}}-{\displaystyle \frac{H_{\text{inf}}}{s}}={\text{Alk}}_{\text{T}}-{\text{Alk}}_{\text{nWCinf}},$$

i.e. as the positive root of the cubic equation:

$$\begin{array}{rl}& {\displaystyle \frac{H^{\mathrm{3}}}{s}}+({\text{Alk}}_{\text{T}}-{\text{Alk}}_{\text{nWCinf}})H^{\mathrm{2}}-\left(K_{\mathrm{1}}\right[{\mathrm{CO}}_{\mathrm{2}}]+K_{\text{W}})H\\ & -\mathrm{2}{K}_{\mathrm{1}}{K}_{\mathrm{2}}\left[{\mathrm{CO}}_{\mathrm{2}}\right]=\mathrm{0}.\end{array}$$

Let us denote this cubic by P(H). It is important to notice that $P\left(\mathrm{0}\right)=-\mathrm{2}{K}_{\mathrm{1}}{K}_{\mathrm{2}}\left[{\mathrm{CO}}_{\mathrm{2}}\right]<\mathrm{0}$ and $P^{\prime }\left(\mathrm{0}\right)=-\left(K_{\mathrm{1}}\right[{\mathrm{CO}}_{\mathrm{2}}]+K_{\text{W}})<\mathrm{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

$$\begin{array}{rl}& \left({\displaystyle \frac{K_{\mathrm{1}}}{H_{\text{sup}}}}+{\displaystyle \frac{\mathrm{2}{K}_{\mathrm{1}}{K}_{\mathrm{2}}}{H_{\text{sup}}^{\mathrm{2}}}}\right)\left[{\mathrm{CO}}_{\mathrm{2}}\right]+{\displaystyle \frac{K_{\text{W}}}{H_{\text{sup}}}}-{\displaystyle \frac{H_{\text{sup}}}{s}}\\ & ={\text{Alk}}_{\text{T}}-{\text{Alk}}_{\text{nWCsup}},\end{array}$$

i.e. as the positive root of the cubic equation

$$\begin{array}{rl}& {\displaystyle \frac{H^{\mathrm{3}}}{s}}+({\text{Alk}}_{\text{T}}-{\text{Alk}}_{\text{nWCsup}})H^{\mathrm{2}}-\left(K_{\mathrm{1}}\right[{\mathrm{CO}}_{\mathrm{2}}]+K_{\text{W}})H\\ & -\mathrm{2}{K}_{\mathrm{1}}{K}_{\mathrm{2}}\left[{\mathrm{CO}}_{\mathrm{2}}\right]=\mathrm{0},\end{array}$$

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.2.2):

1. locate the local minimum of the cubic, in H_min>0;

2. develop the cubic as a quadratic Taylor expansion, Q(H), around that minimum;

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\left(H\right)=Q\left(H\right)+(H-H_{\text{min}}{)}^{\mathrm{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).

2.3 The Alk_T–HCO${}_{\mathrm{3}}^{-}$ problem

2.3.1 Mathematical analysis

For the Alk_T–${\mathrm{HCO}}_{\mathrm{3}}^{-}$ pair, the Alk_C([H⁺]) term in Eq. (2) is written as in Eq. (8):

$$\begin{array}{}\text{(11)}& \begin{array}{rl}& \left(\mathrm{1}+{\displaystyle \frac{\mathrm{2}{K}_{\mathrm{2}}}{\left[{\mathrm{H}}^{+}\right]}}\right)\left[{\mathrm{HCO}}_{\mathrm{3}}^{-}\right]+{\text{Alk}}_{\text{nWC}}\left(\right[{\mathrm{H}}^{+}\left]\right)\\ & +{\displaystyle \frac{K_{\text{W}}}{\left[{\mathrm{H}}^{+}\right]}}-{\displaystyle \frac{\left[{\mathrm{H}}^{+}\right]}{s}}-{\text{Alk}}_{\text{T}}=\mathrm{0}.\end{array}\end{array}$$

The expression on the left-hand side of Eq. (11) decreases monotonously from +∞ to −∞ for $\left[{\mathrm{HCO}}_{\mathrm{3}}^{-}\right]>\mathrm{0}$ fixed as [H⁺] varies from 0⁺ to +∞. Equation (11) thus always has exactly one positive solution.

2.3.2 Root bracketing

The lower bound H_inf can be chosen such that

$$\left(\mathrm{1}+{\displaystyle \frac{\mathrm{2}{K}_{\mathrm{2}}}{H_{\text{inf}}}}\right)\left[{\mathrm{HCO}}_{\mathrm{3}}^{-}\right]+{\displaystyle \frac{K_{\text{W}}}{H_{\text{inf}}}}-{\displaystyle \frac{H_{\text{inf}}}{s}}={\text{Alk}}_{\text{T}}-{\text{Alk}}_{\text{nWCinf}},$$

i.e. as the positive root of the quadratic equation

$$\begin{array}{rl}& {\displaystyle \frac{H^{\mathrm{2}}}{s}}+({\text{Alk}}_{\text{T}}-{\text{Alk}}_{\text{nWCinf}}-[{\mathrm{HCO}}_{\mathrm{3}}^{-}\left]\right)H\\ & -\left(\mathrm{2}{K}_{\mathrm{2}}\right[{\mathrm{HCO}}_{\mathrm{3}}^{-}]+K_{\text{W}})=\mathrm{0}.\end{array}$$

Similarly, the upper bound H_sup can be chosen such that

$$\left(\mathrm{1}+{\displaystyle \frac{\mathrm{2}{K}_{\mathrm{2}}}{H_{\text{sup}}}}\right)\left[{\mathrm{HCO}}_{\mathrm{3}}^{+}\right]+{\displaystyle \frac{K_{\text{W}}}{H_{\text{sup}}}}-{\displaystyle \frac{H_{\text{sup}}}{s}}={\text{Alk}}_{\text{T}}-{\text{Alk}}_{\text{nWCsup}},$$

i.e. as the positive root of the quadratic equation

$$\begin{array}{rl}& {\displaystyle \frac{H^{\mathrm{2}}}{s}}+({\text{Alk}}_{\text{T}}-{\text{Alk}}_{\text{nWCsup}}-[{\mathrm{HCO}}_{\mathrm{3}}^{-}\left]\right)H\\ & -\left(\mathrm{2}{K}_{\mathrm{2}}\right[{\mathrm{HCO}}_{\mathrm{3}}^{-}]+K_{\text{W}})=\mathrm{0}.\end{array}$$

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.

2.4 The Alk_T–CO${}_{\mathrm{3}}^{\mathrm{2}-}$ problem

Whereas any physically meaningful Alk_T–[CO₂] or Alk_T–[${\mathrm{HCO}}_{\mathrm{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–[${\mathrm{CO}}_{\mathrm{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–[${\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}$] pairs. This little-known fact was already documented in the 1960s (see, e.g. Deffeyes, 1965)¹. On the other hand, there are also Alk_T–[${\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}$] 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.

2.4.1 Mathematical analysis and root bracketing

The solution of the Alk_T–${\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}$ problem thus requires a more in-depth mathematical analysis. To start, we write out Eq. (2) with the Alk_C expression for [${\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}$] (Eq. 9):

$$\begin{array}{rl}& {\displaystyle \frac{K_{\mathrm{1}}\left[{\mathrm{H}}^{+}\right]+\mathrm{2}{K}_{\mathrm{1}}{K}_{\mathrm{2}}}{K_{\mathrm{1}}{K}_{\mathrm{2}}}}\left[{\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}\right]+{\text{Alk}}_{\text{nWC}}\left(\right[{\mathrm{H}}^{+}\left]\right)\\ & +{\displaystyle \frac{K_{\text{W}}}{\left[{\mathrm{H}}^{+}\right]}}-{\displaystyle \frac{\left[{\mathrm{H}}^{+}\right]}{s}}-{\text{Alk}}_{\text{T}}=\mathrm{0}.\end{array}$$

Let us collect all the terms that are related to carbonate or water self-ionisation alkalinity on the left-hand side, introduce the shorthand $\mathit{\gamma }=\frac{\left[{\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}\right]}{K_{\mathrm{2}}}-\frac{\mathrm{1}}{s}$ and rewrite the equation as

$$\begin{array}{}\text{(12)}& \mathit{\gamma }\left[{\mathrm{H}}^{+}\right]+{\displaystyle \frac{K_{\text{W}}}{\left[{\mathrm{H}}^{+}\right]}}+\mathrm{2}\left[{\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}\right]={\text{Alk}}_{\text{T}}-{\text{Alk}}_{\text{nWC}}\left(\right[{\mathrm{H}}^{+}\left]\right).\end{array}$$

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 $\left[{\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}\right]<\frac{K_{\mathrm{2}}}{s}$. Since K₂ is of the order of 10⁻⁹ mol (kg SW)⁻¹ and s is of the order of 1, this case is only relevant for ${\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}$ concentrations of the order of 1 nmol (kg SW)⁻¹ and less.

2. If γ=0 (i.e. if $\left[{\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}\right]=\frac{K_{\mathrm{2}}}{s}$), the equation has exactly one root if ${\text{Alk}}_{\text{T}}-\mathrm{2}\left[{\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}\right]-{\text{Alk}}_{\text{nWCinf}}>\mathrm{0}$, and no root otherwise.

3. If γ>0, the left-hand side is not monotonous: it decreases from +∞ in $\left[{\mathrm{H}}^{+}\right]={\mathrm{0}}^{+}$ to a minimum (see below) and then increases back to +∞ as $\left[{\mathrm{H}}^{+}\right]\to +\mathrm{\infty }$. 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.

Figure 2 Schematic representation of the general characteristics of the $L\left(\right[{\mathrm{H}}^{+}];\mathit{\gamma })$ and $R\left(\right[{\mathrm{H}}^{+}];A)$ components of the alkalinity–pH equation for the ${\text{Alk}}_{\text{T}}-{\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}$ pair. The grey band delimits the (monotonous) variations of $R\left(\right[{\mathrm{H}}^{+}];A)$ for a given alkalinity A. The band moves up and down without being distorted as A is increased, respectively, decreased. For a given pair of Alk_T–${\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}$ concentrations, the actual equation to solve is $L\left(\right[{\mathrm{H}}^{+}];\mathit{\gamma })=R\left(\right[{\mathrm{H}}^{+}];{\text{Alk}}_{\text{T}})$, where $\mathit{\gamma }=\frac{\left[{\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}\right]}{K_{\mathrm{2}}}-\frac{\mathrm{1}}{s}$. γ=0 thus corresponds to $\left[{\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}\right]=\frac{K_{\mathrm{2}}}{s}$.

To alleviate notation let us define the two parametric functions:

$$\begin{array}{}\text{(13)}& {\displaystyle }L\left(\right[{\mathrm{H}}^{+}];\mathit{\gamma })& {\displaystyle }=\mathit{\gamma }\left[{\mathrm{H}}^{+}\right]+{\displaystyle \frac{K_{\text{W}}}{\left[{\mathrm{H}}^{+}\right]}}+\mathrm{2}\left[{\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}\right]\text{(14)}& {\displaystyle }R\left(\right[{\mathrm{H}}^{+}];A)& {\displaystyle }=A-{\text{Alk}}_{\text{nWC}}\left(\right[{\mathrm{H}}^{+}\left]\right),\end{array}$$

where [H⁺] is the independent variable and γ and A (alkalinity) are parameters. With these two function definitions, Eq. (12) is then rewritten as $L\left(\right[{\mathrm{H}}^{+}];\mathit{\gamma })=R\left(\right[{\mathrm{H}}^{+}];{\text{Alk}}_{\text{T}})$. Schematic representations of the three γ cases and of the L and R functions are shown in Fig. 2.

Case γ<0

The first case can be handled similarly to the Alk_T–CO₂ and Alk_T–${\mathrm{HCO}}_{\mathrm{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:

$$\begin{array}{}\text{(15)}& \mathit{\gamma }{H}_{\text{inf}}+{\displaystyle \frac{K_{\text{W}}}{H_{\text{inf}}}}+\mathrm{2}\left[{\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}\right]={\text{Alk}}_{\text{T}}-{\text{Alk}}_{\text{nWCinf}}\end{array}$$

for H_inf and

$$\begin{array}{}\text{(16)}& \mathit{\gamma }{H}_{\text{sup}}+{\displaystyle \frac{K_{\text{W}}}{H_{\text{sup}}}}+\mathrm{2}\left[{\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}\right]={\text{Alk}}_{\text{T}}-{\text{Alk}}_{\text{nWCsup}}\end{array}$$

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) ${\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}$ 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 ${\text{Alk}}_{\text{T}}-{\text{Alk}}_{\text{nWCinf}}>\mathrm{2}\left[{\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}\right]$, and no solution otherwise. The root can be easily bracketed from below. It is sufficient to choose H_inf such that

$${\displaystyle \frac{K_{\text{W}}}{H_{\text{inf}}}}={\text{Alk}}_{\text{T}}-\mathrm{2}\left[{\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}\right]-{\text{Alk}}_{\text{nWCinf}}$$

leading to $L(H_{\text{inf}};\mathit{\gamma })-R(H_{\text{inf}};{\text{Alk}}_{\text{T}})>\mathrm{0}$. The analogue equation for H_sup, with Alk_nWCinf replaced by Alk_nWCsup (cf. Eqs. 15 and 16) does not work if ${\text{Alk}}_{\text{T}}-{\text{Alk}}_{\text{nWCsup}}\le \mathrm{2}\left[{\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}\right]$. The newly derived asymptotic approximation for Alk_nWC([H⁺]) as $\left[{\mathrm{H}}^{+}\right]\to +\mathrm{\infty }$ (see the “Mathematical and Technical Details” report in the Supplement) nevertheless provides a means to derive an upper bound. It is sufficient to choose H_sup such that

$${\displaystyle \frac{K_{\text{W}}}{H_{\text{sup}}}}={\text{Alk}}_{\text{T}}-\mathrm{2}\left[{\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}\right]-{\text{Alk}}_{\text{nWCinf}}-{\displaystyle \frac{{\sum }_i\left[\mathrm{\Sigma }{\mathrm{A}}_{\left[i\right]}\right]K_{\mathrm{1},\left[i\right]}}{H_{\text{sup}}}},$$

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

$${\mathrm{Alk}}_{\mathrm{nWC}}\left(\right[{\mathrm{H}}^{+}\left]\right)<{\mathrm{Alk}}_{\mathrm{nWCinf}}+{\displaystyle \frac{{\sum }_i\left[\mathrm{\Sigma }{\mathrm{A}}_{\left[i\right]}\right]K_{\mathrm{1},\left[i\right]}}{\left[{\mathrm{H}}^{+}\right]}},$$

which is valid for $\left[{\mathrm{H}}^{+}\right]>\mathrm{0}$, it is straightforward to show that $L(H_{\text{sup}};\mathit{\gamma })-R(H_{\text{sup}};{\text{Alk}}_{\text{T}})<\mathrm{0}$ with this choice. Equation (12), which is equivalent to $L(H;\mathit{\gamma })-R(H;{\text{Alk}}_{\text{T}})=\mathrm{0}$ thus has a 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\left(\right[{\mathrm{H}}^{+}];\mathit{\gamma })$ 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:

$$H_{\text{min}}=\sqrt{{\displaystyle \frac{K_{\text{W}}}{\mathit{\gamma }}}}\text{and}{L}_{\text{min}}=\mathrm{2}\sqrt{\mathit{\gamma }{K}_{\text{W}}}+\mathrm{2}\left[{\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}\right].$$

There are two ranges of Alk_T values where firm conclusions can be drawn right away.

1. If $R(H_{\text{min}};{\text{Alk}}_{\text{T}})>L_{\text{min}}$, i.e. if ${\text{Alk}}_{\text{T}}>L_{\text{min}}+{\text{Alk}}_{\text{nWC}}\left(H_{\text{min}}\right)$, Eq. (12) has two distinct roots, since R(H;Alk_T) is bounded. Furthermore, the roots – let us provisionally denote the lower one H₁ and the greater one H₂ – are such that H₁<H_min and H₂>H_min. H_min can thus be used as an upper bracket for H₁ and as a lower bracket for H₂. However, if ${\text{Alk}}_{\text{T}}-{\text{Alk}}_{\text{nWCsup}}>L_{\text{min}}$, the abscissae of the intersection points P_LL and P_LR (see Fig. 2), which are solutions of

   $$\mathit{\gamma }H+{\displaystyle \frac{K_{\text{W}}}{H}}={\text{Alk}}_{\text{T}}-\mathrm{2}\left[{\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}\right]-{\text{Alk}}_{\text{nWCsup}},$$

   provide tighter brackets than H_min.

2. If ${\text{Alk}}_{\text{T}}-{\text{Alk}}_{\text{nWCinf}}\le L_{\text{min}}$, i.e. if ${\text{Alk}}_{\text{T}}\le L_{\text{min}}+{\text{Alk}}_{\text{nWCinf}}$, Eq. (12) does not have any roots.

For intermediate values of Alk_T, no firm quantitative statement regarding the root(s) of Eq. (12) can be made a priori. As Alk_T decreases from L_min+Alk_nWC(H_min) to L_min+Alk_nWCinf, Eq. (12) will at first still have two roots, but both are greater than or equal to H_min. At some intermediate value, $L\left(\right[{\mathrm{H}}^{+}];\mathit{\gamma })$ and $R\left(\right[{\mathrm{H}}^{+}];{\text{Alk}}_{\text{T}})$ become tangent. At this point, Eq. (12) has one double root, which is the abscissa of that tangent point, H_tan. H_tan is actually a universally valid separation limit between two roots, if there are any. For lower values of Alk_T, the problem does not have any solutions.

Figure 3 Determination of the A value for which the $L\left(\right[{\mathrm{H}}^{+}];\mathit{\gamma })$ and $R\left(\right[{\mathrm{H}}^{+}];A)$ curves become tangent, or, equivalently, the lowest Alk_T value for which the equation $L\left(\right[{\mathrm{H}}^{+}];\mathit{\gamma })-R\left(\right[{\mathrm{H}}^{+}];{\text{Alk}}_{\text{T}})=\mathrm{0}$ has a solution. Panel (a) shows how relevant characteristic points can be derived by considering the particular $R\left(\right[{\mathrm{H}}^{+}];A)$ curve that intersects $L\left(\right[{\mathrm{H}}^{+}];\mathit{\gamma })$ at its minimum. Panel (b) shows the locus of the solutions of $L\left(\right[{\mathrm{H}}^{+}];\mathit{\gamma })-R\left(\right[{\mathrm{H}}^{+}];{\text{Alk}}_{\text{T}})=\mathrm{0}$ in an $\left[{\mathrm{H}}^{+}\right]-{\text{Alk}}_{\text{T}})$ graph, i.e. the curve ${\text{Alk}}_{\text{T}})=L(\left[{\mathrm{H}}^{+}\right];\mathit{\gamma })+{\text{Alk}}_{\text{ncW}}(\left[{\mathrm{H}}^{+}\right])$. Please notice that $A_{\text{min}}=L_{\text{min}}+{\text{Alk}}_{\text{ncW}}\left(H_{\text{min}}\right)$ denotes the alkalinity value obtained for $\left[{\mathrm{H}}^{+}\right]=H_{\text{min}}$, and not the minimum value of the curve shown in panel (b). See text for details.

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). To start, we reconsider $L\left(\right[{\mathrm{H}}^{+}];\mathit{\gamma })-R\left(\right[{\mathrm{H}}^{+}];A)=\mathrm{0}$ not as an equation in [H⁺] for given parameter values γ (or, equivalently, [${\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}$]) 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:

$$A\left(\right[{\mathrm{H}}^{+}\left]\right)=L\left(\right[{\mathrm{H}}^{+}];\mathit{\gamma })+{\text{Alk}}_{\text{ncW}}\left(\right[{\mathrm{H}}^{+}\left]\right).$$

Figure 3 shows how the two problems are related and which information can be derived from the analysis of $L\left(\right[{\mathrm{H}}^{+}];\mathit{\gamma })$ and $R\left(\right[{\mathrm{H}}^{+}];A)$ to contribute to the solution of the minimisation of A([H⁺]). 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\left(\right[{\mathrm{H}}^{+}];\mathit{\gamma })$ and the horizontal line at Alk_min−Alk_nWCinf, which are the roots of

$$\mathit{\gamma }{H}^{\mathrm{2}}-({\text{Alk}}_{\text{min}}-\mathrm{2}[{\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}]-{\text{Alk}}_{\text{nWCinf}})H+K_{\text{W}}=\mathrm{0}.$$

By construction, ${\text{Alk}}_{\text{min}}-{\text{Alk}}_{\text{nWCinf}}>L_{\text{min}}=\mathrm{2}\sqrt{\mathit{\gamma }{K}_{\text{W}}}+\mathrm{2}\left[{\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}\right]$. The discriminant of this quadratic equation is therefore strictly positive and the equation has two positive roots (their sum and their product are positive) as required. It is possible to show that the second derivative of $R\left(\right[{\mathrm{H}}^{+}];A)$ with respect to [H⁺] is positive provided that the successive dissociation constants K_j,[i] of the different acid systems (denumbered by i) resulting from the dissociation of an acid ${\text{H}}_{n_{\left[i\right]}}{A}_{\left[i\right]}$ are such that $K_{j,\left[i\right]}<\frac{\mathrm{1}}{\mathrm{2}}{K}_{j-\mathrm{1},\left[i\right]},j=\mathrm{2},\mathrm{\dots },n_{\left[i\right]}$ – a very weak constraint as these constants generally differ by a few orders of magnitude. This has been verified to be the case for acid systems with $n_{\left[i\right]}=\mathrm{1},\mathrm{\dots },\mathrm{12}$. The underlying technical developments can be found in the “Mathematical and Technical Details” report in the Supplement. $R\left(\right[{\mathrm{H}}^{+}];A)$ is thus concave, while $L\left(\right[{\mathrm{H}}^{+}];\mathit{\gamma })$ is convex for γ>0. A([H⁺]) thus has only a single local minimum comprised between H_L and H_R.

Once H_tan is known, the root brackets can be completed by the intersection points between $L\left(\right[{\mathrm{H}}^{+}];\mathit{\gamma })$ 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 ${\text{Alk}}_{\text{T}}-{\text{Alk}}_{\text{nWCinf}}>L_{\text{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.

Figure 4 (a) Evolution of the concentrations of the different species composing Alk_T–C_T, as a function of pH, for Alk_T=2.3 mmol kg⁻¹. C_T and the concentrations all of its components reduce to 0 at pH = 11.48 (marked by the long-dashed vertical black line) in this example. The dashed orange line represents the joint contribution of B(OH)${}_{\mathrm{4}}^{-}$ and OH⁻ which are the dominant Alk_T contributors at high pH. (b) $\left[{\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}\right]$ as a function of pH for different Alk_T values (indicated in mmol kg⁻¹ for each curve). Each curve represents a horizontal cross-section at the corresponding Alk_T level through the $\left[{\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}\right]$ distribution depicted in Fig. 1e.

2.4.2 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–${\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}$ 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 quantitative, 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–${\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}$ problem above, we determined the Alk_T ranges that would, respectively, lead to two, one or no roots, for a given γ, i.e. $\left[{\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}\right]$. 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 $\left[{\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}\right]$ evolves as a function of pH for a given value of Alk_T, and determine the $\left[{\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}\right]$ 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⁻¹ and draw the evolution of $\left[{\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}\right]$ as a function of pH following

$$\begin{array}{}\text{(17)}& {\mathrm{CO}}_{\mathrm{3}}\left(\right[{\text{H}}^{+}];{\text{Alk}}_{\text{T}})={\displaystyle \frac{{\text{Alk}}_{\text{T}}-{\text{Alk}}_{\text{nWC}}\left(\right[{\text{H}}^{+}\left]\right)-\frac{K_{\text{W}}}{\left[{\text{H}}^{+}\right]}+\frac{\left[{\text{H}}^{+}\right]}{s}}{\frac{\left[{\text{H}}^{+}\right]}{K_{\mathrm{2}}}+\mathrm{2}}},\end{array}$$

obtained by first using Eq. (9) to express $\left[{\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}\right]$ 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₂ and ${\mathrm{HCO}}_{\mathrm{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⁻¹; Fig. 4b shows the $\left[{\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}\right]$ evolution curves for different Alk_T values, ranging from 0.5 to 2.5 mmol kg⁻¹. Solving the Alk_T–${\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}$ problem for our showcase sample where $\left[{\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}\right]=\mathrm{0.1}$ mmol kg⁻¹ thus means drawing a horizontal line through the 0.1 mmol kg⁻¹ concentration level and locating the intersection points with the ${\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}$ curve, if any. There are actually two of them, located at pH = 8.03 and pH = 11.43. With increasing target values for $\left[{\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}\right]$, i.e. moving the horizontal line upward, the two pH roots will move closer and closer together, until the maximum of the CO₃ curve is touched, at CO_3 max=841.7 µmol kg⁻¹. 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₃ curve any more: there are no roots for $\left[{\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}\right]>{\text{CO}}_{\text{3\hspace{0.17em}max}}$. As illustrated in Fig. 4b, the value of CO_3 max grows as Alk_T increases, thus extending the range of [${\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}$] that allows for roots. Another noteworthy fact in Fig. 4b is that all the displayed ${\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}$ 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 ${\text{Alk}}_{\text{T}}<{\text{Alk}}_{\text{nWCinf}}+\mathrm{2}\frac{K_{\mathrm{2}}}{s}$ are monotonously decreasing), all of them nevertheless go to zero. The equation ${\mathrm{CO}}_{\mathrm{3}}\left(\right[{\text{H}}^{+}];{\text{Alk}}_{\text{T}})=\mathrm{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_{\mathrm{T}}\simeq \left[{\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}\right]$: ${\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}$ 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⁻¹ 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. ${\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}$ 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₃ curves are greater than 0.47 mmol kg⁻¹ for Alk_T≥1.5 mmol kg⁻¹ 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 $\left[{\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}\right]<\mathrm{0.47}$ mmol kg⁻¹ and Alk_T≥1.5 mmol kg⁻¹, the greater of the two pH roots always implies that ${\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}$ 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.

2.5 Initialisation: rationale

Since we have bracketing intervals for each diagnosed root, we may always use the fall-back initial value $H_{\mathrm{0}}=\sqrt{H_{\text{inf}}{H}_{\text{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₂, Alk_T–${\mathrm{HCO}}_{\mathrm{3}}^{-}$ and Alk_T–${\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}$ 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₀=H_inf; if it is greater than H_sup, we set H₀=H_inf. For problems that have two positive [H⁺] solutions (Alk_T–${\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}$ 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 left-hand 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.

3 Numerical experiments

3.1 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₂ 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₂ takes the role of C_T these favourable properties are lost: with [CO₂] 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₂ 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.² 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 ϵ ($\mathit{ϵ}={\mathrm{10}}^{-\mathrm{8}}$ by default). However, iterations for Alk_T–CO₂ and for the greater root of Alk_T–${\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}$ 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_{\text{max}}-H_{\text{min}})<\mathit{ϵ}\frac{\mathrm{1}}{\mathrm{2}}(H_{\text{max}}+H_{\text{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.

Table 1 Ranges of variation for the input variables for the five test cases. Experiments always considered Alk_T and either one of C_T, [CO₂], [${\mathrm{HCO}}_{\mathrm{3}}^{-}$] or [${\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}$].

Finally, as explained above, some Alk_T–${\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}$ 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).

3.2 Results and discussion

3.2.1 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₂, Alk_T–${\mathrm{HCO}}_{\mathrm{3}}^{-}$ and Alk_T–${\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}$ 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₂ and ${\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}$, which are most conveniently 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 [CO₂], [${\mathrm{HCO}}_{\mathrm{3}}^{-}$] or [${\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}$] distributions are adapted. Although the [CO₂], $\left[{\mathrm{HCO}}_{\mathrm{3}}^{-}\right]$ and $\left[{\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}\right]$ ranges for each test case reported on 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₂], $\left[{\mathrm{HCO}}_{\mathrm{3}}^{-}\right]$ or $\left[{\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}\right]$ distribution calculated on the underlying C_T–Alk_T grid. The so-obtained arrays of species concentrations are then used to define the set of Alk_T–CO₂, Alk_T–${\mathrm{HCO}}_{\mathrm{3}}^{-}$ and Alk_T–${\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}$ 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₂]–Alk_T, [${\mathrm{HCO}}_{\mathrm{3}}^{-}$]–Alk_T or [${\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}$]–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.

Figure 5 pH distributions for the SW2-sc test case (SW2 under cold surface conditions, where T=273.15 K, S=35 and P=0 bar), obtained with solve_at_general2_sec: (a) Alk_T–CO₂; (b) Alk_T–${\mathrm{HCO}}_{\mathrm{3}}^{-}$; (c) the lower [H⁺] root (higher pH root) of Alk_T–${\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}$; (d) the greater [H⁺] root (lower pH root) of Alk_T–${\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}$. The thick dashed grey line in panels (c) and (d) shows the critical limit above which the Alk_T–${\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}$ 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.

3.2.2 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 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). 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₂ 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 secant–regula falsi–bisection scheme. With each one of the two solvers, Alk_T–CO₂ 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₂ 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–${\mathrm{HCO}}_{\mathrm{3}}^{-}$ coming closest to Alk_T–C_T. ABW5 shows a few deviations from the other tests cases. Here, solving the Alk_T–CO₂ 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–${\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}$ (the fastest), and Alk_T–CO₂ (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 requires typically about twice as many iterations to solve the Alk_T–C_T problem than solve_at_general.

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–${\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}$ 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).

Figure 7 Number of iterations to convergence required by the various data pairs (separately for the lower and the greater [H⁺] roots of the Alk_T–${\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}$ problem), 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 20, 21, 29 and 27, for SW2, SW3, BW4 and ABW5, respectively, and 20 for SW1 (not shown).

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₂ than with Alk_T–C_T; for Alk_T–${\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}$, 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–${\mathrm{HCO}}_{\mathrm{3}}^{-}$, the difference is only 20 %. With the secant-based method, the picture is completely different: Alk_T–CO₂ takes only about 30 % more time than Alk_T–C_T, Alk_T–${\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}$ twice as much, whereas Alk_T–${\mathrm{HCO}}_{\mathrm{3}}^{-}$ executes even about 5 % faster. For the Alk_T–CO₂ 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.

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–${\mathrm{HCO}}_{\mathrm{3}}^{-}$ input pairs, compared to the standard cubic polynomial one. Similar increases are obtained with a constant uniform pH = 8 initialisation. For Alk_T–CO₂, and Alk_T–${\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}$, 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 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–${\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}$ 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⁻¹. 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_{\text{T}}\simeq \mathrm{28}\mathrm{mmol}{\mathrm{kg}}^{-\mathrm{1}}$). For carbonate ion concentrations below 400 µmol kg⁻¹, i.e. for most of the naturally occurring waters, the Alk_T–${\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}$ problem will always have two roots and the solution of the auxiliary minimisation problem is not required to characterise them.

Figure 8 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 Alk_T–${\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}$ 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).

4 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₂], [${\mathrm{HCO}}_{\mathrm{3}}^{-}$] and [${\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}$] 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₂ and for Alk_T–${\mathrm{HCO}}_{\mathrm{3}}^{-}$ still have one and only one positive solution, while Alk_T–${\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}$ 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₂ and Alk_T–${\mathrm{HCO}}_{\mathrm{3}}^{-}$ problems, as well as for the two roots of Alk_T–${\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}$ problems that may have to be solved for naturally occurring sample compositions. More uncommon but physically realistic Alk_T–${\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}$ 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–${\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}$ 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₂, Alk_T–${\mathrm{HCO}}_{\mathrm{3}}^{-}$, and Alk_T–${\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}$ 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⁻¹⁹–10⁻¹⁸ compared to 10⁻¹³–10⁻¹²).

For carbonate speciation problems posed by Alk_T and either one of [CO₂], [${\mathrm{HCO}}_{\mathrm{3}}^{-}$] or [${\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}$], 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.

Appendix A: The direct cases

For the sake of completeness, I provide here succinct “recipes” to calculate all the different carbonate system related variables, knowing two of them. Many of these were already known in the 1960s (see, e.g. Park, 1969). The Guide to Best Practices for Ocean CO₂ Measurements (Dickson et al., 2007) lists the most commonly used pairs and furthermore includes procedures for selected triplets and quartets, for which not all of the equilibrium constants are required. In the following, we assume that there are direct and invertible relationships between [CO₂] and the fugacity (fCO₂) or the partial pressure (pCO₂) of CO₂ and between pH and [H⁺] on any chosen pH scale. We therefore restrict ourselves to [CO₂] and [H⁺].

The conditions for the existence of a solution are generally that the concentrations of H⁺ and of the DIC species are strictly positive. In some instances, the input data must fulfil additional constraints that are, however, not always straightforward to quantitatively state a priori.

C_T and CO₂, C_T and CO${}_{\mathbf{3}}^{\mathbf{2}-}$

: (1) With these two pairs, the $\left[{\mathrm{CO}}_{\mathrm{2}}\right]/C_{\text{T}}$ fraction (the $\left[{\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}\right]/C_{\text{T}}$ fraction) is fixed and Eq. (4) (Eq. 6) defines a quadratic equation in [H⁺] that always allows for exactly one positive solution; (2) calculate the remaining two species concentrations from their respective species fraction; (3) Alk_T from Eq. (1). In addition to the positivity of the species concentrations, the following constraints must be met: [CO₂]<C_T and $\left[{\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}\right]<C_{\text{T}}$.

C_T and HCO${}_{\mathbf{3}}^{-}$

: With this pair, the $\left[{\mathrm{HCO}}_{\mathrm{3}}^{-}\right]/C_{\text{T}}$ fraction, denoted by b hereafter, is fixed, and Eq. (5) becomes a quadratic equation in [H⁺]. That equation has two positive solutions if $b<\mathrm{1}/(\mathrm{1}+\mathrm{2}\sqrt{K_{\mathrm{2}}/K_{\mathrm{1}}})$, one double root if $b=\mathrm{1}/(\mathrm{1}+\mathrm{2}\sqrt{K_{\mathrm{2}}/K_{\mathrm{1}}})$ and no real solutions if $b>\mathrm{1}/(\mathrm{1}+\mathrm{2}\sqrt{K_{\mathrm{2}}/K_{\mathrm{1}}})$. This is well illustrated in Fig. 1b above: there are C_T–$\left[{\mathrm{HCO}}_{\mathrm{3}}^{-}\right]$ combinations that allow for two different Alk_T values and, equivalently, two pH values, and there are others that do not allow for any. Please notice that the threshold fraction $\mathrm{1}/(\mathrm{1}+\mathrm{2}\sqrt{K_{\mathrm{2}}/K_{\mathrm{1}}})$ is always lower than 1 and the natural a priori constraint requiring that b<1 is thus insufficient to guarantee a solution: for T=275.15 K, S=35 and P=0 bar, the threshold ratio is 94.48 %.
When there are two roots, one faces a similar dilemma as with the Alk_T–${\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}$ problem: which one to choose? Most often, the lower of the two will again be the appropriate one, as that one typically leads to Alk_T>C_T, whereas the greater one leads to Alk_T<C_T. This criterion might be sufficient to discriminate between the two – in seawater it generally is – but in some instances additional information, quantitative or qualitative, might be in order.
In general, (1) solve the quadratic equation and choose the appropriate of the two roots; (2) calculate [CO₂] and $\left[{\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}\right]$ from their respective species fractions; (3) Alk_T from Eq. (1).

CO₂ and HCO${}_{\mathbf{3}}^{-}$

: (1) [H⁺] from K₁; (2) [${\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}$] from K₂; (3) C_T can be calculated from the three carbonate species concentrations; (4) Alk_T from Eq. (1).

CO₂ and CO${}_{\mathbf{3}}^{\mathbf{2}-}$

: (1) [${\mathrm{HCO}}_{\mathrm{3}}^{-}$] from $[{\mathrm{HCO}}_{\mathrm{3}}^{-}{]}^{\mathrm{2}}=K_{\mathrm{1}}/K_{\mathrm{2}}[{\mathrm{CO}}_{\mathrm{2}}\left]\right[{\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}]$; (2) C_T from the three carbonate species concentrations; (3) [H⁺] from K₁ or K₂; (4) Alk_T from Eq. (1).

HCO${}_{\mathbf{3}}^{-}$ and CO${}_{\mathbf{3}}^{\mathbf{2}-}$

: (1) calculate [H⁺] from K₂; (2) [CO₂] from K₁; (3) C_T from the three carbonate species concentrations; (4) Alk_T from Eq. (1).

CO₂ and H⁺

: (1) calculate [${\mathrm{HCO}}_{\mathrm{3}}^{-}$] from K₁; (2) calculate [${\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}$] from K₂; (3) C_T from the three carbonate species concentrations; (4) Alk_T from Eq. (1).

HCO${}_{\mathbf{3}}^{-}$ and H⁺

: (1) calculate [CO₂] from K₁; (2) calculate [${\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}$] from K₂; (3) C_T from the three carbonate species concentrations; (4) Alk_T from Eq. (1).

CO${}_{\mathbf{3}}^{\mathbf{2}-}$ and H⁺

: (1) calculate [${\mathrm{HCO}}_{\mathrm{3}}^{-}$] from K₂; (2) calculate [CO₂] from K₁; (3) C_T from the three carbonate species concentrations; (4) Alk_T from Eq. (1).

Alk_T and H⁺

: (1) C_T from Eq. (2); (2) individual species concentrations from the species fractions.
As illustrated in Fig. 1d–f above, there are Alk_T–[H⁺] combinations that lead to physically unrealistic negative Alk_C. Following Eq. (3), negative Alk_C requires negative C_T, and vice versa.
The shape of the blank area depends on the non-carbonate contributors to the total alkalinity. In practice, such incompatible combinations are unlikely to arise from measurements, except if the adopted set of Alk_T contributors is inappropriate.

C_T and H⁺

: individual species concentrations from the species fractions; Alk_T from Eq. (1).