Mathematics of the total alkalinity – pH equation – pathway to robust and universal solution algorithms : the SolveSAPHE package v 1 . 0 . 1

The total alkalinity–pH equation, which relates total alkalinity and pH for a given set of total concentrations of the acid–base systems that contribute to total alkalinity in a given water sample, is reviewed and its mathematical properties established. We prove that the equation function is strictly monotone and always has exactly one positive root. Different commonly used approximations are discussed and compared. An original method to derive appropriate initial values for the iterative solution of the cubic polynomial equation based upon carbonate-borate-alkalinity is presented. We then review different methods that have been used to solve the total alkalinity–pH equation, with a main focus on biogeochemical models. The shortcomings and limitations of these methods are made out and discussed. We then present two variants of a new, robust and universally convergent algorithm to solve the total alkalinity–pH equation. This algorithm does not require any a priori knowledge of the solution. SolveSAPHE (Solver Suite for Alkalinity-PH Equations) provides reference implementations of several variants of the new algorithm in Fortran 90, together with new implementations of other, previously published solvers. The new iterative procedure is shown to converge from any starting value to the physical solution. The extra computational cost for the convergence security is only 10–15 % compared to the fastest algorithm in our test series.


Introduction
Biogeochemical models have become indispensable tools to improve our understanding of the cycling of the elements in the Earth system. A central and critical component of almost all biogeochemical models is the pH calculation routine. In ocean carbon cycle models, the air-sea exchange of CO 2 is directly linked to the surface ocean [CO 2 ]; the preservation of biogenic carbonates in the surface sediments at the sea floor is closely linked to the deep sea [CO 2− 3 ] (Broecker and Peng, 1982). The fractions of CO 2 , HCO − 3 and CO 2− 3 in the total dissolved inorganic carbon (i.e. the speciation of the carbonate system) are controlled by pH. Hence, pH changes in seawater may directly influence air-sea exchange of CO 2 or the preservation of carbonates in the deep sea. Conversely, the dissociation of acids, such as carbonic acid, also controls pH: when the ocean takes up or releases CO 2 (e.g. as a result of a rise or a decline of the abundance of CO 2 in the atmosphere), its pH changes. The currently ongoing ocean acidification due to the massive release of CO 2 into the atmosphere by human activity is but one example of such an induced pH change.
The nitrogen cycle is another important biogeochemical cycle where pH plays an important role. The speciation of dissolved ammonium is -as that of any acid-base systemdependent on pH, NH 3 being more abundant than NH − 4 at high pH, and less abundant at low pH. At pH > 9, the concentration of NH 3 in seawater may reach toxic levels.

G. Munhoven: Solving the alkalinity-pH equation: SolveSAPHE
The realistic modelling of biologically mediated fluxes (e.g. marine primary or export production) requires the colimitation or even inhibition by different chemical components to be taken into account. The nitrogen and carbon cycles, already mentioned above, strongly interact, both in the ocean and on land. In the ocean, Fe and other metals act as micronutrients and once again, pH plays an important role as the solubility of metals is strongly dependent on pH (Millero et al., 2009). The resulting coupling of the biogeochemical cycles of different elements makes biogeochemical models become more and more complex and pH calculation more and more difficult.
Biogeochemical models are now increasingly used for settings that are strongly different from present day. Typical applications include future ocean acidification (e.g. Caldeira and Wickett, 2003), the Paleocene-Eocene Thermal Maximum (e.g. Ridgwell and Schmidt, 2010), Snowball Earth (e.g. Le Hir et al., 2009), etc. Some commonly used pH solvers may possibly become unstable and produce unreliable results. The convergence properties of currently used solution methods has actually never been systematically tested.
Unfortunately, information on pH solver failures is only seldom published. Zeebe (2012) reports for his LOSCAR model that negative H + concentrations may be obtained when starting with total alkalinity and dissolved inorganic carbon concentrations in a very high ratio, requiring the model run to be restarted with the respective concentrations in a lower ratio. Hofmann et al. (2010) also indicate that the standard pH solving routine in their R modelling environment AquaEnv may fail when trying to calculate the pH for samples with very low or zero dissolved inorganic carbon concentrations. In this case, they resort to a general purpose interval based root finding routine instead, adopting a very large bracketing interval (see Hofmann et al., 2012), possibly leading to a considerable performance loss. Andy Ridgwell, in his editorial comment to the companion discussion paper, mentions convergence problems encountered with the GE-NIE model code (Ridgwell et al., 2007) encountered while studying the effect of an artificial addition of lime (CaO) to the ocean surface (a particular geoengineering method meant to accelerate the uptake of carbon dioxide from the atmosphere) once total alkalinity came to exceed the typical surface ocean concentrations of dissolved inorganic carbon by about a factor of two. As we will show below, the three models use essentially equivalent pH calculation methods, which become divergent under those typical conditions. The speciation of any acid system, i.e. the determination of the concentrations of each one of the undissociated and the different dissociated forms of an acid, is an underdetermined problem if only the total concentration and thermodynamic or stoichiometric constants are known. This underdetermination can be lifted if pH is known. Being dependent on temperature and pressure, neither pH nor [H + ] are, however, well suitable for being used in transport equations, and thus in biogeochemical models. In biogeochemical models, the common way to resolve this underdetermination is to consider another conservative quantity: total alkalinity, also called titration alkalinity. Total alkalinity, which is also an experimentally measurable quantity, ties all the different acid systems present in a water sample together and allows us to solve the speciation problem. In comparison to pH, it has the advantage of being a conservative quantity: it is only controlled by its sources and sinks, and it is independent on temperature and pressure (Zeebe and Wolf-Gladrow, 2001).
In the following section, we provide a comprehensive introduction to the concept of alkalinity. In our exploration of the mathematical properties of the equation that relates [H + ] to total alkalinity start with a detailed presentation of various approximations commonly used for present-day seawater. The analysis of the mathematical properties of these approximations will provide useful hints for the characteristics of the general case. In Sect. 3 of this paper, we present solution methods for deriving pH from each of the various approximations to total alkalinity considered. Complications that might possibly arise from the various pH scales that are in use in marine chemistry are elucidated in Sect. 4. In Sect. 5, we then show that there are intrinsic bounds that bracket the root of the total alkalinity-pH equation, and that can be directly derived from the approximation used to represent total alkalinity. The existence of such bounds makes it possible to define a new, universal algorithm to solve the alkalinity-pH equation, which requires no a priori knowledge of the root. A reference implementation of two variants of the new algorithm is presented in Sect. 6. The algorithms are tested for their efficiency and robustness and their performance compared with that of the most common previously published general solution methods.

Total alkalinity: general definition and approximations
In the following parts of this section, we review a number of aspects of total alkalinity in natural waters. The main focus will be put onto seawater and on the carbonate system, but all the presented developments can be applied to any natural water sample, provided the required thermodynamic constants are known. We briefly recall the different approximations commonly used for calculating pH and the speciation of acid systems. We will then establish a few basic properties of the expressions that relate the various types of alkalinity to total concentrations and pH. Although simple, these properties do not seem to have been previously explored in detail, nor exploited for designing methods of solution of the alkalinity-pH equation.
Although we primarily focus on modelling in the following developments, the calculation procedures are obviously also applicable in experimental set-ups.

