SolveSAPHE-r2 (v2.0.1): revisiting and extending the Solver Suite for Alkalinity-PH Equations for usage with CO2, HCO3− or CO32− input data
- Dépt. d'Astrophysique, Géophysique et Océanographie, Université de Liège, 4000 Liège, Belgium
Correspondence: Guy Munhoven (email@example.com)
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 (AlkT) and dissolved inorganic carbon (CT), and which converges for any physically sensible pair of such data, has been adapted and further developed to work with AlkT–CO2, AlkT–, and AlkT–. The mathematical properties of the three modified alkalinity–pH equations are explored. It is shown that the AlkT–CO2, and AlkT– problems have one and only one positive root for any physically sensible pair of data (i.e. such that [CO2]>0 and ). The space of AlkT– pairs is partitioned into regions where there is either no solution, one solution or where there are two. The numerical solution of the modified alkalinity–pH equations is far more demanding than that for the original AlkT–CT pair as they exhibit strong gradients and are not always monotonous. The two main algorithms used in SolveSAPHE v1 have been revised in depth to reliably process the three additional data input pairs. The AlkT–CO2 pair is numerically the most challenging. With the Newton–Raphson-based solver, it takes about 5 times as long to solve as the companion AlkT–CT pair; the AlkT– pair requires on average about 4 times as much time as the AlkT–CT 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 AlkT– 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.
Among all the aspects of the ongoing global environmental changes (climate change, ocean acidification, etc.), the solution chemistry of carbon dioxide (CO2) 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: [CO2] (or equivalently the partial pressure of CO2, pCO2, or its fugacity, fCO2), ,  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 CO2 and on one hand, and between and 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), ; (2) total alkalinity, AlkT; (3) pH and (4) pCO2 or fCO2. Recently, a procedure to measure 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 CT and AlkT 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 as a sixth independent variable, although currently not (yet) measurable. Most modellers will call upon CT and AlkT 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 pCO2 directly if that variable is required: the uncertainty of calculated pCO2 is always 5 to 10 times as large as that of the directly measured one.  can be calculated with lower uncertainty from AlkT–CT 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 pCO2 and : 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 AlkT–CT calculated from it are, however, about 20 times as large as those directly measured. The currently most attractive pairs are AlkT–pCO2 and CT–pCO2, both of which allow to calculate pH with better and  with only slightly larger but still acceptable uncertainty than the direct measurement would offer. Direct  measurements, which might be most advisable for tracing carbonate mineral saturation states, are best paired with AlkT or CT (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  calculated from AlkT–CT, 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 AlkT–CT pair which was addressed in full detail by Munhoven (2013) these are (1) AlkT–CO2, (2) AlkT– and (3) AlkT–. 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 AlkT–CT pair. Here, we are revisiting that approach, extending and adapting it so that the AlkT–CO2, AlkT–, and AlkT– 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.
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 right-hand side is that part of the total alkalinity that is not related to the water self-ionisation: , 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, AlknW([H+]) is partitioned into carbonate alkalinity, AlkC([H+]), and non-carbonate alkalinity, AlknWC([H+]), since the relevant carbonate system parameters (the concentrations of CO2, and and their sum, CT) are all directly related to AlkC: . Similarly to AlknW, AlknWC 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 AlknWCinf and AlknWCsup, respectively. Equation (1) is thus formally rewritten as
The carbonate alkalinity term writes, as a function of CT
where K1 and K2 are the first and second stoichiometric dissociation constants of carbonic acid. The individual carbonate species fractions of CT can be expressed as a function of [H+]:
Accordingly, AlkC([H+]) may be rewritten in one of the following forms:
In order to get a first idea about the complications that we might encounter for the solution of the three new data pairs, we start with an exploratory analysis using the Solver Suite for Alkalinity-PH Equations (SolveSAPHE) version 1.0.3 (Munhoven, 2013–2021). The three panels in the upper row of Fig. 1 show the pH and the and concentration distributions for a reduced SW3 test case from Munhoven (2013), with the CT range extending from 0 to 4 mmol kg−1 and the AlkT range from −1 to 3 mmol kg−1 only. The pH distribution from Fig. 1a is then used to derive the corresponding CO2 (not shown), and concentration distributions (Fig. 1b and c). Since we intend to solve the alkalinity–pH equation for given AlkT, and either of [CO2], or , we furthermore produce the concentration isolines for the three species on a pH–AlkT graph (Fig. 1d, e and f). For the three latter, we first calculated AlkC from AlkT by using Eq. (2). Positive AlkC values were then used with Eqs. (7), (8) and (9) to derive the corresponding [CO2], and , respectively. Blank areas represent the pH–AlkT combinations that lead to negative AlkC.
The V- or U-shaped isolines for on the CT–AlkT graph and for on the CT–AlkT and on the pH–AlkT graphs show that the CT– and the AlkT– 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 isoline and on the horizontal line drawn through AlkT=2.3 mmol kg−1. For that pair of data values, there are thus two compatible CT and correspondingly two possible pH values. On the other hand, with the same AlkT, a of 1 mmol kg−1 would not provide any solution as the 1 mmol kg−1 isoline has its minimum at AlkT=2.63 mmol kg−1. Similarly, there are pairs of CT and values that are compatible with two AlkT values and thus two pH values, and others with none: a vertical line drawn through CT=2.2 mmol kg−1 crosses the 2.0 mmol kg−1 isoline for twice and so that pair of data values leads to two pH solutions; a vertical line drawn through CT=2.05 mmol kg−1 does not cross that 2.0 mmol kg−1 isoline for 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 AlkT–CO2 and AlkT– pairs. According to the outcome of our preliminary analysis above, the AlkT– 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 CT– 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 CT by the concentration of one of its individual species.
2.2 The AlkT–CO2 problem
2.2.1 Mathematical analysis
Just like the AlkC([H+]) expression from Eq. (3) is monotonously decreasing with [H+] for CT fixed, that from Eq. (7) is monotonously decreasing with [H+] for [CO2] fixed. The expression on the left-hand side of Eq. (10) decreases from +∞ to −∞ for [CO2]>0 as [H+] varies from 0+ to +∞. Equation (10) thus always has exactly one positive solution.
2.2.2 Root bracketing
i.e. as the positive root of the cubic equation:
Let us denote this cubic by P(H). It is important to notice that and . The equation P(H)=0 has therefore one and only one positive root.
Similarly, the upper bound Hsup 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 Hmin>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 Hmin.
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 Hinf, we may actually choose the Hmin of the first cubic which is lower than the positive root and thus sufficient. Regarding Hsup, it should be noticed that . Accordingly, P(H)>Q(H) for H>Hmin 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 Hsup.
Any bracketing root-finding algorithm can then be used to solve the modified alkalinity pH equation (Eq. 10).
2.3 The AlkT–HCO problem
2.3.1 Mathematical analysis
2.3.2 Root bracketing
The lower bound Hinf can be chosen such that
i.e. as the positive root of the quadratic equation
Similarly, the upper bound Hsup can be chosen such that
i.e. as the positive root of the quadratic equation
Both equations always have two roots: one positive and one negative; their product is negative as indicated by the constant term. With the respective positive roots, we have again bounds for the solution of the modified alkalinity–pH equation and any bracketing root-finding algorithm can be used to solve it.
2.4 The AlkT–CO problem
Whereas any physically meaningful AlkT–[CO2] or AlkT– concentration pairs will always provide one and only one [H+] (or equivalently pH) value as demonstrated above, this cannot be the case for every AlkT– pair, as can be deduced from Fig. 1c and f. On one hand, there are two compatible CT values, and equivalently two pH values, for most AlkT– pairs. This little-known fact was already documented in the 1960s (see, e.g. Deffeyes, 1965)1. On the other hand, there are also AlkT– 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
Let us collect all the terms that are related to carbonate or water self-ionisation alkalinity on the left-hand side, introduce the shorthand 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 . Since K2 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 concentrations of the order of 1 nmol (kg SW)−1 and less.
If γ=0 (i.e. if ), the equation has exactly one root if , and no root otherwise.
If γ>0, the left-hand side is not monotonous: it decreases from +∞ in to a minimum (see below) and then increases back to +∞ as . The right-hand side is bounded and strictly increasing over the same interval (Munhoven, 2013). As a result, the equation has no root if the right-hand side is too low, exactly one if the two curves become tangent and two roots if the right-hand side is great enough.
To alleviate notation let us define the two parametric functions:
where [H+] is the independent variable and γ and A (alkalinity) are parameters. With these two function definitions, Eq. (12) is then rewritten as . Schematic representations of the three γ cases and of the L and R functions are shown in Fig. 2.
The first case can be handled similarly to the AlkT–CO2 and AlkT– 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 Hinf and
for Hsup, and retaining the respective positive roots of each.
The second case might be considered to be only mathematically of importance as it only applies for one exact (and thus improbable) 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 , and no solution otherwise. The root can be easily bracketed from below. It is sufficient to choose Hinf such that
leading to . The analogue equation for Hsup, with AlknWCinf replaced by AlknWCsup (cf. Eqs. 15 and 16) does not work if . The newly derived asymptotic approximation for AlknWC([H+]) as (see the “Mathematical and Technical Details” report in the Supplement) nevertheless provides a means to derive an upper bound. It is sufficient to choose Hsup 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 K1,[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 , it is straightforward to show that with this choice. Equation (12), which is equivalent to thus has a single root between Hinf and Hsup.
The third case is the most commonly encountered and the most challenging. With γ>0, 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 Hmin and the value that L takes there by Lmin:
There are two ranges of AlkT values where firm conclusions can be drawn right away.
If , i.e. if , Eq. (12) has two distinct roots, since R(H;AlkT) is bounded. Furthermore, the roots – let us provisionally denote the lower one H1 and the greater one H2 – are such that H1<Hmin and H2>Hmin. Hmin can thus be used as an upper bracket for H1 and as a lower bracket for H2. However, if , the abscissae of the intersection points PLL and PLR (see Fig. 2), which are solutions of
provide tighter brackets than Hmin.
If , i.e. if , Eq. (12) does not have any roots.
For intermediate values of AlkT, no firm quantitative statement regarding the root(s) of Eq. (12) can be made a priori. As AlkT decreases from Lmin+AlknWC(Hmin) to Lmin+AlknWCinf, Eq. (12) will at first still have two roots, but both are greater than or equal to Hmin. At some intermediate value, and become tangent. At this point, Eq. (12) has one double root, which is the abscissa of that tangent point, Htan. Htan is actually a universally valid separation limit between two roots, if there are any. For lower values of AlkT, the problem does not have any solutions.
The limiting AlkT value for which the two curves are tangential and the corresponding Htan 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 not as an equation in [H+] for given parameter values γ (or, equivalently, ) 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 and to contribute to the solution of the minimisation of A([H+]). The determination of Htan 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 AlkT values for which the exact knowledge of Htan is not indispensable. In these situations Hmin 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, Htan 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 Htan, we can use the three characteristic [H+] values from Fig. 3 as initial conditions. These are Hmin together with the abscissae HL and HR of the intersection points between and the horizontal line at Alkmin−AlknWCinf, which are the roots of
By construction, . 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 with respect to [H+] is positive provided that the successive dissociation constants Kj,[i] of the different acid systems (denumbered by i) resulting from the dissociation of an acid are such that – 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 . The underlying technical developments can be found in the “Mathematical and Technical Details” report in the Supplement. is thus concave, while is convex for γ>0. A([H+]) thus has only a single local minimum comprised between HL and HR.
Once Htan is known, the root brackets can be completed by the intersection points between and the horizontal line at AlkT−AlknWCinf – corresponding to the PUL and PUR points in Fig. 2 with the grey band shifted down to include the minimum – i.e. by solving the same quadratic equation than for HL and HR, with Alkmin replaced by AlkT. We have again , 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 AlkT– 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 CT 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 AlkT– problem above, we determined the AlkT ranges that would, respectively, lead to two, one or no roots, for a given γ, i.e. . 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 evolves as a function of pH for a given value of AlkT, and determine the 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 AlkT fixed at 2.3 mmol kg−1 and draw the evolution of as a function of pH following
obtained by first using Eq. (9) to express as a function of AlkC and [H+], and then Eq. (2) to calculate AlkC as a function of AlkT and [H+] (and the total concentrations of all the other contributing acid–base systems, which we assume to be known, as stated initially). CO2 and evolution curves can be derived similarly, by using Eqs. (7) and (8), respectively, instead of Eq. (9). The concentration evolution for the other AlkT 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 AlkT=2.3 mmol kg−1; Fig. 4b shows the evolution curves for different AlkT values, ranging from 0.5 to 2.5 mmol kg−1. Solving the AlkT– problem for our showcase sample where 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 curve, if any. There are actually two of them, located at pH = 8.03 and pH = 11.43. With increasing target values for , i.e. moving the horizontal line upward, the two pH roots will move closer and closer together, until the maximum of the CO3 curve is touched, at CO3 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 CO3 curve any more: there are no roots for . As illustrated in Fig. 4b, the value of CO3 max grows as AlkT increases, thus extending the range of  that allows for roots. Another noteworthy fact in Fig. 4b is that all the displayed 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 are monotonously decreasing), all of them nevertheless go to zero. The equation actually always has exactly one root, for any given AlkT, 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 CT=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 AlkT value (Munhoven, 2013). The pH value at that zero-crossing point is furthermore the maximum possible pH for the given AlkT: beyond that value, the ever-growing [OH+] would inevitably make AlkT increase above the fixed value, thus requiring AlkC to become negative, which is not possible.
As can be seen in Fig. 4a, the high-pH solution goes together with : represents only 4.6 % of the CT at pH = 8.03, typical for seawater, but 99.2 % at pH = 11.43. Accordingly CT=0.1008 mmol kg−1 at the high-pH root, which is unrealistically low for many natural samples. In the marine realm, this observation regarding the high-pH root is actually correct in general. 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 CO3 curves are greater than 0.47 mmol kg−1 for AlkT≥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 mmol kg−1 and AlkT≥1.5 mmol kg−1, the greater of the two pH roots always implies that represents more than 80% of CT. 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 . 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 AlkT 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 (AlkT–CO2, AlkT– and AlkT– 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 Hinf, we use H0=Hinf; if it is greater than Hsup, we set H0=Hinf. For problems that have two positive [H+] solutions (AlkT– with γ>0 and sufficiently great AlkT), 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.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 AlkT–CO2 pair, the Newton–Raphson-based algorithm showed a few weaknesses. With the AlkT–CT pair that SolveSAPHE v1 had been designed for, each non-water alkalinity term was bounded, just like its derivative. Once CO2 takes the role of CT these favourable properties are lost: with [CO2] 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 CO2 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 CT 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 ϵ ( by default). However, iterations for AlkT–CO2 and for the greater root of AlkT– were prone to early termination with that stopping criterion, as iterates only slowly changed due to the extreme gradients in the AlkC 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 , where Hmax and Hmin 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 [Hmin,Hmax]. 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 Hmax and Hmin have strongly different magnitudes. Unfortunately, the number of iterations required for the original SolveSAPHE pair AlkT–CT 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 AlkT–CT is generally the relevant pair of data, SolveSAPHE v1 remains the most efficient choice. Tests have shown that the two safeguarded algorithms from SolveSAPHE v1 typically require 40 %–45 % less computing time than their SolveSAPHE-r2 counterparts.
Finally, as explained above, some AlkT– 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 AlkT–CO2, AlkT– and AlkT– 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 AlkT and CT 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 CO2 and , 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 [CO2],  or  distributions are adapted. Although the [CO2], and ranges for each test case reported on Table 1 have been defined on the basis of their respective distributions calculated from the AlkT–CT 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 AlkT–CT pair and the results stored. For the other three pairs, the pH distribution obtained with the AlkT–CT pair for the corresponding set of temperature, salinity and pressure is first read in and the corresponding [CO2], or distribution calculated on the underlying CT–AlkT grid. The so-obtained arrays of species concentrations are then used to define the set of AlkT–CO2, AlkT– and AlkT– 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 [CO2]–AlkT, –AlkT or –AlkT grids, respectively, which is nevertheless irrelevant for the histogram syntheses presented in Figs. 6 and 7. These variants of the test cases are denoted SW1CT, SW2CT, SW3CT, BW4CT and ABW5CT.
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 AlkT–CO2 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, AlkT–CO2 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 AlkT–CO2 samples requires 45 to 55 and more iterations.
In comparison, AlkT–CT 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, AlkT– coming closest to AlkT–CT.
ABW5 shows a few deviations from the other tests cases.
Here, solving the AlkT–CO2 problem with
solve_at_general nearly always takes more than 50 iterations, with
solve_at_general almost always nine.
The solution of AlkT–CT for ABW5 with
solve_at_general_sec takes considerably more iterations than AlkT– (the fastest), and AlkT–CO2 (the second fastest).
Finally, Figs. 6 and 7 demonstrate the superiority of
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 AlkT–CT problem than
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 AlkT–CO2 than with AlkT–CT; for AlkT–, it takes 4 times as much (for both roots though, including the solution of the minimisation problem for part of the domain). For AlkT–, the difference is only 20 %. With the secant-based method, the picture is completely different: AlkT–CO2 takes only about 30 % more time than AlkT–CT, AlkT– twice as much, whereas AlkT– executes even about 5 % faster. For the AlkT–CO2 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 AlkT–CT and the AlkT– input pairs, compared to the standard cubic polynomial one. Similar increases are obtained with a constant uniform pH = 8 initialisation. For AlkT–CO2, and AlkT–, 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 AlkT have been made out for γ>0: an upper one at Lmin+AlknWC(Hmin), above which the problem always has two [H+] solutions, and a lower one at Lmin+AlknWCinf, below which the problem does not have any solution at all. For intermediate values of AlkT, it is necessary to determine Htan and Alktan 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 Htan 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 AlkT– pair and the solution of the minimisation problem is required for 173 445 samples (8.89 %), because Htan 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 AlknWC minimum in seawater, because of the high total sulfate concentration in seawater (). For carbonate ion concentrations below 400 µmol kg−1, i.e. for most of the naturally occurring waters, the AlkT– 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 [CO2],  and  instead of the total inorganic carbon concentration, CT. The rationale can be entirely transposed to these three pairs: (1) the amended alkalinity–pH equations for AlkT–CO2 and for AlkT– still have one and only one positive solution, while AlkT– 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 AlkT–CO2 and AlkT– problems, as well as for the two roots of AlkT– problems that may have to be solved for naturally occurring sample compositions. More uncommon but physically realistic AlkT– problems may additionally require the solution of an auxiliary minimisation problem to determine the threshold AlkT 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 AlkT– 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 AlkT–CO2, AlkT–, and AlkT– problems. Most bisection steps have been replaced by regula falsi steps for increased convergence speed. The secant–bisection solver only required minimal adaptations. A Fortran 90 reference implementation, SolveSAPHE-r2, was prepared and used to evaluate the performance of the different methods for solving four benchmark problems. While the secant–bisection method was already slightly superior to the Newton–Raphson–bisection method in SolveSAPHE v1, that advantage has now become overwhelming: in SolveSAPHE-r2, it typically requires 2 to 4 times fewer iterations, and for the newly handled pairs, the equation residuals are orders of magnitude lower than the Newton–Raphson–regula falsi–bisection-based solver (typically of the order of 10−19–10−18 compared to 10−13–10−12).
For carbonate speciation problems posed by AlkT and either one of [CO2],  or , the secant-based routine from SolveSAPHE-r2,
solve_at_general2_sec, is thus clearly the method of choice; for calculations on the basis of AlkT–CT, both
solve_at_general_sec from SolveSAPHE v1 will perform better, although the secant-based 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 CO2 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 [CO2] and the fugacity (fCO2) or the partial pressure (pCO2) of CO2 and between pH and [H+] on any chosen pH scale. We therefore restrict ourselves to [CO2] 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.
- CT and CO2, CT and CO
: (1) With these two pairs, the fraction (the 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) AlkT from Eq. (1). In addition to the positivity of the species concentrations, the following constraints must be met: [CO2]<CT and .
- CT and HCO
: With this pair, the fraction, denoted by b hereafter, is fixed, and Eq. (5) becomes a quadratic equation in [H+]. That equation has two positive solutions if , one double root if and no real solutions if . This is well illustrated in Fig. 1b above: there are CT– combinations that allow for two different AlkT values and, equivalently, two pH values, and there are others that do not allow for any. Please notice that the threshold fraction 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 AlkT– problem: which one to choose? Most often, the lower of the two will again be the appropriate one, as that one typically leads to AlkT>CT, whereas the greater one leads to AlkT<CT. 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 [CO2] and from their respective species fractions; (3) AlkT from Eq. (1).
- CO2 and HCO
: (1) [H+] from K1; (2)  from K2; (3) CT can be calculated from the three carbonate species concentrations; (4) AlkT from Eq. (1).
- CO2 and CO
: (1)  from ; (2) CT from the three carbonate species concentrations; (3) [H+] from K1 or K2; (4) AlkT from Eq. (1).
- HCO and CO
: (1) calculate [H+] from K2; (2) [CO2] from K1; (3) CT from the three carbonate species concentrations; (4) AlkT from Eq. (1).
- CO2 and H+
: (1) calculate  from K1; (2) calculate  from K2; (3) CT from the three carbonate species concentrations; (4) AlkT from Eq. (1).
- HCO and H+
: (1) calculate [CO2] from K1; (2) calculate  from K2; (3) CT from the three carbonate species concentrations; (4) AlkT from Eq. (1).
- CO and H+
: (1) calculate  from K2; (2) calculate [CO2] from K1; (3) CT from the three carbonate species concentrations; (4) AlkT from Eq. (1).
- AlkT and H+
: (1) CT from Eq. (2); (2) individual species concentrations from the species fractions.
As illustrated in Fig. 1d–f above, there are AlkT–[H+] combinations that lead to physically unrealistic negative AlkC. Following Eq. (3), negative AlkC requires negative CT, 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 AlkT contributors is inappropriate.
- CT and H+
: individual species concentrations from the species fractions; AlkT 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 SolveSAPHE-r2 (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 bug-fix releases and updates will also be archived there.
The supplement related to this article is available online at: https://doi.org/10.5194/gmd-14-4225-2021-supplement.
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 kick-off momentum for me to reconsider SolveSAPHE and to complete it in order to combine AlkT with the individual carbonate system species instead of CT – thank you. Thanks are also due to Jean-Pierre Gattuso for pointing out issues with the colour schemes adopted in the originally submitted version of the manuscript and for guidance about how to address them. The comments and recommendations by the two anonymous referees were most helpful to improve the framing of the paper and the scope of the results. The analysis based upon Fig. 4 was suggested by Anonymous Referee #1 and recommended for inclusion in the manuscript by Anonymous Referee #2.
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.
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
Dickson, A. G., Sabine, C. L., and Christian, J. R. (Eds.): Guide to Best Practices for Ocean CO2 Measurements, vol. 3, in: PICES Special Publication, Carbon Dioxide Information and Analysis Center, Oak Ridge (TN), available at: https://cdiac.ess-dive.lbl.gov/ftp/oceans/Handbook_2007/Guide_all_in_one.pdf (last access: 24 June 2021), 2007. a, b, c, d
Munhoven, G.: Future CCD and CSH variations: Deep-sea 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/gmd-6-1367-2013, 2013. a, b, c, d, e, f, g, h, i, j, k, l, m, n, o, p, q
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
Zeebe and Wolf-Gladrow (2001) appear to be aware of it. In their recipe for given AlkT and  (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, CT (see discussion in Sect. 2.4.2 below).