The mathematics of the total alkalinity – p H equation : pathway to robust and universal solution algorithms

Introduction Conclusions References

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 system -dependent 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.
The realistic modelling of biologically mediated fluxes (e.g.marine primary or export production) requires the co-limitation 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.
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, that 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.Figures

Back Close
Full 2 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.

Total alkalinity: 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.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.
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 2092 Figures

Back Close
Full  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.

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, that 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 : Upon substitution of the concentrations of the species by their fractional expressions where C T is the total concentration of dissolved inorganic carbon ), K 1 and K 2 are the first and second dissociation constant for carbonic acid, we get For constant C T , the right-hand side is a strictly decreasing function of

Carbonate and borate alkalinity
The second most important component of natural present-day 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.

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 2095 Introduction

Conclusions References
Tables Figures

Back Close
Full reactions, characterised 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.
If we denote the total concentration of dissolved acid the fractions of undissociated acid and of the various dissociated forms H n−1 A − , Introduction

Conclusions References
Tables Figures

Back Close
Full  where m is an integer constant, which is dependant on the so-called 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 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.
From the previous expressions for the species fractions, we then find that where we have defined  There are two corollaries of this monotonic behaviour worth emphasizing: 1.For any acid system H n A-. . .-A n− , Alk A is bounded: it has a supremum which is equal to (n−m)[ΣA] (i.e. the limit for [H + ] → 0, not actually reachable though), and an infimum, which is equal to −m[ΣA] (i.e. its limit for [H + ] → +∞, also not actually reachable); both of these could, theoretically, be negative if m is sufficiently large.
2. For a water sample that contains a set of acids H n i A [i ] , (i = 1, . ..) of respective known total concentrations [ΣA [i ] ] and with zero proton levels respectively characterised by m i , the total alkalinity-pH equation has exactly one positive root [H + ], for any given value of Alk T : the sum of the respective alkalinity contributions over the set . .} 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 2099 Introduction

Conclusions References
Tables Figures

Back Close
Full 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.2pH 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 [H + ] is Following our discussion in Sect.2.2.1, Eq. ( 6) has a positive root if and only if 0 < Alk C < 2C T ; if there is a positive root, it is unique.Equation ( 6) can be directly solved after conversion to the quadratic equation where

Conclusions References
Tables Figures

Back Close
Full 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. ( 7) 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 relationship , 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. ( 6).
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 2101 Introduction

Conclusions References
Tables Figures

Back Close
Full a direct access to carbonate speciation (at least in part).That 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 stabilisation 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 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. ( 8) and ( 9) 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 Introduction

Conclusions References
Tables Figures