General definition
Total alkalinity, also called titration alkalinity, denoted here Alk T , reflects the excess of chemical bases of the solution relative to an arbitrary specified zero level of protons, or equivalence point. Ideally, Alk T represents the amount of bases contained in a sample of seawater that will accept a proton when the sample is titrated with a strong acid (e.g. hydrochloric acid) to the carbonic acid endpoint. That endpoint is located at the pH below which H + ions get more abundant in solution than HCO − 3 ions; its value is close to 4.5. H + added to water at this pH by adding strong acid will remain as such in solution. Please notice that, for the sake of a simpler notation, we follow here the common usage of denoting protons in solution by H + , although free H + ions sensu stricto do only exist in insignificantly small amounts in aqueous solutions. Each proton is rather bound to a water molecule to form an H 3 O + ion, and each of these H 3 O + in turn is furthermore generally hydrogen bonded to three other H 2 O molecules to form an H 9 O + 4 ion (Dickson, 1984). Rigorously speaking, Alk T is defined as the number of moles of H + ions equivalent to the excess of "proton acceptors", i.e. bases formed from acids characterized by a pK A ≥ 4.5 in a solution of zero ionic strength at 25 • C, over "proton donors", i.e. acids with pK A < 4.5 under the same conditions, per kilogram of sample (Dickson, 1981).
With emphasis on the most important contributors, a rather complete expression for Alk T in a seawater sample is where the ellipses refer to other potential proton donors and acceptors generally present at negligible concentrations only. All of the concentrations are total concentrations (which include free, hydrated and complexed forms of the individual species), except for [H + ] f , which only includes the free and hydrated forms. There are alternative definitions that can be found in the literature, which lead to similar, although not necessarily exactly the same, expressions. However, the above definition is the one that reflects the titration procedure used to measure alkalinity the most accurately. We will therefore base the following developments upon it.
In other natural water samples (lake, river, or brines) the constituent list in Eq. (1) needs to be adapted: some constituents may be neglected and bases of other acid systems have to be included (e.g. bases derived from organic acids, from dissolved metals, etc.). While total alkalinity in seawater samples typically ranges between about 2 and 2.6 meq kg −1 , acid mine drainage samples may even present negative alkalinity, representing the fact that a strong base instead of a strong acid must be added to reach the equivalence pH point of 4.5. Interested readers may refer, e.g. to Kirby and Cravotta III (2005) and references therein for such -from a marine chemist's point of view -exotic samples.

The pH-total alkalinity equation
Total alkalinity as defined above is a conservative quantity with respect to mixing, changes in temperature and pressure (Wolf-Gladrow et al., 2007). It is therefore a cornerstone in biogeochemical cycle models which are most conveniently formulated on the basis of conservation equations. In such models, definition/Eq. (1) above, or an adequate variant, is used to solve the inverse problem for [H + ]. All of the individual species concentrations appearing in Eq. (1) can be expressed in terms of the total concentrations of the acid systems that they respectively belong to and of [H + ]. Given the evolutions of the total concentrations of all the acid systems considered (dissociated and non-dissociated forms) and of Alk T -all of which can be derived from appropriate conservation equations -expression (1) is interpreted as an equation for [H + ] or, equivalently, pH. We will therefore call that equation the total alkalinity-pH equation.
We might actually have called our equation simply the pH equation. Alk T does indeed not play any special or more important role than any of the total concentrations of the other acid systems considered. We do, however, feel that this name would have been too general and thus prefer to include "total alkalinity" in the name to reflect that the overall structure of the equation derives from the definition of total alkalinity.

Typical applications in biogeochemical models
In a typical global ocean carbon cycle model, total alkalinity may commonly be approximated by (2) provides a relationship between C T , B T , Alk T and [H + ] (i.e. pH). The model provides conservation equations for C T and Alk T ; B T can generally be taken proportional to salinity, whose evolution either follows a prescribed scenario or may also derived from a conservation equation. Relationship (2) thus reduces to an equation in [H + ]. The solution of that equation finally provides a means to calculate the complete speciation of the carbonate and the borate systems.
In other biogeochemical studies where other systems than the carbonate system are of interest (such as ammonium, sulphides, etc.), the procedure is entirely analogue. Each one of the individual species concentrations that need to be considered in Eq. (1) for that particular application is expressed in terms of the total concentration of the acid system that it belongs to and conservation equations, scenarios or measurements that are used to evaluate all of the total concentrations, including total alkalinity. These steps again reduce Eq. (1) into an equation in [H + ], whose solution provides a direct means to calculate the speciations of all the systems considered.

Common approximations for total alkalinity in seawater
Here we first analyse the forward problem for a few specific approximations used for seawater: for given total concentrations of dissolved inorganic carbon, total borate, etc., we analyse how the expressions for the different types of alkalinity change as a function of [H + ]. This simple analysis will already provide valuable insight into the overall mathematical properties of the total alkalinity-pH equation and its subcomponents, which we can exploit later for the most general case.

Carbonate alkalinity
The contribution of the carbonic acid system (or carbonate system) to total alkalinity is called carbonate alkalinity and we denote it by Alk C :

Carbonate and borate alkalinity
The second most important component of natural presentday seawater alkalinity is borate alkalinity, Alk B . Together with the carbonate alkalinity we have Upon substitution of the individual species concentrations by their fractional expressions as a function of [H + ], we get where B T is the total concentration of dissolved borates and K B is the dissociation constant for boric acid. For constant B T , Alk B is again a strictly decreasing function with [H + ], similarly to Alk C . Hence, for constant C T and B T , Alk CB is a strictly decreasing function with [H + ] and, as a consequence, 0 < Alk CB < 2C T + B T as long as C T + B T = 0.

Carbonate, borate and water self-ionization alkalinity
In a third stage, we may consider the alkalinity that arises from the dissociation of the solvent water itself (by self-ionization) in addition to carbonate and borate alkalinity and get the next important approximation for natural present-day seawater, called practical alkalinity by Zeebe and Wolf-Gladrow (2001): Upon substitution by the respective speciation relationships, we get

Contribution of a generic acid system to total alkalinity
In common seawater, Alk CBW is entirely sufficient even for applications that require high accuracy. However, in some cases other systems than the carbonate and borate systems need to be considered. This is especially the case in suboxic and anoxic waters, such as semi-closed fjords (e.g. Framvaren Fjord in Norway studied by Yao and Millero, 1995) or at a larger scale, the Black Sea (e.g. Dyrssen, 1999), where, e.g. the contribution from sulphides cannot be neglected. In order to generalize our analysis of the total alkalinity-pH equation, let us consider a generic acid, denoted by H n A, that may potentially lead to n successive dissociation reactions, characterized by stoichiometric dissociation constants K 1 , K 2 , . . ., K n , respectively: For simplicity, we omit the " * " superscript commonly used elsewhere to differentiate stoichiometric from thermodynamic dissociation constants (i.e. elsewhere stoichiometric constants generally write K * i instead of K i ). Throughout this paper, the constants used will relate concentrations and not activities. As such, they include the effect of activity coefficients that differ from unity. The values of such constants not only depend on temperature and pressure but also on the ionic strength of the solution. Everything developed here furthermore applies to all kinds of acids, be they of Arrhenius, Brønsted-Lowry, Lewis or any other type, even if the adopted notation could possibly suggest that our developments only apply to Arrhenius-type acids.
The joint contribution of all the different dissociated and nondissociated forms of H n A to alkalinity, proton donors and proton acceptors alike, is then equal to where m is an integer constant, which is dependant on the socalled zero proton level of the system under consideration: m is such that pK m < 4.5 < pK m+1 if pK 1 < 4.5 and pK n > 4.5 Since pK m < 4.5, all of the H n−j A j − in the H n A -. . . -A n− system for j = 0, . . ., m − 1 are proton donors: the last one (j = m − 1) has a strength of 1 eq mol −1 , the second to last one (j = m − 2) of 2 eq mol −1 , etc. Since pK m+1 > 4.5, the dissociation products H n−j A j − for j = m + 1, . . ., n are proton acceptors, the first one (j = m + 1) with a strength of 1 eq mol −1 , the second one (j = m + 2) with a strength of 2 eq mol −1 , etc. For the carbonic acid system, e.g. n = 2 and m = 0; for the boric acid system, n = 1 and m = 0; for the phosphoric acid system, n = 3 and m = 1.

