Development and technical paper 06 Jul 2021
Development and technical paper  06 Jul 2021
SolveSAPHEr2 (v2.0.1): revisiting and extending the Solver Suite for AlkalinityPH Equations for usage with CO_{2}, HCO_{3}^{−} or CO_{3}^{2−} input data
 Dépt. d'Astrophysique, Géophysique et Océanographie, Université de Liège, 4000 Liège, Belgium
 Dépt. d'Astrophysique, Géophysique et Océanographie, Université de Liège, 4000 Liège, Belgium
Correspondence: Guy Munhoven (guy.munhoven@uliege.be)
Hide author detailsCorrespondence: Guy Munhoven (guy.munhoven@uliege.be)
The successful and efficient approach at the basis of the Solver Suite for AlkalinityPH Equations (SolveSAPHE) (Munhoven, 2013), which determines the carbonate system speciation by calculating pH from total alkalinity (Alk_{T}) and dissolved inorganic carbon (C_{T}), and which converges for any physically sensible pair of such data, has been adapted and further developed to work with Alk_{T}–CO_{2}, Alk_{T}–${\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_{2}, 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_{2}]>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_{2} pair is numerically the most challenging. With the Newton–Raphsonbased 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 secantbased solver offers the best performance. It outperforms the Newton–Raphsonbased 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, SolveSAPHEr2 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 nonoverlapping bracketing intervals. An opensource 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.
 Article
(3598 KB)  Companion paper

Supplement
(14322 KB)  BibTeX
 EndNote
Among all the aspects of the ongoing global environmental changes (climate change, ocean acidification, etc.), the solution chemistry of carbon dioxide (CO_{2}) is one of the best known. The related chemistry of the carbonate system in the oceans and other aqueous environments is well understood and routinely monitored and modelled. The equilibria between the carbonate system species involves four variables: [CO_{2}] (or equivalently the partial pressure of CO_{2}, pCO_{2}, or its fugacity, fCO_{2}), [${\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_{2} 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_{2} or fCO_{2}. 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_{2} directly if that variable is required: the uncertainty of calculated pCO_{2} 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_{2} 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_{2} and C_{T}–pCO_{2}, 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}, (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_{2}, 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 WolfGladrow (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.
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:
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 righthand side is that part of the total alkalinity that is not related to the water selfionisation: ${\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 noncarbonate alkalinity, Alk_{nWC}([H^{+}]), since the relevant carbonate system parameters (the concentrations of CO_{2}, ${\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
The carbonate alkalinity term writes, as a function of C_{T}
where K_{1} and K_{2} are the first and second stoichiometric dissociation constants of carbonic acid. The individual carbonate species fractions of C_{T} can be expressed as a function of [H^{+}]:
Accordingly, Alk_{C}([H^{+}]) may be rewritten in one of the following forms:
In order to get a first idea about the complications that we might encounter for the solution of the three new data pairs, we start with an exploratory analysis using the Solver Suite for AlkalinityPH 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^{−1} and the Alk_{T} range from −1 to 3 mmol kg^{−1} only. The pH distribution from Fig. 1a is then used to derive the corresponding CO_{2} (not shown), ${\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_{2}], $\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_{2}], $\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 Ushaped 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^{−1} $\left[{\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}}\right]$ isoline and on the horizontal line drawn through Alk_{T}=2.3 mmol kg^{−1}. 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^{−1} would not provide any solution as the 1 mmol kg^{−1} isoline has its minimum at Alk_{T}=2.63 mmol kg^{−1}. 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^{−1} crosses the 2.0 mmol kg^{−1} 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^{−1} does not cross that 2.0 mmol kg^{−1} 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_{2} 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 indepth 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_{2} problem
2.2.1 Mathematical analysis
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). Equation (2) then becomes
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_{2}] fixed. The expression on the lefthand side of Eq. (10) decreases from +∞ to −∞ for [CO_{2}]>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
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\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
i.e. as the positive root of the cubic equation
which also has one and only one positive root for the same reasons as above.
The positive roots of these cubic equations can be found by adopting a strategy similar to that used for the cubic initialisation of the iterative solution in SolveSAPHE (Munhoven, 2013, Sect. 3.2.2):

locate the local minimum of the cubic, in H_{min}>0;

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

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 rootfinding 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):
The expression on the lefthand 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
i.e. as the positive root of the quadratic equation
Similarly, the upper bound H_{sup} can be chosen such that
i.e. as the positive root of the quadratic equation
Both equations always have two roots: one positive and one negative; their product is negative as indicated by the constant term. With the respective positive roots, we have again bounds for the solution of the modified alkalinity–pH equation and any bracketing rootfinding 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_{2}] 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 littleknown fact was already documented in the 1960s (see, e.g. Deffeyes, 1965)^{1}. 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 indepth mathematical analysis. To start, we write out Eq. (2) with the Alk_{C} expression for [${\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}}$] (Eq. 9):
Let us collect all the terms that are related to carbonate or water selfionisation alkalinity on the lefthand 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
The value of γ is one of the main controls on the number of roots that this equation has.

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_{2} is of the order of 10^{−9} mol (kg SW)^{−1} and s is of the order of 1, this case is only relevant for ${\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}}$ concentrations of the order of 1 nmol (kg SW)^{−1} and less.

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.