Back Close
Full  (no units given).The method is further used in LOVECLIM (Anne Mouchet, personal communication, 2012) with Alk CBW as an approximation for total alkalinity (Goosse et al., 2010) and most probably in still some others, that do, unfortunately not provide details about the calculation procedures adopted.Bacastow (1981) proposed a variant to improve the rate of convergence of fixedpoint iterations.That variant only uses the recurrence described above for the first two iterates.From the third iteration on, Bacastow (1981)  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 Franc ¸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 Full with Following our discussion in Sect.2.2.2, Eq. ( 11) 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. ( 12).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 arc-cosine, a cosine and a square root).When adopted, the cubic Eq. ( 12) 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; The Taylor expansion to second order in H min , thus intersects the H-axis on the righthand 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 rootfinding 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 ( 14) and ( 13) have the same unique positive root: when Eq. ( 13) is multiplied by the product of all the denominators of the fractions included -a product that does not change sign for [H + ] > 0 -to transform it into Eq.( 14) no new sign changes can be obtained for [H + ] > 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. ( 14), that 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) and Bacastow and Keeling (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 Introduction

Conclusions References
Tables Figures

Back Close
Full

More complete approximations: rational function based solvers
When additional components in total alkalinity need to be considered besides carbonate, borate and water self-ionization, 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 possibly be used to address this problem.However, they bear some potential pitfalls: despite having a solution, the underlying fixed-point equation may be difficult to solve numerically; intermediate estimates of Alk C may go out of bounds (remember that Alk C may only take values between 0 and 2C T ).ICAC methods can therefore not be guaranteed to find the solution.
The only commonly used carbonate chemistry routine that directly solves the rational function form of the equation is that from the Ocean Carbon Cycle Model Intercomparison Project (OCMIP).For the purpose of that project, Orr et al. (2000) prepared standard carbonate speciation routines.Total alkalinity is approximated by Alk T [HCO

pH-scale considerations
As shortly mentioned above, there are a few subtleties related to pH scales that still need to be clarified.In Eq. ( 18) above, [H + ] may be expressed on any chosen pH scale (free, total, seawater) as long as all of the stoichiometric constants K A definition, respectively proportional to the concentration on the free scale (Hansson, 1973;Dickson and Riley, 1979): where S T and F T are, respectively, the total sulphate and fluoride concentrations in solution, while K HSO 4 and K HF are the dissociation constants of HSO − 4 and of HF, respectively, here necessarily expressed on the free pH scale.Accordingly, [H + ] f can be expressed as a simple function of the adequate [H + ] in Eq. ( 18): where [H + ]=10 −pH , for a chosen pH scale, s is the corresponding scale conversion factor derived from either Eq. ( 15) or ( 16).Notice that s ≥ 1 and remind that s is independent of [H + ].Hansson's (1973)  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. ( 18) -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 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 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 2111 Figures

Back Close
Full 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.

Intrinsic bracketing bounds for the root
Our first aim is now to determine H inf > 0 such that R T (H inf ) > 0 and H sup > 0 such that R T (H sup ) < 0. We have previously established that Alk nW ([H + ]) is a strictly decreasing rational function for [H + ] > 0 and that it has the infimum Alk as in this case R T (H inf ) = Alk nW (H inf ) − Alk nWinf > 0. 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 where 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 H inf and H sup define a universal bracketing interval for the root of Eq. ( 18).They only require information that can be directly derived from the nature of the acid systems considered for Alk T .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.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  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 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 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, . . .), 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.Introduction

Conclusions References
Tables Figures

Back Close
Full

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, . . .).The module also hosts the Π j values (Eq.4) 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 initialisation 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 setup after startup); 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 Introduction

Conclusions References
Tables Figures

Back Close
Full 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.

Test case definitions
The pH SWS scale was adopted for all of the calculations below.With each method, iterations were stopped once the relative change of an iterate compared to its predecessor felt 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.15K 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; 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 instrinsic 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, . . . ) 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 stabilisation scenario over the next 50 000 yr (derived from simulation experiments with the coupled carbon cyclesediment model MBM-MEDUSA (Munhoven, 2007(Munhoven, , 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-686-bigmem)running on a 2.53 GHz Intel Core2Duo T9400 CPU; all of the source codes were compiled with gfortran 4.4.5, without any optimisation flags set.Introduction

Conclusions References
Tables Figures

Back Close
Full 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 normalised to the execution time of the general sec routine with cubic initialisation (shown in italics).general sec was the fastest of the routines that successfully passed the complete test series.
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 to 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' initialisation procedure.However, with this safe initialisation procedure, the execution for ocmip exceed those of general sec with the cubic initialisation 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 on 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 initialisation procedure developed above for the cubic polynomial leads to 30-60 % shorter execution times than the constant initialisation ("pH8") which may even lead to divergence (e.g.ocmip).
The reasons for the strong performance loss of icacfp in SW1 become obvious from Fig. 2. That 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. 10) has a solution, When the derivative of Q are 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 on Fig. 3c at low C T values and Alk T values greater than C T .

Quality of the cubic equation based initialisation
As shown above, the initialisation scheme especially designed for the iterative solution of the cubic Eq. ( 12) by the Newton-Raphson method is highly attractive even for more complex approximations to total alkalinity than Alk CB .This is quantified on 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. ( 7) gives rise to errors that are about ten times as large.Introduction

Conclusions References
Tables Figures

Back Close
Full

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 selfionization 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 fail-safe 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 safeguarded iterations) were compared with some common existing ones: (1) the fixedpoint 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) license.
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 Introduction