G. Munhoven: Solving the alkalinity-pH equation: SolveSAPHE
From the previous expressions for the species fractions, we then find that where we have defined ] → +∞, also not actually reachable); both of these could, theoretically, be negative if m is sufficiently large. i] ] and with zero proton levels respectively characterized by m i , the total alkalinity-pH equation,

For a water sample that contains a set of acids H
has exactly one positive root [H + ], for any given value of Alk T : the sum of the respective alkalinity contributions over the set {H n i A [i] |i=1, . . .} of all the acid systems active in the sample is a strictly decreasing function of [H + ]; the contribution from the dissociation of water is also strictly decreasing with [H + ], and may theoretically take any value between +∞ and −∞.

Alkalinity-pH equation in biogeochemical models: approximations and methods of solution
In this section, we are going to review the most common approximations used in ocean carbon and biogeochemical cycle models, focusing on how the corresponding equation is solved.

Carbonate alkalinity based solutions
The straight approximation Alk T Alk C is often used in textbooks (e.g. Broecker and Peng, 1982). There are only a few models (e.g. Opdyke and Walker, 1992;Walker and Opdyke, 1995) that use it directly for their carbonate chemistry speciation. For numerical modelling purposes, its usage is indeed somewhat problematic. [H + ] calculated from Alk T and C T data, by assuming that Alk C = Alk T are typically 30-40 % too low (i.e. 0.15-0.2 pH units too high) for present-day seawater samples. Furthermore, the sensitivity of the C T -Alk C system to perturbations is stronger than that of the C T -Alk CBW system: equilibrium pCO 2 changes, e.g. are of the order of 20 % larger (Munhoven, 1997). The calculation of [H + ] from C T -Alk C remains nevertheless important, as more advanced methods such as those proposed by Bacastow (1981), Peng et al. (1987) or Follows et al. (2006), where Alk C is iteratively recalculated from more complete approximations to Alk T (ICAC methods -see below), rely on it.

Fundamental solution
For given Alk C and C T (C T > 0), the equation to solve for Following our discussion in Sect. 2.2.1, Eq. (7) has a positive root if and only if 0 < Alk C < 2C T ; if there is a positive root, it is unique. Equation (7) can be directly solved after conversion to the quadratic equation: where For valid Alk C values (i.e. for 0 < Alk C < 2C T ), this quadratic equation has two real roots, a positive and a negative one. The positive root is where For Alk C values that are out of range Eq. (8) either has two negative or two complex roots.

Alternative methods
There are other methods to derive [H + ] from Alk C and C T . All of them ultimately seem to rely on the formulae of Park (1969) for deriving the complete speciation of the carbonate system directly from Alk C and C T , without explicitly using [H + ]. Antoine and Morel (1995) first calculate [CO 2 ] from C T and Alk C (which involves the solution of a first parabolic equation), and then derive [H + ] from the relation- , which requires the solution of a second parabolic equation. Ridgwell (2001) first determines the complete speciation of the carbonate system, referring for the adopted procedure to Millero and Sohn (1992), who actually only report the formulae of Park (1969). He then derives two different estimates for [H + ], based upon the definitions of the first and second dissociation constants of carbonic acid, and finally uses the geometric mean of these two estimates as a solution for Eq. (7).
There are no obvious advantages for calling upon these methods instead of the direct quadratic solution above. Even if carefully implemented, both require a significantly higher number of operations than the solution outlined above. Those methods offer a direct access to carbonate speciation (at least in part), which can, however, also be calculated at little extra cost from [H + ].

Iterative carbonate alkalinity correction methods
In most common natural settings, the difference between Alk C and Alk T , albeit small, leads to significant errors on [H + ], if Alk T is used in place of Alk C and one of the procedures above is used to calculate it from C T . To overcome this problem, Alk C can be estimated from Alk T , and then iteratively corrected until stabilization occurs. Such a procedure, which we call here iterative carbonate alkalinity correction (ICAC) can a priori be used with arbitrary chemical compositions, provided Alk C represents a significant fraction of Alk T . If Alk C makes up only a small fraction of Alk T , the method is likely to exhibit unstable behaviour.
In the most straightforward ICAC method, one starts from a trial value H 0 for [H + ], a first estimate Alk C,0 is obtained by subtracting the concentrations of all non-carbonate components from Alk T . That Alk C,0 is then used to calculate a new (improved) estimate H 1 for [H + ] from Eqs. (9) and (10) or one of the alternative methods. H 1 is then used to calculate a new estimate Alk C,1 from Alk T as above and the procedure is iterated until some predefined convergence criterion is fulfilled. This procedure is a classical fixed-point iteration: In this recurrence, Alk C (Alk T , H n ) is the estimate of Alk C obtained from Alk T by subtracting all the non-carbonate components estimated by using H n . Pure fixed-point iterative schemes may be prone to convergence problems (slow convergence or no convergence at all). If the procedure is convergent, the rate of convergence is linear. This plain fixed-point-iteration ICAC method was recently made popular again by Follows et al. (2006). These authors argue that in carbon cycle model simulation experiments, where there is little change in pH from one time step to the next, a single iteration may already provide a sufficiently accurate estimate of [H + ] to derive acceptable pCO 2 estimates, for any chosen approximation of total alkalinity. Follows et al. (2006) suggest, if necessary, to repeat the fixedpoint iteration until a sufficiently accurate estimate is found.
There are a number of models that rely on the ICAC approach for their pH determination. Peng et al. (1987) consider Alk CBW plus the contributions from silicic and phosphoric acid systems in their representation of total alkalinity. 1 They use an initial value of 10 −8 and stop their iterations once |( H )/H | < 0.005 %. They report that less than ten iterations are generally sufficient. Antoine and Morel (1995) adopt Alk CBW as an approximation to Alk T . At each step, they derive [H + ] from C T and Alk C by using their special procedure described above. They iterate until two successive Alk C estimates differ by less than 10 −8 (no units given). Ridgwell (2001) as an approximation to total alkalinity. He calculates [H + ] at each step from C T and Alk C by using his own procedure described above. GENIE (Ridgwell et al., 2007) initially used the same procedure as Ridgwell (2001) ] as an approximation to total alkalinity in GEOCLIM reloaded. They continue to iterate until |Alk CBW +[HS − ]−Alk T | < 10 −6 (no units given). The method is further used in LOVECLIM (A. Mouchet, personal communication, 2012) with Alk CBW as an approximation for total alkalinity (Goosse et al., 2010) and most probably still in some others that, unfortunately, do not provide details about the calculation procedures adopted. Bacastow (1981) proposed a variant to improve the rate of convergence of fixed-point iterations. That variant only uses the recurrence described above for the first two iterates. From the third iteration on, Bacastow (1981) switches to a secant method to solve the fixed-point equation H − Q(Alk C (Alk T , H )) = 0. 2 Fixed-point iterations are thus only used to provide starting values for the solution of the fixedpoint equation by the secant method. The rate of convergence of the method is strongly increased by this approach (and the domain of convergence slightly enlarged -see numerical experiments below). However, for some C T -Alk T combinations the underlying fixed-point equation may still give rise to convergence problems, even with the secant method. However as will be shown below, the method of Bacastow (1981) is strongly preferable over the pure fixed-point scheme.
The Hadley Centre Ocean Carbon Cycle (HadOCC) model (Palmer and Totterdell, 2001) uses Bacastow's method for its carbonate speciation calculation, with the Alk CBW approximation.