If γ>0, the lefthand 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 righthand side is bounded and strictly increasing over the same interval (Munhoven, 2013). As a result, the equation has no root if the righthand side is too low, exactly one if the two curves become tangent and two roots if the righthand side is great enough.
To alleviate notation let us define the two parametric functions:
where [H^{+}] is the independent variable and γ and A (alkalinity) are parameters. With these two function definitions, Eq. (12) 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_{2} 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:
for H_{inf} and
for H_{sup}, and retaining the respective positive roots of each.
Case γ=0
The second case might be considered to be only mathematically of importance as it only applies for one exact (and thus improbable) ${\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
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
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
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}:
There are two ranges of Alk_{T} values where firm conclusions can be drawn right away.

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_{1} and the greater one H_{2} – are such that H_{1}<H_{min} and H_{2}>H_{min}. H_{min} can thus be used as an upper bracket for H_{1} and as a lower bracket for H_{2}. 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}.

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.
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:
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 suboptimal 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
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.
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. SolveSAPHEr2 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 WolfGladrow (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^{−1} and draw the evolution of $\left[{\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}}\right]$ as a function of pH following
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_{2} 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^{−1}; 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^{−1}. 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^{−1} thus means drawing a horizontal line through the 0.1 mmol kg^{−1} 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_{3} curve is touched, at CO_{3 max}=841.7 µmol kg^{−1}. At this exact value, there will only be one root: pH = 10.1943. Positioning the line higher up does not allow any intersection with the CO_{3} curve any more: there are no roots for $\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 righthand 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 zerocrossing point is furthermore the maximum possible pH for the given Alk_{T}: beyond that value, the evergrowing [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 highpH 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^{−1} at the highpH root, which is unrealistically low for many natural samples. In the marine realm, this observation regarding the highpH 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_{3} curves are greater than 0.47 mmol kg^{−1} for Alk_{T}≥1.5 mmol kg^{−1} and that they are located at pH >10. Since the larger of the two solutions is always at a pH greater than or equal to the maximum of the curve that it must intersect, we may conclude that for $\left[{\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}}\right]<\mathrm{0.47}$ mmol kg^{−1} and Alk_{T}≥1.5 mmol kg^{−1}, 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 fallback initial value ${H}_{\mathrm{0}}=\sqrt{{H}_{\text{inf}}\phantom{\rule{0.125em}{0ex}}{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 fallback 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 fallback initial value is used. The roots for that quadratic are then determined. For problems that have only one positive [H^{+}] solution (Alk_{T}–CO_{2}, Alk_{T}–${\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_{0}=H_{inf}; if it is greater than H_{sup}, we set H_{0}=H_{inf}. For problems that have two positive [H^{+}] solutions (Alk_{T}–${\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 righthand side extremum is too close to 0, the estimated root of the cubic may be negative. In this case, the quadratic fitted to lefthand extremum should be considered as well and the greater of its roots tested. Because of the symmetries of a cubic, that root can be calculated with a few extra additions only.
The developments for each of the three input pairs are presented in full detail in the “Mathematical and Technical Details” report in the Supplement.
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 SolveSAPHEr2 are, however, only wrappers to the newly added Newton–Raphsonbased solve_at_general2 and secant (or more precisely regula falsi) based solve_at_general2_sec both of which are able to process problems that have two roots. They return the number of roots of the problem, as well as their actual values, if any.
In the course of the developments related to the Alk_{T}–CO_{2} pair, the Newton–Raphsonbased algorithm showed a few weaknesses. With the Alk_{T}–C_{T} pair that SolveSAPHE v1 had been designed for, each nonwater alkalinity term was bounded, just like its derivative. Once CO_{2} takes the role of C_{T} these favourable properties are lost: with [CO_{2}] fixed, the carbonate alkalinity term and its derivative with respect to [H^{+}] become unbounded. Newton iterates can then change by large amounts and floating point over and underflow errors on the exponential correction became common. The rate of change for Newton–Raphson iterates during each step was therefore limited to a factor of 100. With high CO_{2} concentration values prescribed, there was another loss of control on the iteration sequence that had not been encountered before. At some iterations, most often at the first one, it happened that one of the two root brackets, say the upper one, was reduced to the iteration value. In the next iteration, that same bound was exceeded by the trial Newton–Raphson iterate, which was then rejected and replaced by a bisection iterate on the interval delimited by the previous iterate and the upper bracket. Since both were identical, the bisection actually produced no variation and falsely led to convergence diagnosis. This has been fixed by changing the interval whereon the bisection step is performed to that delimited by the lower and the upper brackets of the root, which are always different.^{2} The unbounded variations of the carbonate alkalinity term when one of the individual species was used instead of C_{T} furthermore required to modify the stopping criterion for the iterations: in SolveSAPHE v1, iterations are stopped as soon as the relative difference between successive iterates falls below a set tolerance ϵ ($\mathit{\u03f5}={\mathrm{10}}^{\mathrm{8}}$ by default). However, iterations for Alk_{T}–CO_{2} 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{\u03f5}\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 SolveSAPHEr2 counterparts.
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://mathspeople.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_{2}, 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_{2} 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_{2}], [${\mathrm{HCO}}_{\mathrm{3}}^{}$] or [${\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}}$] distributions are adapted. Although the [CO_{2}], $\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_{2}], $\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 soobtained arrays of species concentrations are then used to define the set of Alk_{T}–CO_{2}, 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_{2}]–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}.
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 MBMMEDUSA (Munhoven, 2009). A wider selection of results also for the other cases is presented in the “Additional Results” in the Supplement. pH distributions for the SW2 test case are shown in Fig. 5.
The difficulties posed by Alk_{T}–CO_{2} that were at the origin of most of the amendments to the solver algorithms show up in the histograms for the number of iterations required to reach convergence shown in Fig. 6 for solve_at_general
which uses the hybrid Newton–Raphson–regula falsi–bisection scheme and in Fig. 7 for solve_at_general_sec
which uses the hybrid secant–regula falsi–bisection scheme.
With each one of the two solvers, Alk_{T}–CO_{2} problems require in general more iterations to conclude than the other three pairs.
This is especially pronounced with solve_at_general
(Fig. 6), where a considerable fraction of the Alk_{T}–CO_{2} samples requires 45 to 55 and more iterations.
In comparison, Alk_{T}–C_{T} samples typically require about four to eight iterations for naturally occurring compositions, and only in some rare instances more than 20 for the extreme SW3.
The other pairs range between these two, Alk_{T}–${\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_{2} problem with solve_at_general
nearly always takes more than 50 iterations, with solve_at_general
almost always nine.
The solution of Alk_{T}–C_{T} for ABW5 with solve_at_general_sec
takes considerably more iterations than Alk_{T}–${\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}}$ (the fastest), and Alk_{T}–CO_{2} (the second fastest).
Finally, Figs. 6 and 7 demonstrate the superiority of solve_at_general_sec
over solve_at_general
.
All in all, the former requires only onefourth to onehalf 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
.
All these observations are also reflected in the execution times of the two solvers. The Newton–Raphsonbased solver takes more than 5 times as much time for the SW2 test case with Alk_{T}–CO_{2} than with Alk_{T}–C_{T}; for Alk_{T}–${\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 secantbased method, the picture is completely different: Alk_{T}–CO_{2} 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_{2} pair of input data, the difference between the two solvers is greatest: the secantbased one takes less than onefourth of the time taken by the Newton–Raphsonbased 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 fallback 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_{2}, 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 SW2sc). 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 skewsymmetric 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 secantbased 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^{−1}. This is due to the hydrogen sulfate acid system which strongly dominates the Alk_{nWC} minimum in seawater, because of the high total sulfate concentration in seawater (${S}_{\text{T}}\simeq \mathrm{28}\phantom{\rule{0.125em}{0ex}}\mathrm{mmol}\phantom{\rule{0.125em}{0ex}}{\mathrm{kg}}^{\mathrm{1}}$). For carbonate ion concentrations below 400 µmol kg^{−1}, 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.
The approach adopted in SolveSAPHE (Munhoven, 2013) to safely determine carbonate speciation in particular, and speciation calculations of mixtures of acids in aqueous solution in general, knowing only the total concentrations of the different acid systems and the total alkalinity of the system was adapted and extended here to use [CO_{2}], [${\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_{2} 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_{2} 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–bisectionbased solver required extensive modifications for the reliable solution of the numerically far more challenging Alk_{T}–CO_{2}, 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, SolveSAPHEr2, 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 SolveSAPHEr2, 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–bisectionbased solver (typically of the order of 10^{−19}–10^{−18} compared to 10^{−13}–10^{−12}).
For carbonate speciation problems posed by Alk_{T} and either one of [CO_{2}], [${\mathrm{HCO}}_{\mathrm{3}}^{}$] or [${\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}}$], the secantbased routine from SolveSAPHEr2, 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 secantbased solver is marginally faster once again.
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_{2} 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_{2}] and the fugacity (fCO_{2}) or the partial pressure (pCO_{2}) of CO_{2} and between pH and [H^{+}] on any chosen pH scale. We therefore restrict ourselves to [CO_{2}] 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_{2}, 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_{2}]<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_{2}] and $\left[{\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}}\right]$ from their respective species fractions; (3) Alk_{T} from Eq. (1).  CO_{2} and HCO${}_{\mathbf{3}}^{}$

: (1) [H^{+}] from K_{1}; (2) [${\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}}$] from K_{2}; (3) C_{T} can be calculated from the three carbonate species concentrations; (4) Alk_{T} from Eq. (1).
 CO_{2} 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]\phantom{\rule{0.125em}{0ex}}\right[{\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}}]$; (2) C_{T} from the three carbonate species concentrations; (3) [H^{+}] from K_{1} or K_{2}; (4) Alk_{T} from Eq. (1).
 HCO${}_{\mathbf{3}}^{}$ and CO${}_{\mathbf{3}}^{\mathbf{2}}$