Conclusions References
Tables Figures

Back Close
Full 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 deep-sea surface sedimentary carbonates; 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 rootfinding 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 initialisation scheme, at the price of much longer execution times though (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 to save 30-60 % of calculation time compared to the standard pH = 8 initialisation.
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 proposed new algorithms thus offer almost convergence security over an extremely wide range of total concentrations for the contributions of the various acid systems to total alkalinity, at a marginal additional computational cost only.Introduction

Conclusions References
Tables Figures

Back Close
Full (x i y j − x j y i ) 2 . With we have n j =0 x 2 j = D, n j =0 y 2 j = D 2 and n j =0 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 Full 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 cy-Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Each proton is rather bound to a water molecule to form an H 3 O + ion, and each of these H 3 Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | 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. ( Discussion Paper | Discussion Paper | Discussion Paper | For constant B T , Alk B is again a strictly decreasing function with [H + ], similarly to Alk C .Hence, for constant C T and Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper |

.
The joint contribution of all the different dissociated and non dissociated forms of H n A to alkalinity, proton donors and proton acceptors alike, is then equal to Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | switches to a secant method to Figures Back Close Full Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | 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 fixed-point equation by the secant method.
Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | 2. develop P CB ([H+ ]) to second order around that minimum;3.determine the greatest root of the resulting parabola and use it as a starting value.FiguresBack CloseFull Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | [HF]−[H 3 PO 4 ].The different species concentrations were, as above, expressed as a function of the total concentrations of their respective acid systems and of [H + ].The resulting equation was then solved for pH by a hybrid Newtonbisection method, based upon the rtsafe solver from Press et al. (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.Discussion Paper | Discussion Paper | Discussion Paper | n i for all of the acid systems are expressed on this same pH scale.The same holds for [H + ] f , the free concentration of H + , which must also be expressed on (or converted to) this same scale.The concentrations of H + on the total and on the sea-water scale are, by Discussion Paper | Discussion Paper | Discussion Paper | original definitions of [H + ] T (and of [H + ] SWS for seawater containing fluoride) were based upon the total analytical concentration of the hydrogen ion in solution: [H + ] SWS(Hansson) =[H + ] f + [HSO − 4 ] + [HF].Differences are negligible in common present-day seawater(Munhoven, 1997).However, since our aim here is to set up a generally valid algorithm, we are not adopting the approximation [H + ] SWS =[H + ] SWS(Hansson)  a priori.Instead, we consider that the effects of Discussion Paper | Discussion Paper | Discussion Paper | 0 is the discriminant of the quadratic.Similarly, because Alk nW ([H+ ]) has the supremum Alk nWsup = i (n i − m i )[ΣA [i ] ], it is sufficient to chose H sup such that K W H sup − H sup s = Alk T −Alk nWsup Discussion Paper | Discussion Paper | Discussion Paper | 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.Discussion Paper | Discussion Paper | Discussion Paper | 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 Discussion Paper | Discussion Paper | Discussion Paper | ) drtsafe iterations are rooted in the [H + ]-Alk T space.(2) drtsafe requires Introduction Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper |

Figure 1
Figure1shows 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 Table1.The times for each of the three test series have been normalised to the execution time of the general sec routine with cubic initialisation (shown in italics).general sec was Fig. 1, SW3 includes a large number of C T -Alk T Discussion Paper | Discussion Paper | Discussion Paper | 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 on Fig. 3.The values were calculated from the H + concentrations shown on Fig. 1.The derivative values are negative over the whole domain and fixed-point iterations thus oscillate around the solution (i.e. the fixed-point of Q).Fixed-point iterations can only converge to the fixed point of Q where the derivative is strictly 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.They are even somewhat larger.Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | acid system H n A -H n−1 A − -. . .-A n− , i.e.for constant [ΣA], at equilibrium.For notational convenience, we writeAlk A = [ΣA] H + ] n−j and D 1 = n j =0 j Π j [H + ] n−j .Since Π j > 0 for any j ≥ 0, we know that D > 0 and D 1 > 0. We may then 2 Π j [H + ] n−j , Discussion Paper | Discussion Paper | Discussion Paper |which is positive, similarly to D and D 1 .We hence find Discussion Paper | Discussion Paper | Discussion Paper | Peng, T.-H., Takahashi, T., Broecker, W. S., and Olafsson, J.: Seasonal variability of carbon dioxide, nutrients and oxygen in the northern North Atlantic surface water: observations and a model, Tellus B, 39, 439-458, 1987.2100, 2103 Press, W. H., Flannery, B. P., Teukolsky, S. A., and Vetterling, W. T.: Numerical Recipes (FOR-TRAN Version), Cambridge University Press, Cambridge, 1989.2108Discussion Paper | Discussion Paper | Discussion Paper |

Fig. 1 .Fig. 2 .
Fig. 1. pH SWS values obtained with the new universal algorithm (general) for test cases (a) SW1, (c) SW2 and (b) SW3 -please notice the extended colour scale.(d) Absolute value of the equation residual at the adopted root value, derived with that same algorithm started with the cubic-based initialisation to solve test case SW3.Applied convergence criterion: |H n+1 − H n |/H n < 10 −8 .

Fig. 3 .
Fig. 3. (a) Derivative with respect to H of the function underlying the ICAC methods, i.e. of the function Q given by Eqs.(8) and (9) that defines the recurrence H n+1 = Q(Alk C (Alk T , H n ), C T ).The white line indicates where the derivative is equal to −1; in the stippled area, the derivative is strictly lower than −1.Also shown are the numbers of iterations required to meet the convergence criterion for (b) general, (c) icacfp and (d) bacastow with secant iterations on [H + ].White areas indicate no convergence or an excessive number of iterations (n > 50).
Alk C < 2C T if C T = 0.Both bounds are strict (i.e. they cannot be reached) and represent the limits of Alk C (C T ; [H + ]) for [H + ] → +∞ (lower bound) and [H + ] → 0 (upper bound), for C T fixed.
Ridgwell (2001)el (1995)ovide 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 fixed-point 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)adopts Alk CB +[OH Arndt et al. (2011)007) T and Alk C by using his own procedure described above.GENIE(Ridgwell et al., 2007)initially used the same procedure asRidgwell (2001); in more recent versions of GENIE, a complete representation of the phosphoric acid component is used(Ridgwell, personal communication, 2012).Arndt et al. (2011)use Alk CBW +[HS − ] as an approximation to total alkalinity in GEOCLIM reloaded.They continue to iterate until |Alk CBW +[HS − ]−Alk T | < 10 −6 That local minimum, if it exists (i.e. if c

3.5 Other approaches
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 byLuff 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.There are still some other fine pH solvers, such as CO2SYS of Lewis and Wal- (Lavigne and Gattuso, 2012)1)preadsheet versions, MATLAB versions, etc.-see http://cdiac.ornl.gov/oceans/co2rprt.html for more information), the MATLAB routines fromZeebe and Wolf-Gladrow (2001)or the SEACARB package for R(Lavigne and Gattuso, 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.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 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 setH k+1 = 10 −pH k+1 = H k ×exp(−R(H k )/(H k dR/dH| H k )) to completethe Newton-Raphson step.5.ConstrainH 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.
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.

Table 1 .
Execution times for the SW1, SW2 and SW3 tests, each one with the three initialisation 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.