Carbonate and borate alkalinity based solution
Only a few models appear to use pH calculation routines based upon Alk CB . MBM-MEDUSA (Munhoven and François, 1996;Munhoven, 1997Munhoven, , 2007 is one of them, the model of Marchal et al. (1998) is another one.

Basic formulation and solution methods
The equation to solve for [H + ] is, for given Alk CB , C T and B T , This equation may be converted into the polynomial equation: with Following our discussion in Sect. 2.2.2, Eq. (12) has a positive root if and only if 0 < Alk CB < 2C T + B T ; if there is a positive root, it is unique. The same holds for the cubic Eq. (13). The cubic equation could possibly be solved with closed formulae, such as Cardano's formulae (which may, however, suffer from precision problems, require numerically expensive cubic root evaluations or possibly complex arithmetic) or Viète's trigonometric formulae (which require a combination of an arccosine, a cosine and a square root). When adopted, the cubic Eq. (13) is therefore generally solved numerically with a Newton-Raphson scheme. In this case, determining an adequate starting value is the main problem to address in order to design a robust and fast solution algorithm.

Efficient starting value for iterative methods
An excellent initial value for the Newton-Raphson scheme can be found by adopting the following procedure: 1. locate the local minimum closest to the largest root -if it exists, it is the extremum; to second order around that minimum; and 3. determine the greatest root of the resulting parabola and use it as a starting value.
The Taylor expansion to second order in H min , thus intersects the H axis on the right-hand side of H min at provided P CB (H min ) < 0. By completing the Taylor expansion to third order, it is straightforward to show that H 0 is greater than the root. The so-defined H 0 provides an excellent starting value not only for solving the cubic polynomial equation, but also for other iterative methods.

Carbonate, borate and water self-ionization alkalinity
With Alk CBW , C T and B T given, the equation to solve is One may either solve this equation in that rational fraction form with some iterative root-finding method or by one of the ICAC methods described above, or one may transform it into a quintic polynomial equation: with The polynomial equation can then be solved with appropriate standard root finding techniques, selecting the positive root found. Equations (15) and (14)  ] > 0. Alk CBW is probably the most commonly used approximation for total alkalinity in global carbon cycle models of all kinds of complexity. It was already adopted by Bacastow and Keeling (1973), who based their pH calculation on the quintic Eq. (15), which they solve by Newton's method, with a stopping criterion |( H )/H | < 10 −10 . Hoffert et al. (1979) adopt the same procedure (for which they refer to Keeling, 1973 andBacastow andKeeling, 1973), but with a less stringent stopping criterion |( H )/H | < 10 −6 . Keeling (1973) uses a variant, where C T is replaced by an equivalent term in pCO 2 .
As already mentioned above, LOVECLIM (Goosse et al., 2010) and HadOCC (Palmer and Totterdell, 2001) use Alk CBW as an approximation for total alkalinity. Alk CBW is also used in the PISCES model (Aumont and Bopp, 2006), following a simplified version of the OCMIP standard protocol (see next section). PISCES is included in NEMO and in some versions of the Bern3D model (Gangstø et al., 2011). Other models that base their pH calculation on Alk CBW include the Hamburg Model of the Ocean Carbon Cycle (HaMOCC) family (Maier-Reimer and Hasselmann, 1987;Heinze et al., 1991;Maier-Reimer, 1993;Maier-Reimer et al., 2005), the models of Bolin et al. (1983) and Shaffer et al. (2008). No details regarding the adopted solution algorithms are provided, though.

More complete approximations: rational function based solvers
When additional components in total alkalinity need to be considered besides carbonate, borate and water selfionization, converting the resulting rational function equation to an equivalent polynomial form becomes more and more tedious and the rational function form becomes the preferred basis for finding the solution. ICAC methods are the only ones that we have encountered so far that could pos-sibly be used to address this problem.  (1989). All of the models that participated in OCMIP had to use the provided routines for a set of well defined experiments. A number of models still routinely use these OCMIP routines for their pH calculations. These include some versions of the Bern3D model (Müller et al., 2008) and the NCAR global coupled carbon cycle-climate model CSM1.4-carbon (Doney et al., 2006). As mentioned above, PISCES (Aumont and Bopp, 2006) includes a version of the OCMIP solver trimmed down to Alk CBW only. Other models still offer the OCMIP solvers as an option. Luff et al. (2001) have provided a suite of pH calculation routines mainly meant to be used in reactive transport models, but suitable for general speciation calculations as well. The methods proposed by Luff et al. (2001) solve the complete system of equations that control the chemical equilibria between the individual species considered in the total alkalinity approximation. These are required for grid-based reactive transport models where different species are diffusing at different diffusivities. For common applications in biogeochemical carbon cycle models, this approach is nevertheless unnecessarily complex.