: (1) calculate [H^{+}] from K_{2}; (2) [CO_{2}] from K_{1}; (3) C_{T} from the three carbonate species concentrations; (4) Alk_{T} from Eq. (1).
 CO_{2} and H^{+}

: (1) calculate [${\mathrm{HCO}}_{\mathrm{3}}^{}$] from K_{1}; (2) calculate [${\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}}$] from K_{2}; (3) C_{T} from the three carbonate species concentrations; (4) Alk_{T} from Eq. (1).
 HCO${}_{\mathbf{3}}^{}$ and H^{+}

: (1) calculate [CO_{2}] from K_{1}; (2) calculate [${\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}}$] from K_{2}; (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}; (2) calculate [CO_{2}] from K_{1}; (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 noncarbonate 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).
All the Fortran 90 codes of SolveSAPHE version 1 series (of which v1.0.3 was used to derive the results presented in Fig. 1) are available on Zenodo from Munhoven (2013–2021) for use under the GNU Lesser General Public Licence version 3 (LGPLv3) or later. The codes for SolveSAPHEr2 (v2.0.1) that are described in this paper are included in the Supplement and made available for use under the same licence. They are also archived on Zenodo (Munhoven, 2021). Future bugfix releases and updates will also be archived there.
Epitalon et al. (2021) have ported SolveSAPHEr2 to R (not used here) for usage under the GNU General Public License version 2 (GPL2) or 3 (GPL3).
The supplement related to this article is available online at: https://doi.org/10.5194/gmd1442252021supplement.
The author declares that there is no conflict of interest.
Publisher’s note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
James Orr provided the kickoff momentum for me to reconsider SolveSAPHE and to complete it in order to combine Alk_{T} with the individual carbonate system species instead of C_{T} – thank you. Thanks are also due to JeanPierre Gattuso for pointing out issues with the colour schemes adopted in the originally submitted version of the manuscript and for guidance about how to address them. The comments and recommendations by the two anonymous referees were most helpful to improve the framing of the paper and the scope of the results. The analysis based upon Fig. 4 was suggested by Anonymous Referee #1 and recommended for inclusion in the manuscript by Anonymous Referee #2.
Financial support for this work was provided by the Belgian Fund for Scientific Research – F.R.S.FNRS (project SERENATA, grant CDR J.0123.19). Guy Munhoven is a research associate with the Belgian Fund for Scientific Research – F.R.S.FNRS.
This paper was edited by Sandra Arndt and reviewed by two anonymous referees.
Brent, R. P.: Algorithms for minimization without derivatives, PrenticeHall, Englewood Cliffs, NJ, 1973. a, b
Byrne, R. H. and Yao, W.: Procedures for measurement of carbonate ion concentrations in seawater by direct spectrophotometric observations of Pb(II) complexation, Mar. Chem., 112, 128–135, https://doi.org/10.1016/j.marchem.2008.07.009, 2008. a
Deffeyes, K. S.: Carbonate Equilibria : A Graphic and Algebraic Approach, Limnol. Oceanogr., 10, 412–426, https://doi.org/10.4319/lo.1965.10.3.0412, 1965. a, b
Dickson, A. G., Sabine, C. L., and Christian, J. R. (Eds.): Guide to Best Practices for Ocean CO_{2} Measurements, vol. 3, in: PICES Special Publication, Carbon Dioxide Information and Analysis Center, Oak Ridge (TN), available at: https://cdiac.essdive.lbl.gov/ftp/oceans/Handbook_2007/Guide_all_in_one.pdf (last access: 24 June 2021), 2007. a, b, c, d
Epitalon, J.M., Gattuso, J.P., and Munhoven, G.: SolveSAPHE: Solver Suite for AlkalinityPH Equations, available at: https://CRAN.Rproject.org/package=SolveSAPHE, last access: 24 June 2021. a
Munhoven, G.: Future CCD and CSH variations: Deepsea impact of ocean acidification, Geochim. Cosmochim. Ac., 73, A917, 2009. a
Munhoven, G.: Mathematics of the total alkalinity–pH equation – pathway to robust and universal solution algorithms: the SolveSAPHE package v1.0.1, Geosci. Model Dev., 6, 1367–1388, https://doi.org/10.5194/gmd613672013, 2013. a, b, c, d, e, f, g, h, i, j, k, l, m, n, o, p, q
Munhoven, G.: SolveSAPHE (Solver Suite for AlkalinityPH Equations) [software], Zenodo, https://doi.org/10.5281/zenodo.3752250, 2013–2021. a, b, c
Munhoven, G.: SolveSAPHEr2 (Solver Suite for AlkalinityPH Equations–Release 2) [software], Zenodo, https://doi.org/10.5281/zenodo.4771132, 2021. a
Orr, J. C., Epitalon, J.M., and Gattuso, J.P.: Comparison of ten packages that compute ocean carbonate chemistry, Biogeosciences, 12, 1483–1510, https://doi.org/10.5194/bg1214832015, 2015. a
Park, P. K.: Oceanic CO_{2} System: An Evaluation of ten Methods of Investigation, Limnol. Oceanogr., 14, 179–186, https://doi.org/10.4319/lo.1969.14.2.0179, 1969. a
Patsavas, M. C., Byrne, R. H., Yang, B., Easley, R. A., Wanninkhof, R., and Liu, X.: Procedures for direct spectrophotometric determination of carbonate ion concentrations : Measurements in US Gulf of Mexico and East Coast waters, Mar. Chem., 168, 80–85, https://doi.org/10.1016/j.marchem.2014.10.015, 2015. a
Sharp, J. D. and Byrne, R. H.: Carbonate ion concentrations in seawater: Spectrophotometric determination at ambient temperatures and evaluation of propagated calculation uncertainties, Mar. Chem., 209, 70–80, https://doi.org/10.1016/j.marchem.2018.12.001, 2019. a, b, c, d
Yao, W. and Millero, F. J.: The Chemistry Of the Anoxic Waters in the Framvaren Fjord, Norway, Aquat. Geochem., 1, 53–88, https://doi.org/10.1007/BF01025231, 1995. a
Zeebe, R. E. and WolfGladrow, D.: CO_{2} in seawater : Equilibrium, kinetics, isotopes, vol. 65, in: Elsevier Oceanography Series, Elsevier, Amsterdam (NL), ISBN 9780444505798, 2001. a, b, c
Zeebe and WolfGladrow (2001) appear to be aware of it. In their recipe for given Alk_{T} and [${\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}}$] (on pp. 276–277), they indicate that the quintic equation to solve with their practical alkalinity approximation has two positive and three negative roots and that the larger positive one should be used (without any further justification, though). As shown here, this statement is not universally true – there are instances where that equation has only one positive or no positive roots. It is nevertheless true for typical seawater and the lower of the two positive roots actually implies unrealistically low, yet physically sensible, C_{T} (see discussion in Sect. 2.4.2 below).