EGU Journal Logos (RGB) O O O p p p Advances in e Annales e Nonlinear Processes e n n n A A A Geosciencesc Geophysicaec in Geophysicsc c c c e e e s s s s s s Natural Hazards O Natural Hazards O p p e e and Earth System n A and Earth System n A c c c c Sciencese Sciencese s s s s Discussions Atmospheric O Atmospheric O p p e e n n Chemistry A Chemistry A c c c c and Physicse and Physicse s s s s Discussions Atmospheric O Atmospheric O p p e e n n Measurement A Measurement A c c c c Techniquese Techniquese s s s s Discussions O O p p e Biogeosciencese n n Biogeosciences A A c c c Discussions c e e s s s s O O p Climate p Climate e e n n A of the Past A of the Pastc c c c e e s Discussionss s s O O p Earth System p Earth System e e n n A Dynamics A Dynamicsc c c c e e s Discussionss s s Geoscientific Geoscientific O O p p Instrumentation e Instrumentation e n n A A Methods andc Methods andc c c e e Data Systemsss Data Systemsss Discussions Geosci. ModelDev.,6,1367–1388,2013 O GeoscientificO p p www.geosci-model-dev.net/6/1367/2013/ Geoscientificen en doi:10.5194/gmd-6-1367-2013 Model Development Ac Model Development Ac c c ©Author(s)2013. CCAttribution3.0License. es Discussionses s s Hydrology and O Hydrology and O p p e e Earth Systemn A Earth Systemn A c c c c Sciencese Sciencese s s s s Mathematics of the total alkalinity–pH equation – Discussions O O pathway to robust and universal solution algorithms: p p e Ocean Sciencee n n Ocean Science A A the SolveSAPHE package v1.0.1 cc Discussions cc e e s s s s G.Munhoven De´partementd’Astrophysique,Ge´ophysiqueetOce´anographie,Universite´ deLie`ge,4000Lie`ge,Belgium O O p p e Solid Earthe n n Correspondenceto:G.Munhoven([email protected]) Solid Earth A A cc Discussions cc e e s s Received:21February2013–PublishedinGeosci.ModelDev.Discuss.:12March2013 s s Revised:17June2013–Accepted:8July2013–Published:30August2013 O O p p e The Cryospheree n n Abstract.Thetotalalkalinity–pHequation,whichrelatesto- 1 Introduction The Cryosphere A A c c tal alkalinity and pH for a given set of total concentrations ce Discussions ce s s of the acid–base systems that contribute to total alkalinity Biogeochemicalmodelshavebecomeindispenssabletoolsto s in a given water sample, is reviewed and its mathematical improveourunderstandingofthecyclingoftheelementsin propertiesestablished.Weprovethattheequationfunctionis theEarthsystem.Acentralandcriticalcomponentofalmost strictly monotone and always has exactly one positive root. all biogeochemical models is the pH calculation routine. In Differentcommonlyusedapproximationsarediscussedand ocean carbon cycle models, the air–sea exchange of CO2 is compared. An original method to derive appropriate initial directly linked to the surface ocean [CO2]; the preservation valuesfortheiterativesolutionofthecubicpolynomialequa- of biogenic carbonates in the surface sediments at the sea tionbaseduponcarbonate-borate-alkalinityispresented.We flooriscloselylinkedtothedeepsea[CO2−](Broeckerand 3 then review different methods that have been used to solve Peng,1982).ThefractionsofCO ,HCO− andCO2− inthe 2 3 3 the total alkalinity–pH equation, with a main focus on bio- totaldissolvedinorganiccarbon(i.e.thespeciationofthecar- geochemical models. The shortcomings and limitations of bonatesystem)arecontrolledbypH.Hence,pHchangesin thesemethodsaremadeoutanddiscussed.Wethenpresent seawatermaydirectlyinfluenceair–seaexchangeofCO or 2 twovariantsofanew,robustanduniversallyconvergental- the preservation of carbonates in the deep sea. Conversely, gorithm to solve the total alkalinity–pH equation. This al- thedissociationofacids,suchascarbonicacid,alsocontrols gorithm does not require any a priori knowledge of the so- pH: when the ocean takes up or releases CO (e.g. as a re- 2 lution. SolveSAPHE (Solver Suite for Alkalinity-PH Equa- sult of a rise or a decline of the abundance of CO in the 2 tions)providesreferenceimplementationsofseveralvariants atmosphere), its pH changes. The currently ongoing ocean ofthenewalgorithminFortran90,togetherwithnewimple- acidification due to the massive release of CO into the at- 2 mentations of other, previously published solvers. The new mosphere by human activity is but one example of such an iterative procedure is shown to converge from any starting inducedpHchange. valuetothephysicalsolution.Theextracomputationalcost The nitrogen cycle is another important biogeochemical for the convergence security is only 10–15% compared to cycle where pH plays an important role. The speciation of thefastestalgorithminourtestseries. dissolvedammoniumis–asthatofanyacid–basesystem– − dependent on pH, NH being more abundant than NH at 3 4 high pH, and less abundant at low pH. At pH>9, the con- centrationofNH inseawatermayreachtoxiclevels. 3 PublishedbyCopernicusPublicationsonbehalfoftheEuropeanGeosciencesUnion. 1368 G.Munhoven: Solvingthealkalinity–pHequation: SolveSAPHE The realistic modelling of biologically mediated fluxes the common way to resolve this underdetermination is to (e.g. marine primary or export production) requires the co- consider another conservative quantity: total alkalinity, also limitation or even inhibition by different chemical compo- called titration alkalinity. Total alkalinity, which is also an nentstobetakenintoaccount.Thenitrogenandcarboncy- experimentallymeasurablequantity,tiesallthedifferentacid cles,alreadymentionedabove,stronglyinteract,bothinthe systemspresentinawatersampletogetherandallowsusto ocean and on land. In the ocean, Fe and other metals act as solvethespeciationproblem.IncomparisontopH,ithasthe micronutrientsandonceagain,pHplaysanimportantroleas advantage of being a conservative quantity: it is only con- thesolubilityofmetalsisstronglydependentonpH(Millero trolledbyitssourcesandsinks,anditisindependentontem- et al., 2009). The resulting coupling of the biogeochemical peratureandpressure(ZeebeandWolf-Gladrow,2001). cycles of different elements makes biogeochemical models Inthefollowingsection,weprovideacomprehensivein- become more and more complex and pH calculation more troductiontotheconceptofalkalinity.Inourexplorationof andmoredifficult. themathematicalpropertiesoftheequationthatrelates[H+] Biogeochemicalmodelsarenowincreasinglyusedforset- tototalalkalinitystartwithadetailedpresentationofvarious tingsthatarestronglydifferentfrompresentday.Typicalap- approximations commonly used for present-day seawater. plications include future ocean acidification (e.g. Caldeira Theanalysisofthemathematicalpropertiesoftheseapprox- and Wickett, 2003), the Paleocene–Eocene Thermal Max- imations will provide useful hints for the characteristics of imum (e.g. Ridgwell and Schmidt, 2010), Snowball Earth thegeneralcase.InSect.3ofthispaper,wepresentsolution (e.g. Le Hir et al., 2009), etc. Some commonly used pH methods for deriving pH from each of the various approx- solvers may possibly become unstable and produce unreli- imations to total alkalinity considered. Complications that ableresults.Theconvergencepropertiesofcurrentlyusedso- might possibly arise from the various pH scales that are in lutionmethodshasactuallyneverbeensystematicallytested. useinmarinechemistryareelucidatedinSect.4.InSect.5, Unfortunately, information on pH solver failures is only wethenshowthatthereareintrinsicboundsthatbracketthe seldom published. Zeebe (2012) reports for his LOSCAR root of the total alkalinity–pH equation, and that can be di- model that negative H+ concentrations may be obtained rectlyderivedfromtheapproximationusedtorepresenttotal when starting with total alkalinity and dissolved inorganic alkalinity.Theexistenceofsuchboundsmakesitpossibleto carbon concentrations in a very high ratio, requiring the defineanew,universalalgorithmtosolvethealkalinity–pH modelruntoberestartedwiththerespectiveconcentrations equation, which requires no a priori knowledge of the root. inalowerratio.Hofmannetal.(2010)alsoindicatethatthe Areferenceimplementationoftwovariantsofthenewalgo- standard pH solving routine in their R modelling environ- rithm is presented in Sect. 6. The algorithms are tested for mentAquaEnvmayfailwhentryingtocalculatethepHfor their efficiency and robustness and their performance com- samples with very low or zero dissolved inorganic carbon pared with that of the most common previously published concentrations.Inthiscase,theyresorttoageneralpurpose generalsolutionmethods. interval based root finding routine instead, adopting a very largebracketinginterval(seeHofmannetal.,2012),possibly 2 Totalalkalinity:generaldefinitionand leadingtoaconsiderableperformanceloss.AndyRidgwell, approximations inhiseditorialcommenttothecompaniondiscussionpaper, mentions convergence problems encountered with the GE- Inthefollowingpartsofthissection,wereviewanumberof NIE model code (Ridgwell et al., 2007) encountered while aspects of total alkalinity in natural waters. The main focus studyingtheeffectofanartificialadditionoflime(CaO)to will be put onto seawater and on the carbonate system, but theoceansurface(aparticulargeoengineeringmethodmeant all the presented developments can be applied to any natu- to accelerate the uptake of carbon dioxide from the atmo- ralwatersample,providedtherequiredthermodynamiccon- sphere)oncetotalalkalinitycametoexceedthetypicalsur- stantsareknown.Webrieflyrecallthedifferentapproxima- face ocean concentrations of dissolved inorganic carbon by tions commonly used for calculating pH and the speciation aboutafactoroftwo.Aswewillshowbelow,thethreemod- of acid systems. We will then establish a few basic proper- elsuseessentiallyequivalentpHcalculationmethods,which ties of the expressions that relate the various types of alka- becomedivergentunderthosetypicalconditions. linitytototalconcentrationsandpH.Althoughsimple,these The speciation of any acid system, i.e. the determination properties do not seem to have been previously explored in of the concentrations of each one of the undissociated and detail,norexploitedfordesigningmethodsofsolutionofthe the different dissociated forms of an acid, is an underdeter- alkalinity–pHequation. minedproblemifonlythetotalconcentrationandthermody- Althoughweprimarilyfocusonmodellinginthefollow- namicorstoichiometricconstantsareknown.Thisunderde- ing developments, the calculation procedures are obviously termination can be lifted if pH is known. Being dependent alsoapplicableinexperimentalset-ups. ontemperatureandpressure,neitherpHnor[H+]are,how- ever,wellsuitableforbeingusedintransportequations,and thus in biogeochemical models. In biogeochemical models, Geosci. ModelDev.,6,1367–1388,2013 www.geosci-model-dev.net/6/1367/2013/ G.Munhoven: Solvingthealkalinity–pHequation: SolveSAPHE 1369 2.1 Totalalkalinity instead of a strong acid must be added to reach the equiv- alence pH point of 4.5. Interested readers may refer, e.g. to 2.1.1 Generaldefinition KirbyandCravottaIII(2005)andreferencesthereinforsuch –fromamarinechemist’spointofview–exoticsamples. Total alkalinity, also called titration alkalinity, denoted here AlkT,reflectstheexcessofchemicalbasesofthesolutionrel- 2.1.2 ThepH–totalalkalinityequation ativetoanarbitraryspecifiedzerolevelofprotons,orequiva- lencepoint.Ideally,Alk representstheamountofbasescon- Total alkalinity as defined above is a conservative quantity T tainedinasampleofseawaterthatwillacceptaprotonwhen withrespecttomixing,changesintemperatureandpressure the sample is titrated with a strong acid (e.g. hydrochloric (Wolf-Gladrow et al., 2007). It is therefore a cornerstone in acid)tothecarbonicacidendpoint.Thatendpointislocated biogeochemical cycle models which are most conveniently atthepHbelowwhichH+ionsgetmoreabundantinsolution formulated on the basis of conservation equations. In such thanHCO−ions;itsvalueiscloseto4.5.H+addedtowater models, definition/Eq. (1) above, or an adequate variant, is 3 atthispHbyaddingstrongacidwillremainassuchinsolu- used to solve the inverse problem for [H+]. All of the in- tion.Pleasenoticethat,forthesakeofasimplernotation,we dividual species concentrations appearing in Eq. (1) can be follow here the common usage of denoting protons in solu- expressedintermsofthetotalconcentrationsoftheacidsys- tionbyH+,althoughfreeH+ionssensustrictodoonlyexist temsthattheyrespectivelybelongtoandof[H+].Giventhe in insignificantly small amounts in aqueous solutions. Each evolutionsofthetotalconcentrationsofalltheacidsystems protonisratherboundtoawatermoleculetoformanH O+ considered (dissociated and non-dissociated forms) and of 3 ion,andeachoftheseH O+inturnisfurthermoregenerally Alk –allofwhichcanbederivedfromappropriateconser- 3 T hydrogen bonded to three other H O molecules to form an vation equations – expression (1) is interpreted as an equa- 2 H O+ion(Dickson,1984). tionfor[H+]or,equivalently,pH.Wewillthereforecallthat 9 4 Rigorously speaking, Alk is defined as the number of equationthetotalalkalinity–pHequation. T molesofH+ ionsequivalenttotheexcessof“protonaccep- WemightactuallyhavecalledourequationsimplythepH tors”,i.e.basesformedfromacidscharacterizedbyapK ≥ equation.Alk doesindeednotplayanyspecialormoreim- A T 4.5inasolutionofzeroionicstrengthat25◦C,over“proton portantrolethananyofthetotalconcentrationsoftheother donors”, i.e. acids with pK <4.5 under the same condi- acidsystemsconsidered.Wedo,however,feelthatthisname A tions,perkilogramofsample(Dickson,1981). wouldhavebeentoogeneralandthusprefertoinclude“total Withemphasisonthemostimportantcontributors,arather alkalinity”inthenametoreflectthattheoverallstructureof completeexpressionforAlk inaseawatersampleis theequationderivesfromthedefinitionoftotalalkalinity. T AlkT=[HCO−3]+2×[CO23−]+[B(OH)−4]+[OH−] 2.1.3 Typicalapplicationsinbiogeochemicalmodels +[HPO2−]+2×[PO3−]+[H SiO−] 4 4 3 4 Inatypicalglobaloceancarboncyclemodel,totalalkalinity +[NH3]+[HS−]+2×[S2−]+... maycommonlybeapproximatedby −[H+] −[HSO−]−[HF]−[H PO ]−..., (1) f 4 3 4 Alk ’[HCO−]+2×[CO2−]+[B(OH)−]+[OH−]−[H+], (2) T 3 3 4 wheretheellipsesrefertootherpotentialprotondonorsand acceptorsgenerallypresentatnegligibleconcentrationsonly. where[H+]’[H+] +[HSO−+[HF].[HCO−]and[CO2−] f 4 3 3 Alloftheconcentrationsaretotalconcentrations(whichin- can be expressed as a function of the total concentration of clude free, hydrated and complexed forms of the individ- dissolvedinorganiccarbon,C ,and[H+](seeSect.2.2.1for T ual species), except for [H+] , which only includes the free details) while [B(OH)−] can be expressed as a function of f 4 andhydratedforms.Therearealternativedefinitionsthatcan thetotalborateconcentration,B ,and[H+](seeSect.2.2.2 T be found in the literature, which lead to similar, although for details); [OH−] is directly linked to [H+] via the equi- not necessarily exactly the same, expressions. However, the librium constant for the dissociation of water. Accordingly, above definition is the one that reflects the titration proce- Eq. (2) provides a relationship between C , B , Alk and T T T dureusedtomeasurealkalinitythemostaccurately.Wewill [H+] (i.e. pH). The model provides conservation equations thereforebasethefollowingdevelopmentsuponit. forC andAlk ;B cangenerallybetakenproportionalto T T T In other natural water samples (lake, river, or brines) the salinity,whoseevolutioneitherfollowsaprescribedscenario constituent list in Eq. (1) needs to be adapted: some con- ormayalsoderivedfromaconservationequation.Relation- stituents may be neglected and bases of other acid sys- ship(2)thusreducestoanequationin[H+].Thesolutionof tems have to be included (e.g. bases derived from organic thatequationfinallyprovidesameanstocalculatethecom- acids, from dissolved metals, etc.). While total alkalinity pletespeciationofthecarbonateandtheboratesystems. in seawater samples typically ranges between about 2 and Inotherbiogeochemicalstudieswhereothersystemsthan 2.6meqkg−1,acidminedrainagesamplesmayevenpresent thecarbonatesystemareofinterest(suchasammonium,sul- negative alkalinity, representing the fact that a strong base phides, etc.), the procedure is entirely analogue. Each one www.geosci-model-dev.net/6/1367/2013/ Geosci. ModelDev.,6,1367–1388,2013 1370 G.Munhoven: Solvingthealkalinity–pHequation: SolveSAPHE oftheindividualspeciesconcentrationsthatneedtobecon- 2.2.2 Carbonateandboratealkalinity sideredinEq.(1)forthatparticularapplicationisexpressed in terms of the total concentration of the acid system that it The second most important component of natural present- belongstoandconservationequations,scenariosormeasure- day seawater alkalinity is borate alkalinity, AlkB. Together mentsthatareusedtoevaluateallofthetotalconcentrations, withthecarbonatealkalinitywehave including total alkalinity. These steps again reduce Eq. (1) into an equation in [H+], whose solution provides a direct AlkCB=AlkC+AlkB=[HCO−3]+2[CO23−]+[B(OH)−4]. meanstocalculatethespeciationsofallthesystemsconsid- Uponsubstitutionoftheindividualspeciesconcentrationsby ered. theirfractionalexpressionsasafunctionof[H+],weget 2.2 Commonapproximationsfortotalalkalinity K [H+]+2K K K inseawater Alk =C 1 1 2 +B B , CB T[H+]2+K [H+]+K K T[H+]+K 1 1 2 B Here we first analyse the forward problem for a few spe- whereB isthetotalconcentrationofdissolvedboratesand T cific approximations used for seawater: for given total con- K is the dissociation constant for boric acid. For constant B centrations of dissolved inorganic carbon, total borate, etc., B , Alk is again a strictly decreasing function with [H+], T B weanalysehowtheexpressionsforthedifferenttypesofal- similarly to Alk . Hence, for constant C and B , Alk kalinity change as a function of [H+]. This simple analysis is a strictly decreCasing function with [H+]Tand, asTa consCeB- will already provide valuable insight into the overall math- quence,0<Alk <2C +B aslongasC +B 6=0. CB T T T T ematical properties of the total alkalinity–pH equation and its subcomponents, which we can exploit later for the most generalcase. 2.2.3 Carbonate,borateandwaterself-ionization alkalinity 2.2.1 Carbonatealkalinity In a third stage, we may consider the alkalinity that The contribution of the carbonic acid system (or carbonate arises from the dissociation of the solvent water itself (by system) to total alkalinity is called carbonate alkalinity and self-ionization) in addition to carbonate and borate alkalin- wedenoteitbyAlkC: ity and get the next important approximation for natural present-day seawater, called practical alkalinity by Zeebe Alk =[HCO−]+2[CO2−]. C 3 3 andWolf-Gladrow(2001): Upon substitution of the concentrations of the species by Alk theirfractionalexpressionsasafunctionof[H+], CBW =Alk +[OH−]−[H−] CB [HCO−3]=CT[H+]2+KK1[[HH++]]+K K =[HCO−3]+2[CO23−]+[B(OH)−4]+[OH−]−[H+]. 1 1 2 Uponsubstitutionbytherespectivespeciationrelationships, and weget [CO23−]=CT[H+]2+KK11[HK+2]+K1K2, AlkCBW=CT[H+K]21+[HK+1][+H+2K]+1KK21K2 (3) K K where C is the total concentration of dissolved inorganic +B B + W −[H+], carbon (CT =[CO ]+[HCO−]+[CO2−]), K and K are T[H+]+KB [H+] T 2 3 3 1 2 the first and second dissociation constant for carbonic acid, whereK isthedissociationconstantofwaterinseawater. W weget Atthisstage,wedonotwanttoinsistonsubtletiesrelatedto pHscales.Normally,thelastterm[H+]inthetwoprevious AlkC=CT[H+K]21+[HK+1][+H+2K]+1KK21K2. edqifufeartieonncsesbheotwuledenac[tHu+al]lyanrdea[dH+[H]+i]nf.SWecet.w4ibllelaodwd.ress the f SinceAlk isdecreasingwith[H+],forconstantC and For constant C , the right-hand side is a strictly decreasing CB T functionof[H+T]:itsderivativewithrespectto[H+]isstrictly BT, the same holds for AlkCBW, because KW/[H+]−[H+] negative for positive [H+]. As a consequence, 0<Alk < is again decreasing with [H+]. However, unlike AlkCB, C 2C if C 6=0. Both bounds are strict (i.e. they cannot be AlkCBW is unbounded and it can take arbitrarily low values reaTched) Tand represent the limits of Alk (C ;[H+]) for (for[H+](cid:29))andarbitrarilygreatvalues(for[H+](cid:28)). C T [H+]→+∞ (lower bound) and [H+]→0 (upper bound), forC fixed. T Geosci. ModelDev.,6,1367–1388,2013 www.geosci-model-dev.net/6/1367/2013/ G.Munhoven: Solvingthealkalinity–pHequation: SolveSAPHE 1371 2.2.4 Contributionofagenericacidsystemtototal IfwedenotethetotalconcentrationofdissolvedacidH A n alkalinity by [6A]=[H A]+...+[An−], the fractions of undissoci- n ated acid and of the various dissociated forms Hn−1A−, In common seawater, AlkCBW is entirely sufficient even for Hn−2A2−,...,An−are applications that require high accuracy. However, in some [H A] [H+]n cases other systems than the carbonate and borate systems n = needtobeconsidered.Thisisespeciallythecaseinsuboxic [6A] [H+]n+K1[H+]n−1+K1K2[H+]n−2 +...+K1K2···Kn and anoxic waters, such as semi-closed fjords (e.g. Fram- varenFjordinNorwaystudiedbyYaoandMillero,1995)or [H+]n = , atalargerscale,theBlackSea(e.g. Dyrssen,1999),where, n j [H+]n+ P[H+]n−j QK e.g.thecontributionfromsulphidescannotbeneglected. i j=1 i=1 Inordertogeneralizeouranalysisofthetotalalkalinity– pHequation,letusconsideragenericacid,denotedbyH A, that may potentially lead to n successive dissociation renac- [Hn−1A−] = K1[H+]n−1 , [6A] n j tions,characterizedbystoichiometricdissociationconstants [H+]n+ P[H+]n−j QK i K1,K2,...,Kn,respectively: j=1 i=1 . HnA(cid:10)H++Hn−1A−, K1= [H+][[HHnA−]1A−] .. j n (QK )[H+]n−j Hn−1A−(cid:10)H++Hn−2A2−, K2= [H+[H][nH−n1−A2−A]2−] [Hn[−6jAA]j−] = [H+]ni+=1Pni[H+]n−kQk Ki, k=1 i=1 . . . . . . . . . n Q K HA(n−1)−(cid:10)H++An−, Kn= [[HHA+(]n[−A1n)−−]]. [[A6nA−]] = [H+]n+ Pni=[1H+i]n−j Qj Ki. j=1 i=1 For simplicity, we omit the “∗” superscript commonly used Thejointcontributionofallthedifferentdissociatedandnon- elsewhere to differentiate stoichiometric from thermody- dissociated forms of H A to alkalinity, proton donors and n namic dissociation constants (i.e. elsewhere stoichiometric protonacceptorsalike,isthenequalto constantsgenerallywriteK∗insteadofK ).Throughoutthis i i pacatpiveri,tiethse. Acosnssutacnht,sthueseydinwcilluldreeltahtee ceoffneccetnotrfaaticotnivsitayndconeof-t AlkA=Xn (j−m)[Hn−jAj−], ficients that differ from unity. The values of such constants j=0 not only depend on temperature and pressure but also on wheremisanintegerconstant,whichisdependantontheso- theionicstrengthofthesolution.Everythingdevelopedhere calledzeroprotonlevelofthesystemunderconsideration: furthermore applies to all kinds of acids, be they of Arrhe- nius, Brønsted–Lowry, Lewis or any other type, even if the – missuchthatpKm<4.5<pKm+1 if pK1<4.5 adopted notation could possibly suggest that our develop- andpKn>4.5 mentsonlyapplytoArrhenius-typeacids. – m=0ifpK >4.5 1 – m=nifpK <4.5 n Since pKm<4.5, all of the Hn−jAj− in the HnA - ... - An− system for j =0,...,m−1 are proton donors: the last one (j =m−1) has a strength of 1eqmol−1, the second to lastone(j =m−2)of2eqmol−1,etc.SincepKm+1>4.5, thedissociationproductsHn−jAj− forj =m+1,...,nare proton acceptors, the first one (j =m+1) with a strength of1eqmol−1,thesecondone(j =m+2)withastrengthof 2eqmol−1,etc.Forthecarbonicacidsystem,e.g.n=2and m=0; for the boric acid system, n=1 and m=0; for the phosphoricacidsystem,n=3andm=1. www.geosci-model-dev.net/6/1367/2013/ Geosci. ModelDev.,6,1367–1388,2013 1372 G.Munhoven: Solvingthealkalinity–pHequation: SolveSAPHE From the previous expressions for the species fractions, 3.1 Carbonatealkalinitybasedsolutions wethenfindthat The straight approximation Alk ’Alk is often used in n T C P(j−m)5j[H+]n−j textbooks (e.g. Broecker and Peng, 1982). There are only Alk ([H+])=[6A] j=0 a few models (e.g. Opdyke and Walker, 1992; Walker and A Pn 5 [H+]n−j Opdyke,1995)thatuseitdirectlyfortheircarbonatechem- j istry speciation. For numerical modelling purposes, its us- j=0  n  age is indeed somewhat problematic. [H+] calculated from Pj5j[H+]n−j AlkT and CT data, by assuming that AlkC=AlkT are typi- =[6A] j=0 −m, (4) cally 30–40% too low (i.e. 0.15–0.2 pH units too high) for  Pn 5 [H+]n−j  present-dayseawatersamples.Furthermore,thesensitivityof j theC -Alk systemtoperturbationsisstrongerthanthatof j=0 T C theC -Alk system:equilibriumpCO changes,e.g.are T CBW 2 wherewehavedefined oftheorderof20%larger(Munhoven,1997). j Thecalculationof[H+]fromC -Alk remainsneverthe- Y T C 5j = Ki, j =1,...,n and 50=1 (5) lessimportant,asmoreadvancedmethodssuchasthosepro- i=1 posed by Bacastow (1981), Peng et al. (1987) or Follows tosimplifythenotation. et al. (2006), where AlkC is iteratively recalculated from Similartothecarbonateandboratesystemsabove,AlkAis morecompleteapproximationstoAlkT(ICACmethods–see strictlydecreasingwith[H+],for[6A]fixed.Amathemati- below),relyonit. callyrigorousdemonstrationofthisbehaviourforthegeneral 3.1.1 Fundamentalsolution caseisprovidedinAppendixA. There are two corollaries of this monotonic behaviour For given Alk and C (C >0), the equation to solve for worthemphasizing. C T T [H+]is 1. ForanyacidsystemH A-...-An−,Alk isbounded:it n A hasasupremumwhichisequalto(n−m)[6A](i.e.the K [H+]+2K K limitfor[H+]→0,notactuallyreachablethough),and RC([H+])≡CT[H+]21+K [H+]+1 K2 K −AlkC=0. (7) 1 1 2 an infimum, which is equal to −m[6A] (i.e. its limit for [H+]→+∞, also not actually reachable); both of FollowingourdiscussioninSect.2.2.1,Eq.(7)hasapositive thesecould,theoretically,benegativeifmissufficiently rootifandonlyif0<AlkC<2CT;ifthereisapositiveroot, large. itisunique. Equation(7)canbedirectlysolvedafterconversiontothe 2. ForawatersamplethatcontainsasetofacidsHniA[i], quadraticequation: (i=1,...) of respective known total concentrations [6A[i]] and with zero proton levels respectively char- P ([H+])≡[H+]2+a [H+]+a =0, (8) C 1 0 acterizedbym ,thetotalalkalinity–pHequation, i XAlk ([H+])+ Kw −[H+] −Alk =0, (6) where A[i] [H+] f T (cid:18) (cid:19) (cid:18) (cid:19) i a =K 1− CT and a =K K 1− 2CT . hasexactlyonepositiveroot[H+],foranygivenvalue 1 1 AlkC 0 1 2 AlkC of Alk : the sum of the respective alkalinity contribu- T tions over the set {HniA[i]|i=1,...} of all the acid sys- FquoardrvaatliicdeqAulaktiConvhaalusetswo(ir.eea.lrfooorts0,a<pAoslkitCiv<ea2nCdTa),netghais- tems active in the sample is a strictly decreasing func- tion of [H+]; the contribution from the dissociation of tiveone.Thepositiverootis waterisalsostrictlydecreasingwith[H+],andmaythe- (cid:18) (cid:19) oreticallytakeanyvaluebetween+∞and−∞. [H+]=Q(Alk ,C )≡ K1 CT −1+p1 , (9) C T C 2 Alk C 3 Alkalinity–pH equation in biogeochemical models: where approximationsandmethodsofsolution (cid:18) C (cid:19)2 K (cid:18)2C (cid:19) 1 = 1− T +4 2 T −1 . (10) Inthissection,wearegoingtoreviewthemostcommonap- C Alk K Alk C 1 C proximations used in ocean carbon and biogeochemical cy- cle models, focusing on how the corresponding equation is ForAlk valuesthatareoutofrangeEq.(8)eitherhastwo C solved. negativeortwocomplexroots. Geosci. ModelDev.,6,1367–1388,2013 www.geosci-model-dev.net/6/1367/2013/ G.Munhoven: Solvingthealkalinity–pHequation: SolveSAPHE 1373 3.1.2 Alternativemethods convergence or no convergence at all). If the procedure is convergent,therateofconvergenceislinear. There are other methods to derive [H+] from AlkC and CT. Thisplainfixed-point-iterationICACmethodwasrecently AllofthemultimatelyseemtorelyontheformulaeofPark madepopularagainbyFollowsetal.(2006).Theseauthors (1969) for deriving the complete speciation of the carbon- argue that in carbon cycle model simulation experiments, atesystemdirectlyfromAlkCandCT,withoutexplicitlyus- where there is little change in pH from one time step to the ing [H+]. Antoine and Morel (1995) first calculate [CO2] next, a single iteration may already provide a sufficiently from CT and AlkC (which involves the solution of a first accurate estimate of [H+] to derive acceptable pCO2 esti- parabolicequation),andthenderive[H+]fromtherelation- mates,foranychosenapproximationoftotalalkalinity.Fol- ship[CO2]=AlkC[H+]2/(K1[H+]+2K1K2),whichrequires lows et al. (2006) suggest, if necessary, to repeat the fixed- thesolutionofasecondparabolicequation.Ridgwell(2001) pointiterationuntilasufficientlyaccurateestimateisfound. firstdeterminesthecompletespeciationofthecarbonatesys- There are a number of models that rely on the ICAC ap- tem,referringfortheadoptedproceduretoMilleroandSohn proach for their pH determination. Peng et al. (1987) con- (1992),whoactuallyonlyreporttheformulaeofPark(1969). sider Alk plus the contributions from silicic and phos- CBW Hethenderivestwodifferentestimatesfor[H+],basedupon phoric acid systems in their representation of total alkalin- thedefinitionsofthefirstandseconddissociationconstants ity.1Theyuseaninitialvalueof10−8andstoptheiriterations ofcarbonicacid,andfinallyusesthegeometricmeanofthese once|(1H)/H|<0.005%.Theyreportthatlessthantenit- twoestimatesasasolutionforEq.(7). erations are generally sufficient. Antoine and Morel (1995) There are no obvious advantages for calling upon these adopt Alk as an approximation to Alk . At each step, CBW T methodsinsteadofthedirectquadraticsolutionabove.Even they derive [H+] from C and Alk by using their special T C ifcarefullyimplemented,bothrequireasignificantlyhigher procedure described above. They iterate until two succes- numberofoperationsthanthesolutionoutlinedabove.Those siveAlk estimatesdifferbylessthan10−8(nounitsgiven). C methodsofferadirectaccesstocarbonatespeciation(atleast Ridgwell (2001) adopts Alk +[OH−]+1.1[PO3−] as an CB 4 inpart),whichcan,however,alsobecalculatedatlittleextra approximationtototalalkalinity.Hecalculates[H+]ateach costfrom[H+]. stepfromC andAlk byusinghisownproceduredescribed T C above.GENIE(Ridgwelletal.,2007)initiallyusedthesame 3.1.3 Iterativecarbonatealkalinitycorrectionmethods procedure as Ridgwell (2001); in more recent versions of GENIE, a complete representation of the phosphoric acid In most common natural settings, the difference between component is used (A. Ridgwell, personal communication, Alk and Alk , albeit small, leads to significant errors on C T 2012). Arndt et al. (2011) use Alk +[HS−] as an ap- [H+],ifAlk isusedinplaceofAlk andoneoftheproce- CBW T C proximationtototalalkalinityinGEOCLIMreloaded.They duresaboveisusedtocalculateitfromC .Toovercomethis T continuetoiterateuntil|Alk +[HS−]−Alk |<10−6(no problem, Alk can be estimated from Alk , and then itera- CBW T C T units given). The method is further used in LOVECLIM tivelycorrecteduntilstabilizationoccurs.Suchaprocedure, (A. Mouchet, personal communication, 2012) with Alk which we call here iterative carbonate alkalinity correction CBW asanapproximationfortotalalkalinity(Goosseetal.,2010) (ICAC) can a priori be used with arbitrary chemical com- andmostprobablystillinsomeothersthat,unfortunately,do positions, provided Alk represents a significant fraction of C notprovidedetailsaboutthecalculationproceduresadopted. Alk . If Alk makes up only a small fraction of Alk , the T C T Bacastow (1981) proposed a variant to improve the rate methodislikelytoexhibitunstablebehaviour. of convergence of fixed-point iterations. That variant only InthemoststraightforwardICACmethod,onestartsfrom usestherecurrencedescribedaboveforthefirsttwoiterates. atrialvalueH for[H+],afirstestimateAlk isobtained 0 C,0 From the third iteration on, Bacastow (1981) switches to bysubtractingtheconcentrationsofallnon-carbonatecom- a secant method to solve the fixed-point equation H− ponents from Alk . That Alk is then used to calculate T C,0 Q(Alk (Alk ,H))=0.2Fixed-pointiterationsarethusonly a new (improved) estimate H for [H+] from Eqs. (9) and C T 1 used to provide starting values for the solution of the fixed- (10) or one of the alternative methods. H is then used to 1 pointequationbythesecantmethod.Therateofconvergence calculateanewestimateAlk fromAlk asaboveandthe C,1 T ofthemethodisstronglyincreasedbythisapproach(andthe procedureisiterateduntilsomepredefinedconvergencecri- terion is fulfilled. This procedure is a classical fixed-point 1Peng et al. (1987) adopt, however, a slightly different defini- iteration: tionoftotalalkalinitybysystematicallyweightingspeciesbytheir respectivecharge.Thisleadstodifferenceswiththephosphoricacid Hn+1=Q(AlkC(AlkT,Hn),CT). (11) system: e.g. the definition of Peng et al. (1987) is equivalent to adoptingm=0forthephosphoricacidsystem. In this recurrence, AlkC(AlkT,Hn) is the estimate of AlkC 2Bacastow(1981)actuallysolvesthealkalinityequationforthe obtained from AlkT by subtracting all the non-carbonate scaledinverseof[H+].Weprovidecodesforthetwoapproaches, components estimated by using Hn. Pure fixed-point itera- althoughweonlybaseourdiscussionsontheversionwithsecant tive schemes may be prone to convergence problems (slow iterationson[H+]. www.geosci-model-dev.net/6/1367/2013/ Geosci. ModelDev.,6,1367–1388,2013 1374 G.Munhoven: Solvingthealkalinity–pHequation: SolveSAPHE domainofconvergenceslightlyenlarged–seenumericalex- 3.2.2 Efficientstartingvalueforiterativemethods perimentsbelow).However,forsomeC -Alk combinations T T theunderlyingfixed-pointequationmaystillgiverisetocon- An excellent initial value for the Newton–Raphson scheme vergence problems, even with the secant method. However canbefoundbyadoptingthefollowingprocedure: as will be shown below, the method of Bacastow (1981) is stronglypreferableoverthepurefixed-pointscheme. 1. locatethelocalminimumclosesttothelargestroot–if TheHadleyCentreOceanCarbonCycle(HadOCC)model itexists,itistheextremum; (PalmerandTotterdell,2001)usesBacastow’smethodforits carbonate speciation calculation, with the AlkCBW approxi- 2. develop PCB([H+]) to second order around that mini- mation. mum;and 3.2 Carbonateandboratealkalinitybasedsolution 3. determinethegreatestrootoftheresultingparabolaand useitasastartingvalue. Only a few models appear to use pH calculation rou- tines based upon Alk . MBM-MEDUSA (Munhoven and That local minimum, if it exists (i.e. if c2−3c >0), is lo- CB 2 1 Franc¸ois,1996;Munhoven,1997,2007)isoneofthem,the catedat modelofMarchaletal.(1998)isanotherone. q 3.2.1 Basicformulationandsolutionmethods H = −c2+ c22−3c1 = −c1 . min q 3 c + c2−3c The equation to solve for [H+] is, for given Alk , C and 2 2 1 CB T B , T TheTaylorexpansiontosecondorderinH ,thusintersects min R ([H+])≡C K1[H+]+2K1K2 theH axisontheright-handsideofHminat CB T[H+]2+K [H+]+K K 1 1 2 v K u P (H ) +B B −Alk H =H +u− CB min , T[H+]+K CB 0 min t q B c2−3c =0. (12) 2 1 This equation may be converted into the polynomial equa- provided P (H )<0. By completing the Taylor expan- CB min tion: sion to third order, it is straightforward to show that H is 0 P (cid:0)[H+](cid:1)≡[H+]3+c [H+]2+c [H+]+c =0, (13) greaterthantheroot. CB 2 1 0 Theso-definedH providesanexcellentstartingvaluenot 0 with onlyforsolvingthecubicpolynomialequation,butalsofor (cid:18) B (cid:19) (cid:18) C (cid:19) otheriterativemethods. c =K 1− T +K 1− T , 2 B 1 Alk Alk CB CB 3.3 Carbonate,borateandwaterself-ionization (cid:18) (cid:18) (cid:19) (cid:18) (cid:19)(cid:19) B C C c =K K 1− T − T +K 1−2 T , alkalinity 1 1 B 2 Alk Alk Alk CB CB CB c =K K K (cid:0)1−2CT+BT(cid:1). WithAlkCBW,CTandBTgiven,theequationtosolveis 0 1 2 B Alk CB R ([H+]) Following our discussion in Sect. 2.2.2, Eq. (12) has a pos- CBW itive root if and only if 0<AlkCB<2CT+BT; if there is ≡C K1[H+]+2K1K2 +B KB a positive root, it is unique. The same holds for the cubic T[H+]2+K [H+]+K K T[H+]+K 1 1 2 B Eq.(13). K + W −[H+]−Alk =0. (14) The cubic equation could possibly be solved with closed [H+] CBW formulae,suchasCardano’sformulae(whichmay,however, suffer from precision problems, require numerically expen- One may either solve this equation in that rational fraction sive cubic root evaluations or possibly complex arithmetic) form with some iterative root-finding method or by one of or Vie`te’s trigonometric formulae (which require a combi- theICACmethodsdescribedabove,oronemaytransformit nation of an arccosine, a cosine and a square root). When intoaquinticpolynomialequation: adopted,thecubicEq.(13)isthereforegenerallysolvednu- mericallywithaNewton–Raphsonscheme.Inthiscase,de- P ([H+]) CBW termininganadequatestartingvalueisthemainproblemto ≡[H+]5+q4[H+]4+q3[H+]3+q2[H+]2+q1[H+]+q0 (15) address in order to design a robust and fast solution algo- rithm. =0 Geosci. ModelDev.,6,1367–1388,2013 www.geosci-model-dev.net/6/1367/2013/ G.Munhoven: Solvingthealkalinity–pHequation: SolveSAPHE 1375 with sibly be used to address this problem. However, they bear somepotentialpitfalls:despitehavingasolution,theunder- q4=AlkCBW+K1+KB, lying fixed-point equation may be difficult to solve numeri- q3=(AlkCBW−CT+KB)K1 cally;intermediateestimatesofAlkC maygooutofbounds +(Alk −B )K +K K −K , (remember that AlkC may only take values between 0 and CBW T B 1 2 W 2C ).ICACmethodscanthereforenotbeguaranteedtofind q =(Alk −2C +K )K K T 2 CBW T B 1 2 thesolution. +(AlkCBW−CT−BT)K1KB−K1KW−KBKW, Theonlycommonlyusedcarbonatechemistryroutinethat q =(Alk −2C −B )K K K directly solves the rational function form of the equation 1 CBW T T 1 2 B −K K K −K K K , is that from the Ocean Carbon Cycle Model Intercompari- 1 2 W 1 B W son Project (OCMIP). For the purpose of that project, Orr q =−K K K K . 0 1 2 B W et al. (2000) prepared standard carbonate speciation rou- tines.TotalalkalinityisapproximatedbyAlk ’[HCO−]+ Thepolynomialequationcanthenbesolvedwithappropriate T 3 standard root finding techniques, selecting the positive root 2×[CO2−]+[B(OH)−]+[OH−]+[HPO2−]+2×[PO3−]+ 3 4 4 4 found.Equations(15)and(14)havethesameuniquepositive [H3SiO−4]−[H+]f−[HSO−4]−[HF]−[H3PO4]. The different root: when Eq. (14) is multiplied by the product of all the species concentrations were, as above, expressed as a func- denominatorsofthefractionsincluded–aproductthatdoes tion of the total concentrations of their respective acid sys- notchangesignfor[H+]>0–totransformitintoEq.(15) temsandof[H+].Theresultingequationwasthensolvedfor nonewsignchangescanbeobtainedfor[H+]>0. pH by a hybrid Newton-bisection method, based upon the Alk isprobablythemostcommonlyusedapproxima- rtsafe solver from Press et al. (1989). All of the models CBW tion for total alkalinity in global carbon cycle models of all thatparticipatedinOCMIPhadtousetheprovidedroutines kinds of complexity. It was already adopted by Bacastow for a set of well defined experiments. A number of models and Keeling (1973), who based their pH calculation on the stillroutinelyusetheseOCMIProutinesfortheirpHcalcu- quinticEq.(15),whichtheysolvebyNewton’smethod,with lations. These include some versions of the Bern3D model astoppingcriterion|(1H)/H|<10−10.Hoffertetal.(1979) (Mu¨ller et al., 2008) and the NCAR global coupled carbon adopt the same procedure (for which they refer to Keeling, cycle–climate model CSM1.4-carbon (Doney et al., 2006). 1973andBacastowandKeeling,1973),butwithalessstrin- As mentioned above, PISCES (Aumont and Bopp, 2006) gent stopping criterion |(1H)/H|<10−6. Keeling (1973) includes a version of the OCMIP solver trimmed down to usesavariant,whereCTisreplacedbyanequivalenttermin AlkCBW only.OthermodelsstilloffertheOCMIPsolversas pCO . anoption. 2 Asalreadymentionedabove,LOVECLIM(Goosseetal., 2010) and HadOCC (Palmer and Totterdell, 2001) use 3.5 Otherapproaches Alk as an approximation for total alkalinity. Alk is CBW CBW Luff et al. (2001) have provided a suite of pH calculation also used in the PISCES model (Aumont and Bopp, 2006), routinesmainlymeanttobeusedinreactivetransportmod- followingasimplifiedversionoftheOCMIPstandardproto- els, but suitable for general speciation calculations as well. col(seenextsection).PISCESisincludedinNEMOandin The methods proposed by Luff et al. (2001) solve the com- some versions of the Bern3D model (Gangstø et al., 2011). pletesystemofequationsthatcontrolthechemicalequilibria Other models that base their pH calculation on Alk CBW between the individual species considered in the total alka- include the Hamburg Model of the Ocean Carbon Cycle linity approximation. These are required for grid-based re- (HaMOCC) family (Maier-Reimer and Hasselmann, 1987; activetransportmodelswheredifferentspeciesarediffusing Heinze et al., 1991; Maier-Reimer, 1993; Maier-Reimer atdifferentdiffusivities.Forcommonapplicationsinbiogeo- et al., 2005), the models of Bolin et al. (1983) and Shaffer chemicalcarboncyclemodels,thisapproachisnevertheless etal.(2008).Nodetailsregardingtheadoptedsolutionalgo- unnecessarilycomplex. rithmsareprovided,though. There are still some other fine pH solvers, such as 3.4 Morecompleteapproximations:rationalfunction CO2SYS of Lewis and Wallace (1998) and derivatives basedsolvers (spreadsheet versions, MATLAB versions, etc. – see http: //cdiac.ornl.gov/oceans/co2rprt.html for more information), When additional components in total alkalinity need to the MATLAB routines from Zeebe and Wolf-Gladrow be considered besides carbonate, borate and water self- (2001), or the R packages seacarb (Lavigne and Gattuso, ionization, converting the resulting rational function equa- 2012) and AquaEnv (Hofmann et al., 2010, 2012). These tion to an equivalent polynomial form becomes more and are, however, generally not suitable for inclusion in global moretediousandtherationalfunctionformbecomesthepre- biogeochemicalmodels,astheyweredevelopedwithspecial ferred basis for finding the solution. ICAC methods are the programmingenvironmentsinmind.Theirfocusismoreon only ones that we have encountered so far that could pos- data processing or modelling with the special programming www.geosci-model-dev.net/6/1367/2013/ Geosci. ModelDev.,6,1367–1388,2013 1376 G.Munhoven: Solvingthealkalinity–pHequation: SolveSAPHE environmenttheyweredesignedfor.Astheirnamesalready difference between the values is about ±3.5%; in log units, suggest, they are mainly aimed at carbonate speciation cal- thevaluesdifferby±0.016. culations. They also often offer the possibility to chose any There is abundant literature on pH scales for seawater. twoamongpH,[CO ](orpCO ),[HCO−],[CO2−],C ,or Besides the original fundamental papers by, e.g. Hansson 2 2 3 2 T Alk tocalculatealltheothers. (1973), Bates and Culberson (1977), Khoo et al. (1977), T DicksonandRiley(1979),Bates(1982)orDickson(1990), theclassicalreviewpapersbyDickson(1984,1993),orstan- 4 pH-scaleconsiderations dard textbooks (e.g. Zeebe and Wolf-Gladrow, 2001), there arealsoseveralrecentpapersonthesubject,suchasthere- Asshortlymentionedabove,thereareafewsubtletiesrelated viewsbyDickson(2010)andMarionetal.(2011)orthere- topHscalesthatstillneedtobeclarified.Themereexistence searchpaperbyWatersandMillero(2013).Inthefollowing ofmorethanonepHscalereflectsthedifficultiestoapplythe sections,wewillthereforeonlygiveacomparativelygeneral fundamentaldefinitionofpH(whichinvolvesanimmeasur- overview, which we have nevertheless tried to keep as self- ablequantity–seenextsection)totheexperimentaldetermi- consistentaspossible. nation of acidity in seawater. All of our calculations never- theless rely on the availability of equilibrium constants that 4.1 FundamentaldefinitionofpHandstandard have to be experimentally derived and we therefore have to potentiometricdeterminationofpH care about differences arising from the usage of various pH scales. While pH as a measure of the acidity of a solution may Let us, similarly to Bates and Culberson (1977) consider appear as a straightforward concept, its experimental deter- theequilibriumrelationship(themassactionlaw)foranacid mination and interpretation are not. The fundamental def- dissociation reaction. Without loss of generality, we may inition of pH recommended by the International Union of write that relationship for the first dissociation reaction of Pure and Applied Chemistry (IUPAC, Buck et al., 2002) ourgenericacidfromSect.2.2.4: states that pH=−log(aH+), where aH+ denotes the activity K1= [H+][[HHnnA−]1A−]. otaiHfo+nth=eofHγHH+++i[oHtnh+sro]i.unTgshhoeltuhateciotianvc.ittaiyvHict+yoeiscffiorceeiflefiantcetideonfttoHγt+Hhe+d,ceospuneccnhednstthroaan-t the exact chemical composition of the solution. The more When the dependency of K on temperature and salinity is 1 experimentallydetermined,thefraction[Hn−1A−]/[HnA]is dilutethesolutionis,thecloserthevaluesofactivitycoeffi- measuredorcalculatedforeachexperiment.[H+]cannotbe cients come to one. The activity of an individual ion in so- lution cannot be measured by any thermodynamically valid directly measured, but gets assigned a value from some pH measurementviathereverserelationship[H+]=10−pH.Tak- methodandthemeasurementofpHthereforerequiresanop- erationalconvention(Bucketal.,2002).Thereasonsforthe ing the negative logarithm (antilogarithm) of the previous equationandwritingpK =−logK ,weget existence of several pH scales in seawater then also simply 1 1 “[...] reflect the gradual gradual refinement of the experi- pK =pH−log(cid:18)[Hn−1A−](cid:19). mentallyconvenientpotentiometricdeterminationofacidity 1 [H A] in order that the numbers obtained might be usefully inter- n preted as a property of hydrogen ion in solution” (Dickson, In a given setting (i.e. for given temperature, salinity, pres- 1984). sure,solutionchemistry,etc.),theratio[Hn−1A−]/[HnA]is ThepotentiometricmethodmentionedbyDickson(1984) set and different calibrations of the pH-meter used, i.e. dif- is the classical method used for the quantitative determina- ferentscaleschosenforthepH-meter,willthusproducedif- tion of acidity in an aqueous solution. It is based upon the ferent pK values. Any experimentally derived parameteri- useofelectrochemicalcellsandhasbeenusedformorethan 1 zationforK canthereforeonlybeusedinconjunctionwith 100yr. Potentiometric pH measuring devices for practical 1 aH+ concentrationscalethatisconsistentwiththepHscale use are made up by two electrodes: an H+ sensitive glass- that was used to derive it. Before a particular empirical pa- electrode and a well reproducible second electrode, a so- rameterizationforK canbeusedwithadifferentscaleofpH calledreferenceelectrode.Bothelectrodesareimmersedinto 1 (e.g.duetoadifferentconventionaloroperationaldefinition thesamplesolutiontoformanelectrochemicalcell.Thepo- ofpH),itmustbeconverted. tential difference between the two electrodes, i.e. the emf Additional conversion may be required because of (electromotive force) of the cell, is linked to the logarithm the usage of different concentration units: both mo- of the activity of the H+ ions in solution. The total emf of lal units (mol/kg-H O) and mol/kg-seawater are com- the cell, E, can be separated into three major contributions 2 mon. They can be converted according to [mol/kg-SW]= (Dickson, 1984). The first one, which we denote here as m(mol/kg-H O)×(1−0.001005S) (Dickson et al., 2007, E ,isduetothepotentialdifferenceacrossthemembrane 2 gem chapter 5, p. 13), where S denotes salinity. For S=35, the oftheglasselectrode,whichisassumedtobehavefollowing Geosci. ModelDev.,6,1367–1388,2013 www.geosci-model-dev.net/6/1367/2013/