Other approaches
There are still some other fine pH solvers, such as CO2SYS of Lewis and Wallace (1998) and derivatives (spreadsheet versions, MATLAB versions, etc. -see http: //cdiac.ornl.gov/oceans/co2rprt.html for more information), the MATLAB routines from Zeebe and Wolf-Gladrow (2001), or the R packages seacarb (Lavigne and Gattuso, 2012) and AquaEnv (Hofmann et al., 2010(Hofmann et al., , 2012. These are, however, generally not suitable for inclusion in global biogeochemical models, as they were developed with special programming environments in mind. Their focus is more on data processing or modelling with the special programming environment they were designed for. As their names already suggest, they are mainly aimed at carbonate speciation calculations. They also often offer the possibility to chose any two among pH, [CO 2 ] (or pCO 2 ), [HCO − 3 ], [CO 2− 2 ], C T , or Alk T to calculate all the others.

pH-scale considerations
As shortly mentioned above, there are a few subtleties related to pH scales that still need to be clarified. The mere existence of more than one pH scale reflects the difficulties to apply the fundamental definition of pH (which involves an immeasurable quantity -see next section) to the experimental determination of acidity in seawater. All of our calculations nevertheless rely on the availability of equilibrium constants that have to be experimentally derived and we therefore have to care about differences arising from the usage of various pH scales.
Let us, similarly to Bates and Culberson (1977) consider the equilibrium relationship (the mass action law) for an acid dissociation reaction. Without loss of generality, we may write that relationship for the first dissociation reaction of our generic acid from Sect. 2.2.4: . When the dependency of K 1 on temperature and salinity is experimentally determined, the fraction [H n−1 A − ]/[H n A] is measured or calculated for each experiment. [H + ] cannot be directly measured, but gets assigned a value from some pH measurement via the reverse relationship [H + ]=10 −pH . Taking the negative logarithm (antilogarithm) of the previous equation and writing pK 1 = − log K 1 , we get . In a given setting (i.e. for given temperature, salinity, pressure, solution chemistry, etc.), the ratio [H n−1 A − ]/[H n A] is set and different calibrations of the pH-meter used, i.e. different scales chosen for the pH-meter, will thus produce different pK 1 values. Any experimentally derived parameterization for K 1 can therefore only be used in conjunction with a H + concentration scale that is consistent with the pH scale that was used to derive it. Before a particular empirical parameterization for K 1 can be used with a different scale of pH (e.g. due to a different conventional or operational definition of pH), it must be converted.
Additional conversion may be required because of the usage of different concentration units: both molal units (mol/kg-H 2 O) and mol/kg-seawater are common. They can be converted according to [mol/kg-SW] = m(mol/kg-H 2 O) × (1 − 0.001005S) (Dickson et al., 2007, chapter 5, p. 13), where S denotes salinity. For S = 35, the difference between the values is about ±3.5%; in log units, the values differ by ±0.016.
There is abundant literature on pH scales for seawater. Besides the original fundamental papers by, e.g. Hansson (1973), Bates and Culberson (1977), Khoo et al. (1977), Dickson and Riley (1979), Bates (1982) or Dickson (1990), the classical review papers by Dickson (1984Dickson ( , 1993, or standard textbooks (e.g. Zeebe and Wolf-Gladrow, 2001), there are also several recent papers on the subject, such as the reviews by Dickson (2010) and Marion et al. (2011) or the research paper by Waters and Millero (2013). In the following sections, we will therefore only give a comparatively general overview, which we have nevertheless tried to keep as selfconsistent as possible.

Fundamental definition of pH and standard potentiometric determination of pH
While pH as a measure of the acidity of a solution may appear as a straightforward concept, its experimental determination and interpretation are not. The fundamental definition of pH recommended by the International Union of Pure and Applied Chemistry (IUPAC, Buck et al., 2002) states that pH=− log(a H + ), where a H + denotes the activity of the H + ions in solution. a H + is related to the concentration of H + through the activity coefficient γ H + , such that a H + = γ H + [H + ]. The activity coefficient of H + depends on the exact chemical composition of the solution. The more dilute the solution is, the closer the values of activity coefficients come to one. The activity of an individual ion in solution cannot be measured by any thermodynamically valid method and the measurement of pH therefore requires an operational convention (Buck et al., 2002). The reasons for the existence of several pH scales in seawater then also simply "[. . . ] reflect the gradual gradual refinement of the experimentally convenient potentiometric determination of acidity in order that the numbers obtained might be usefully interpreted as a property of hydrogen ion in solution" (Dickson, 1984).
The potentiometric method mentioned by Dickson (1984) is the classical method used for the quantitative determination of acidity in an aqueous solution. It is based upon the use of electrochemical cells and has been used for more than 100 yr. Potentiometric pH measuring devices for practical use are made up by two electrodes: an H + sensitive glasselectrode and a well reproducible second electrode, a socalled reference electrode. Both electrodes are immersed into the sample solution to form an electrochemical cell. The potential difference between the two electrodes, i.e. the emf (electromotive force) of the cell, is linked to the logarithm of the activity of the H + ions in solution. The total emf of the cell, E, can be separated into three major contributions (Dickson, 1984). The first one, which we denote here as E gem , is due to the potential difference across the membrane of the glass electrode, which is assumed to behave following Nernst's law, i.e.

E gem = (RT ln 10/F ) log(a H + ),
where R is the universal gas constant, T the absolute temperature and F Faraday's constant. The second contribution, E J , is due to the potential difference across the liquid junction that is required to bring the filling solution of the reference electrode into contact with the sample solution. The third one, E • , takes into account all of the other potential differences arising from the characteristics of the internal electrolytes and the design of the two cells, and it can be assumed that this is an invariant of the instrument. The measured emf is thus such that By using this cell for sequentially measuring the emf of the standard buffer solution S, E S , and of the sample solution X, E X , (both at the same temperature T ) we have log(a H + (X))= log(a H + (S)) + E S −E X RT ln 10/F + E JX −E JS RT ln 10/F . This property is used (Buck et al., 2002) to operationally define the pH of the sample X from its deviation from the known or assigned pH of the standard buffer solution S, as thus implicitly assuming that the residual liquid junction potential E JX − E JS can be neglected. Buck et al. (2002) propose a number of primary standards (buffer solutions) that have to meet some fundamental metrological qualities. These primary standards are commonly known as NBS buffers (where NBS stands for the US National Bureau of Standards, now National Institute of Standards and Technology, NIST). One of the constraints imposed upon these standards is that their ionic strengths must not exceed 0.1 molal: the calibration procedure of the standards rests to some extent on the Debye-Hückel theory for ionic interactions, which is applicable only for ionic strengths < 0.1 molal (Buck et al., 2002).

Complications and simplifications for seawater pH
Seawater has a much higher ionic strength of 0.7 molal (for S = 35) than the standard NBS buffers. The use of such dilute buffers for the determination of pH in seawater samples is therefore problematic: the concentration gradients across the reference electrode's liquid junction will be significantly different between any seawater sample and the standard buffers. Hansson (1973) therefore developed and calibrated buffers on the basis of artificial seawater. In addition to this, Hansson (1973) also devised new pH scales for seawater. His scales are based upon the peculiar composition of natural seawater. Of all the solutes present, six alone make up 99.3% of the total dissolved salts: Cl − 55.0%, Na + 30.7%, SO 2− 4 7.7%, Mg 2+ 3.6%, Ca 2+ 1.2% and K + 1.1% (Millero et al., 2008); and their respective ratios remain essentially constant throughout the oceans. All of the minor constituents, which include carbonate, bicarbonate, all the nutrient salts, H + itself, etc., make up less than 0.7%. Basically, natural seawater can thus be seen as a dilute solution of the minor constituents in a solvent that is an ionic medium of rather high ionic strength, but of constant composition for a given salinity, instead of a concentrated solution of all the solutes in the pure water solvent.
As a result, the activity coefficient of H + remains fairly constant over a large concentration range (and even close to one, since the solution is dilute with the ionic medium as the solvent, Hansson, 1973). Accordingly, the concentrations and activities of H + become equivalent and it is possible to define special pH scales for seawater that are directly expressed in terms of concentrations. Solutions of well-known concentrations of H + can then easily be used as standard buffer solutions.

Towards the definition of seawater pH scales
Similarly to dilute aqueous solutions, where it is impossible to distinguish between the free and the different hydrated forms of dissolved H + since water activity is not noticeably changed during experiments, in ionic solutions where the activities of medium ions are constant, it is not possible to distinguish between the free or hydrated H + ions and those formed by protonation of medium ions (Hansson, 1973). In the seawater, H + ions may interact with SO 2− 4 and, to a lesser extent, F − ions, both of which are present in the solvent ionic medium: Accordingly, Hansson (1973) In order to make [H + ] SWS strictly proportional to [H + ] f , Dickson and Riley (1979) slightly modified Hansson's (1973) original definition, thereby simplifying the conversions between different scales and leading to the currently adopted definitions of the three most important pH scales in seawater (given here with their current common denominations and notations): the free scale (pH f =− log([H + ] f )), the total scale (pH T =− log([H + ] T )) and the seawater scale  (Hansson) are negligible in common present-day seawater (Munhoven, 1997;Dickson, 2010), those between pH f and pH T or pH SWS are not. At S = 35 and T = 298.15 K pH T and pH SWS are respectively about 0.11 and 0.12 lower than pH f (Zeebe and Wolf-Gladrow, 2001), leading to differences of the order of 30% between the corresponding concentrations of H + .

Implications for the alkalinity-pH equations
In all of the alkalinity-pH equations and equilibrium relationships from Sects. 2 and 3, as well as in the general form of the total alkalinity-pH below, Eq. (21) , since our aim here is to set up a generally valid algorithm. Instead, we consider for the time being that the effects of HSO − 4 and of HF on total alkalinity are taken into account together with the components of all the other acid systems in the sample, with K HSO 4 and K HF being expressed on the same pH scale as the constants for those components.

Development of a universal and robust algorithm
Our ultimate goal here is to develop a universal algorithm to solve the equation where collects the contributions from all the acid systems to total alkalinity, system by system, proton donors and acceptors combined for each one -except for the contribution that results from the self-ionization of water which we keep explicit in Eq. (21) -and where Alk T and all the total concentrations [ A [i] ] are given. We will use a hybrid method that combines the speed of convergence of super-linear and higher-order methods (such as the secant or the Newton-Raphson methods) with the global convergence security of the bisection or the regula falsi method. A similar method is used by the OCMIP carbonate speciation routine. Such methods are standard in root finding for non-linear equations (e.g. Dowell and Jarrett, 1971;Anderson and Björk, 1973;Bus and Dekker, 1975). They are particularly suitable for our problem with its advantageous mathematical characteristics, the more so since that problem also has intrinsic a priori root bracketing, as we will show in the next section. Because of the strict monotonicity of the rational function it is sufficient to make sure that iterates remain within bounds. As long as iterates remain within bounds, they will unconditionally improve either one of the two bounds, thus allowing us to tighten the bracketing interval at each step. We can therefore use a high-order numerical root-finding method, such as Newton-Raphson or the secant method as the main iterative scheme. In case the main scheme yields an out-of-bounds iterate at some step, that iterate is rejected and a bisection iterate is used instead. Similarly, if an iteration with the main scheme does not make the absolute value of the equation decrease faster than expected for a bisection step (i.e. by a factor of two) it is replaced by a bisection step. This helps to prevent cyclic iterations.
A bisection step may temporarily slow down the rate of convergence, but it will secure convergence and again unconditionally improve the bracketing. A regula falsi step could be used instead of bisection; bisection has proven to be more economic though.
Geosci. Model Dev., 6, 1367-1388, 2013 www. ] > 0 and that it has the infimum Alk nWinf = − i m i [ A [i] ]. It is therefore sufficient to have H inf such that Equivalently, we require that H 2 inf +s(Alk T −Alk nWinf )H inf − sK W = 0, and H inf > 0. This problem has the unique solution which would lead to R T (H sup ) = Alk nW (H sup )−Alk nWsup < 0 as requested. Equivalently, we require that H 2 sup +s(Alk T − Alk nWsup )H sup − sK W = 0 and H sup > 0. We finally obtain where sup = s 2 (Alk T − Alk nWsup ) 2 + 4sK W > 0. H inf and H sup define a universal bracketing interval for the root of Eq. (21). They only require information that can be directly derived from the nature of the acid systems considered for Alk T , and they can theoretically be used with any numerical scheme to keep iterations bracketed right from the start of the scheme, without any need for manually prescribing them.

Outline of the algorithm
The proposed algorithm is formally set up in the pH-Alk T space. There are several advantages for rooting the algorithm in the pH-Alk T space: (1) the equation's overall appearance is closer to linear in the pH-Alk T space than in the more commonly used [H + ]-Alk T space; (2) physically meaningless negative [H + ] values cannot be produced by the iterative scheme; this is not warranted with methods that are rooted in the [H + ]-Alk T space. There is nevertheless also a potential disadvantage: passing between the two spaces a priori requires costly power and logarithm evaluations at each step. As shown below, these operations can, however, be avoided to a large extent by transposing the actual calculations into the [H + ]-Alk T space and carrying them out there.
The algorithm comes in two variants: one based upon the Newton-Raphson and bisection methods, and one that is based upon the secant and bisection methods. We will first describe the Newton-Raphson/bisection variant.
Let R = R(H ) denote the rational function chosen to approximate Alk T . Before starting we determine the intrinsic lower and upper bounds H inf and H sup , and constrain the initial value H 0 to be within bounds.
Then, at each step k + 1, k = 0, . . .: 1. prepare to carry out a Newton-Raphson step where pH k+1 =pH k + pH, with pH=−R| pH k /(dR/dpH)| pH k : (dR/dpH)| pH k can be calculated from R(H k ) and dR/dH | H k , noticing that (dR/dpH)| pH k =−(dR/dH )| H k × H k × ln(10); 2. adapt the bracketing interval: if R| pH k > 0 then adjust H inf := H k , if R| pH k < 0 then adjust H sup := H k ; 3. require |R(H k )| to decrease faster than one would typically expect from bisection under the same conditions: compare it with min(|R(H j )|, ∀j < k) and if greater than half that value (bisection halves the bracketing interval at each step and is linearly convergent), do not complete the Newton-Raphson step, but adopt a bisection iterate between the current pH inf and pH sup (updated just before) and return to stage 1 for the next step; 4. provisionally set H k+1 = 10 −pH k+1 = H k × exp(−R(H k )/(H k dR/dH | H k )) to complete the Newton-Raphson step; 5. constrain H k+1 to remain within the current bracketing interval: if H k+1 > H sup or H k+1 < H inf , reject the Newton-Raphson iterate, replace it by the bisection iterate as in stage 4 and return to stage 1 for the next step; 6. stop the iterations if either the maximum permissible number of iterations is exceeded or if |(H k+1 − H k )/H k | < ( being a pre-set tolerance), else return to stage 1 for another step.
At most, one exponential has to be evaluated per step (at stage 4). This is, computationally speaking, the most expensive operation in each step. It can, however, often be avoided: when |R T (H k )/(H k dR/dH | H k )| , then H k × exp(−R T (H k )/(H k dR/dH )| H k )) H k − R T (H k )/(dR/dH | H k ) and the iterate can be assimilated to a plain [H + ]-Alk T space Newton-Raphson iterate. Once the argument of the exponential becomes sufficiently small (a threshold value of 1 in absolute value has proven efficient) we may switch to the linear approximation, thereby saving the exponential operation. The relative error |(H k+1 − H k )/H k | from the stopping criterion can be approximated by |R(H k )/(H k dR/dH | H k )|, i.e. we can reuse the argument of the exponential above (no extra operations required). At any stage, bisection between pH inf and pH sup translates to calculating H k+1 as the geometric mean of H inf and H sup : H k+1 = H inf H sup . By construction, any accepted iterate will thus be strictly between the current H inf and H sup , and, because of the strictly decreasing nature of R(H ), will always lead to contribute to improve either the lower or the upper bound.
In a variant of the above we replace the Newton-Raphson scheme by a secant scheme. However, rooting a secant scheme in the pH-Alk T space while carrying out operations in the [H + ]-Alk T space will require non-integer power operations at each step which are even more costly than exponentials. It is therefore preferable to completely root the calculations in the [H + ]-Alk T space with the secant method, despite the potential trade-offs for the convergence efficiency. Secant iterations have the advantage of requiring only one evaluation of the equation at each iteration; in addition to the equation evaluation, Newton-Raphson iterations also require the evaluation of the derivative of the equation. The cheaper iterations of the secant method possibly outweigh its lower order of convergence, which is 1.62, compared to 2 for the Newton-Raphson method.

Discussion: comparison with the OCMIP solver
A similar technique was also used in the drtsafe routine in the OCMIP suite (Orr et al., 2000), which is fundamentally the rtsafe routine of Press et al. (1989) with some essential error trapping removed. That routine also combines the global convergence properties of the bisection with the speed of convergence of the Newton-Raphson method.
The algorithm presented here differs from that used in drtsafe in several significant ways. (1) drtsafe iterations are rooted in the [H + ]-Alk T space. (2) drtsafe requires brackets to be explicitly provided at the subroutine call. In case these are inappropriate (e.g. no sign change of the equation function over the defined bracketing interval), it would simply iterate to the maximum number of iterations allowed because of the absence of validity checks and return some meaningless result (in general one of the two bounds provided). (3) drtsafe always starts its iterations from the midpoint of the provided bracketing interval. It is thus critically dependent on a valid interval (no validity checks performed though) and, because of the rooting in the [H + ]-Alk T space, on a tight bracketing interval for efficiency. The algorithms proposed here only use the bracketing values to secure convergence in case Newton-Raphson iterates are not decreasing fast enough or would go out of bounds, and they use an independent initial value. Because we root our iterations in the pH-Alk T space, even bisection steps may accommodate [H + ] changes over several orders of magnitude during the initial steps in case a far off starting value is used.

Sample implementation of the new algorithms in Fortran 90
A sample implementation of the algorithms realized in standard Fortran 90 is provided in the Supplement to this paper. Together with the drivers that were used to carry out the experiments described below they make up the Solver Suite for Alkalinity-PH Equations (SolveSAPHE). Parts of the code contain C-preprocessor directives for enabling or disabling specific parts (debugging messages, optional code parts, etc.), and to select among the cases treated below. After pre-processing, the source files are strictly standard conforming Fortran 90. The codes are made available under the GNU Lesser Public Licence, version 3. A complete user manual that covers the technical details of SolveSAPHE is included in the archive provided in the Supplement. Here we only give a short overview of the two central modules.

Summary description
The module mod chemconst provides parametric expressions for the stoichiometric constants of the acid systems taken into account (carbonates, borates, hydrogen sulphate, sulphides, phosphates, etc.). The module also hosts the j values (Eq. 5) for the various acid systems.
The module mod phsolvers provides six different solvers: 1. the function solve at general: the new algorithm described above; 2. the function solve at icacfp: a fixed-point only ICAC method; 3. the function solve at bacastow: Bacastow's method, an ICAC method with secant iterations (with secant iterations either on [H + ] or on its scaled inverse); 4. the function solve at general sec: the variant of solve at general that uses secant instead of Newton-Raphson iterations; 5. the function solve at ocmip: a re-implementation of the OCMIP solver with Newton-Raphson/bisection iterations, completed with a bare minimum of error trapping and fitted with the optional initialization scheme common to all of the solvers (the latter was only adapted to use an interval of ±0.5 pH unit interval around an optional initial value to emulate the recommended OCMIP set-up after start-up); 6. the function solve at fast: a simplified version of solve at general without the bracketing control (may not always converge).
Each one of the six solvers takes into account all of the constituents that explicitly appear in Eq. (1), except for S 2− whose concentration is negligible even at high pH values. mod phsolvers logging.F90 is a special version of mod phsolvers.F90 that includes extra bookkeeping regarding the number of iterations required for convergence, the numbers of bisection iterations due to limiting, the initial values adopted, the initial bracketing values (if relevant), the intermediate iterates, etc. mod phsolvers logging.F90 does not include solve at fast though. For more technical details, please refer to the manual that goes with the source codes. The modules mod acb solvers and mod acbw solvers provide simpler and more streamlined solvers, based upon the Alk CB and Alk CBW approximations, respectively. In these, we assume, e.g. that

Test case definitions
The list of species considered in the approximation of Alk T adopted in the six solvers presented here are much more complete than absolutely necessary in global biogeochemical cycle models. Most of these call upon the Alk CBW approximation for Alk T , which requires only Alk CBW , C T and B T to be known, the two former being generally controlled by conservation equations, the latter taken proportional to salinity. A few models, essentially those using the complete OCMIP solver, also take into account the small contributions of phosphate alkalinity ( For our test cases, we focus on the influence of various combinations of C T and Alk T on the convergence properties of the six solvers. For all the other acid dissociation systems that are taken into account in the solver routines, we adopt constant total concentrations (see below).
The pH SWS scale was adopted for all of the calculations. With each method, iterations were stopped once the relative change of an iterate compared to its predecessor fell below 10 −8 ; the maximum number of iterations was set to 50. For all of the three cases, we adopt a temperature of 275.15 K, a salinity of 35 and an applied pressure of 0 bar. Additional results for a temperature of 298.15 K or a pressure of 300 bar can be found in the Supplement.
The convergence properties for each one of the pH solvers are explored for three different (nested) subsets of the C T -Alk T space: SW1 -for C T ranging between 1.85 mmol kg −1 and 2.45 mmol kg −1 , and Alk T between 2.20 meq kg −1 and 2.50 meq kg −1 , on a regular 600 × 300 cell centred grid; SW2 -for C T ranging between 1.85 mmol kg −1 and 3.35 mmol kg −1 , and Alk T between 2.20 meq kg −1 and 3.50 meq kg −1 , on a regular 1500 × 1300 cell centred grid; 3 SW3 -for C T ranging between 0 mmol kg −1 and 6 mmol kg −1 , and Alk T between −1 meq kg −1 and 5 meq kg −1 , on a regular 600 × 600 cell centred grid.
The other concentrations are set as follows: P T = 0.5 µmol kg −1 , Si T = 5 µmol kg −1 , [NH 4 ] T = 0 µmol kg −1 and [H 2 S] T = 0 µmol kg −1 . With each test case, three different schemes are considered to start the iterations: (cub) starting values are derived from the scheme designed for the cubic polynomial in Sect. 3.2.2; (pH8) a constant starting value [H + ]=10 −8 mol kg −1 is used, except for solve at ocmip, for which the recommended "cold-start" brackets corresponding to pH = 6 and pH = 9 are used; (safe) the midpoint of the pH interval defined by the intrinsic brackets H inf and H sup (from Sect. 5.1) is used as a starting value, except for solve at ocmip, for which H inf and H sup are used as initial brackets. Timing information is based upon driver at general.F90, all other information (numbers of iterations, of divergences, errors, etc.) upon driver at logging.F90.
SW1 covers the typical range of present-day seawater samples. Every solver should be able to determine the root of the equation without a single failure. SW2 covers the expected range of sea-water samples under the S750 stabilization scenario over the next 50 000 yr (derived from simulation experiments with the coupled carbon cycle-sediment model MBM-MEDUSA, Munhoven, 2007Munhoven, , 2009. SW3 is of more theoretical nature and is meant to analyse the performance of the solvers with extremely low alkalinity or C T values. It will nevertheless also provide important information about the convergence domains of solve at icacfp and solve at bacastow, as we will see below.

Results
We here only present results obtained with the solvers from mod phsolvers (for the timings) and mod phsolvers logging. To simplify the presentation, we leave out the prefix "solve at " when referring to the different solver functions below. The testing platform had a Debian 6.0.6 operating system (32 bit kernel 2.6.32-5-686bigmem) running on a 2.53 GHz Intel Core2Duo T9400 CPU; all of the source codes were compiled with gfortran 4.4.5, without any optimization flags set. Table 1. Execution times for the SW1, SW2 and SW3 tests, each one with the three initialization types (see text). Crosses (×) indicate test series that were affected by divergences and could not be considered for time measurements (notice one exception); dashes (-) indicate that the experiment was not carried out. Within each one of the groups SW1, SW2 and SW3, figures were normalized to the execution time of the respective "cub" run with general sec and rounded to the nearest multiple of 0.05 (i.e. the order of the the largest standard deviation). "cub" results for general sec are reported in italics as they have been implicitly set to exactly 1 by the normalization procedure and are not affected by the rounding procedure.  Figure 1 shows the distributions of pH for test cases SW1, SW2 and SW3, as calculated by the new algorithm with Newton-Raphson iterations. Also shown is the equation residual for SW3 (which encompasses the two others). Residual values smaller than 10 −21 mol kg −1 in absolute value are shown as 10 −21 mol kg −1 . The residual is at least five orders of magnitude smaller than the actual H + concentrations, emphasizing that convergence was significantly reached. Execution times for the SW1, SW2 and SW3 test series are reported in Table 1. The times for each of the three test series have been normalized to the execution time of the general sec routine with cubic initialization (shown in italics). general sec was the fastest of the routines that successfully passed the complete test series.

Comparison of the six solvers
For test cases SW1 and SW2, Bacastow's method is clearly the fastest. It is about 20 % faster than the general and general sec routines developed here, and twice or two and a half times as fast as its closest relative icacfp. general sec is generally about 5-10 % faster than general. The results obtained with fast indicate that the overhead required by the safeguard bracketing requires about 5-10 % extra computing time, if everything works fine. However, the comparison of the SW2-pH8 results for general and fast shows that replacing unacceptable Newton-Raphson steps by safe but inherently slower bisection steps may overall even lead to a gain of time in more critical situations. Neither fast, nor bacastow, nor icacfp were able to complete test case SW3; ocmip was furthermore not able to complete the SW2-pH8 test, because of invalid initial bracketing over parts of the domain. Convergence failures with ocmip can be avoided if we use the intrinsic bracketing bounds obtained in Sect. 5.1, as can be seen from the "safe" initialization procedure. However, with this safe initialization procedure, the execution for ocmip exceed those of general sec with the cubic initialization by a factor of 5.6-5.7 in test case SW1 and SW2, and by a factor of 3.7 for case SW3. As can be seen in Fig. 1, SW3 includes a large number of C T -Alk T combinations that lead to extreme pH values (either lower than 4 or higher than 10), where the intrinsic bounds are comparatively restrictive.
In each test, and with every method used, the initialization procedure developed above for the cubic polynomial leads to 30-60 % shorter execution times than the constant initialization ("pH8") which may even lead to divergence (e.g. ocmip).
The reasons for the strong performance loss of icacfp in SW1 become obvious in Fig. 2. The figure shows the number of iterations required to trigger the stopping criterion for the SW2 test. The SW1 domain is included at the lower left of the SW2 domain: it ends at C T = 2.45 mmol kg −1 and Alk T = 2.50 meq kg −1 . In that area, general and bacastow require at most four iterations, ocmip generally six or seven, but icacfp often fifteen and more.

Shortcomings of ICAC methods
As we have seen, ICAC methods have divergence problems on the SW3 grid. These problems are inherent to the method and can only be alleviated to a limited extent. It should be noticed that the fixed-point equation H = Q(H ) (see Eq. 11) has a solution, i.e. that Q(H ) has a fixed-point, for any C T − Alk T pair, since the total alkalinity-pH equation has a solution.
The divergence pattern of the ICAC methods can easily be explained from the derivative of the underlying function Q with respect to H , shown for SW3 in Fig. 3. The values were calculated from the H + concentrations shown in Fig. 1     smaller than 1 in absolute value, i.e. is strictly greater than −1 here. The thick white line indicates where the derivative of Q is equal to −1. The white areas on the other graphs indicate where the methods did not converge. The white areas for icacfp clearly match the areas where the derivative of Q is lower than −1, and they are even somewhat larger.
When the derivative of Q is just slightly greater than −1, iterations may become "operationally divergent": the pre-set maximum number of iterations is insufficient to meet the stopping criterion as the generated suite converges too slowly there. Bacastow's method, on the other hand, has a slightly larger convergence domain than delimited by the thick white line in the graph of the derivative. The fixed-point equation can indeed be solved by the secant method in some instances where straight fixed-point iterations would produce slowly divergent iteration suites. As the derivative of Q is negative, fixed-point iterations oscillate around the root as long as they can be evaluated, i.e. as long as the Alk C estimates obtained from Alk T with the H iterates range between 0 and 2C T . If they can be successfully calculated, the two first iterates used to initialize the secant iterations in Bacastow's method thus bracket the root and we may expect that the first application of the secant method provides an excellent estimate for the root, even if the two first iterates would generate a fixed-point suite that slowly diverges. This is especially obvious inside the white bulge in Fig. 3c at low C T values and Alk T values greater than C T .

Quality of the cubic equation based initialization
As shown above, the initialization scheme especially designed for the iterative solution of the cubic Eq. (13) by the Newton-Raphson method is highly attractive even for more complex approximations to total alkalinity than Alk CB . This is quantified in Fig. 4, exemplified by SW2 results. The relative error of H 0 , determined as outlined in Sect. 3.2.2, on the actual H + concentration (as calculated with general) is less than 7 % over the SW2 domain. In comparison, over that same domain, the approximation Alk C Alk T and usage of Eq. (8) gives rise to errors that are about ten times as large.

Conclusions
We have explored the mathematical properties of the total alkalinity-pH equation, i.e. the equation that relates [H + ] (or equivalently pH) to total alkalinity and the total concentrations of all the acid systems contributing to total alkalinity.  We have demonstrated that the rational function expression of that equation is strictly monotone. If water self-ionization is considered, the total alkalinity-pH equation has one and only one positive root, for any given value of total alkalinity and for any given non-negative values for the total concentrations of the acid system components of total alkalinity. All other roots of the equation are either negative or complex with non-zero imaginary parts. We have shown that there are intrinsic upper and lower bounds for the positive root of the equation that only depend on information from the list of included acid systems. These seemingly straightforward mathematical properties have apparently not been published before and currently available solvers do not take advantage of them. They actually enabled us to design a universal and failsafe algorithm to solve the total alkalinity-pH equation. We propose two variants, one using the Newton-Raphson, the other the secant scheme. The performances of the two algorithms (plus one simplified version without safe-guarded iterations) were compared with some common existing ones: (1) the fixed-point iterative carbonate alkalinity correction method (ICAC), a classical method recently made popular again by Follows et al. (2006); (2) the method of Bacastow (1981), which is a variant of the previous one using secant instead of fixed-point iterations and (3) the OCMIP-2 standard protocol routines (Orr et al., 2000), re-implemented here. Source code with a reference implementation of the six algorithms discussed in the text is provided in the Solver Suite for Alkalinity-PH Equations (SolveSAPHE) in the Supplement for use under the GNU Lesser General Public Licence version 3 (LGPLv3).
We have defined three test cases for a comparative analysis of the six methods: one for the typical open-ocean concentrations of total dissolved inorganic carbon and total alkalinity of the present-day ocean; another one covering the expected future distributions of these concentrations under progressing ocean acidification and subsequent dissolution of deepsea surface sedimentary carbonates; and a third one covering extremely low concentrations of dissolved inorganic carbon and total alkalinity, and even negative values for total alkalinity. Different approaches for starting iterative root-finding methods have been tested as well for their efficiency.
The two new algorithms are the only ones that successfully complete all of the tests. The same convergence security could be achieved with the OCMIP solver as well after a modification of its initialization scheme, though at the price of much longer execution times (typically by a factor of three to six). Bacastow's method is the fastest of all the tested general methods overall in the common regions of convergence. The two new algorithms are only 10-20 % slower than Bacastow's method and more than 50 % faster than the next best performant ones. The secant variant of our algorithm is about 5-10 % faster than the Newton-Raphson variant. We have developed an original starting scheme for solving the cubic polynomial equation that is to be solved to determine pH from carbonate and borate alkalinity alone. That starting scheme can easily be completed for usage with general total alkalinity-pH equation solvers and we show here that it typically allows one to save 30-60 % of calculation time compared to the standard pH = 8 initialization.
The two proposed algorithms are furthermore extremely robust. As documented in the Supplement, the sample implementation has been successfully used with random values (covering up to six orders of magnitude) for the total concentrations of the acid system components to total alkalinity and of total alkalinity itself, with random pH starting values between 1 and 14, and still ensured convergence in 100 % of the cases.
The two new algorithms proposed thus practically offer convergence security over an extremely wide range of total concentrations for the contributions of the various acid systems to total alkalinity, while only at a marginal additional computational cost.  The result now follows from Lagrange's identity: x j y j = D 1 .
To conclude, it is then sufficient to notice that which is strictly positive if n > 0. Alternatively, the result also follows from the Cauchy-Schwarz inequality in n + 1 dimensions, noticing that the conditions for equality are not met. Accordingly, Alk A ([H + ]) is strictly decreasing as a function of [H + ].