Interactive comment on “ Comparison of seven packages that compute ocean carbonate chemistry ”

One or two pHScale conversion factor(s) are passed as optional parameter(s) to speed-up computation: (i) kSWS2scale to convert from the seawater scale (SWS) to the pH scale selected at the hydrostatic pressure value indicated and (ii) ktotal2SWS_P0 to convert from the total scale to the SWS at an hydrostatic pressure of 0.These conversion factors are calculated using function kconv() if they are not given. Computations are vectorized rather than using loops. Warning messages for out-of-validity-domain: only one message per constant (instead of one per data entry).

However, the practical user may be more interested in evaluating the package comparison in light of realistic measurement uncertainties.For example, Dickson (2010) provides the following estimated uncertainties for a single measurement on a sample of surface seawater: TA: 2-3 umol/kg DIC: 2-3 umol/kg pH: 0.005 pCO2: 2 uatm This applies to state-of-the-art methods using reference materials and: "... performed by an experienced laboratory with well-trained analysts, and with a good quality assurance program in place."In fact, all differences between relevant surface variables from the various packages shown in Orr 's Figs. 2,3,5,6*,7,8,9,10,and 11 are smaller than the measurement uncertainties cited above.This even includes "pCO2" from csys (which is actually fCO2 and agrees with CO2SYS' fCO2 to within ∼0.1 uatm, see below).*One exception is Fig. 6 for salinities < 10, where K1 and K2 from Lueker et al. (2000) cannot be applied.Note that Dickson's (2010) Table 1.5 referred to in Orr's manuscript gives uncertainties for *Reference Materials* distributed by Dickson's laboratory.These are *not* uncertainties for typical measurements performed in the user's laboratory.
We agree.In regard to Dickson's Table 1.5, our Discussion paper referred to the "best measurement uncertainty", not the uncertainty for typical measurement.For numerical packages not to add substantially to total uncertainty, they should agree to within much less than the measurement uncertainty (as detailed in our section 2.3).Since we arbi-C4421 trarily divided the measurement uncertainty by ten to obtain a threshold for numerical uncertainty, changing the best measurement uncertainty to the typical measurement uncertainty values given by Drs.Zeebe and Wolf-Gladrow would not substantially alter our evaluation of when differences between packages should be considered significant.Furthermore, the errors listed above are typically considered as random errors.Those are fundamentally different than systematic errors (biases) such as those due to errors in calculation programs.
In summary, the presented package comparison may be put into perspective as follows.For the vast majority of users dealing with measurement uncertainties in surface samples equal to or larger than those given above, there is virtually no difference as to which carbonate chemistry package they prefer to use.The reason is that in most cases the measurement uncertainties vastly exceed differences between numerical routines.
The comment above is a fair statement for the findings of our comparison when packages used the set of constants recommended for best practices (Dickson et al., 2007).However, with K 1 and K 2 from Millero (2010), packages disagreed by much more than the numbers given above.Furthermore, we are not convinced that all users would be happy to use a package that was known to give a biased answer even if that bias was less than the random errors listed above.
In any case, the purpose of this comparison was to compare packages quantitatively and to report the results.Because of our Discussion paper, the community now knows that the packages, when used with the best-practice constants, do agree closely enough to be used interchangeably for most studies.Yet even small differences among packages (e.g., when all of them use K 1 and K 2 from Lueker et al., 2000) seem to be of interest to some marine chemists as well as to package developers, three of which have responded with updated versions.

C4422
The reader should note that the "pCO2" from csys as plotted in Orr's manuscript (Figs 2,3,5,6,8) is actually fCO2 (fugacity), rather than pCO2.At the time when we (Wolf-Gladrow and Zeebe) started csys (prior to 1993) and later provided csys as a supplement to our book (Zeebe and Wolf-Gladrow, 2001), we never imagined that csys would be used for the purpose of sub-uatm calculations/inter-comparisons.In that case (which applies to Orr's comparison), one needs to take into account the difference between fCO2 and pCO2, which is about 1 uatm at typical surface ocean conditions (Zeebe and Wolf-Gladrow, 2001, Chapter 1.4).In fact, csys' surface fCO2 agrees with CO2SYS within 0.1 uatm or so.
We state clearly in the Discussion paper that the pCO 2 variable in csys was actually f CO 2 , a statement for which we are later criticized for by Dr. Zeebe.That finding was never published previously as far as we are aware.We presume that the mislabeling of f CO 2 as pCO 2 in csys was unintentional.The documentation of csys did not contain any warning about when not to use the mislabeled pCO 2 variable.The authors were well aware of the fact that csys' "pCO2" is actually fCO2 (see their Abstract and Conclusion section), so it remains unclear why they compare this variable to pCO2 from other programs on a sub-uatm scale.If a revised version of the ms will be invited, this should be corrected.
For the comparison, we simply compared all variables that were labeled pCO 2 in all packages.Because that variable from csys differed slightly, we thought it might actually be f CO 2 , and we confirmed that by studying the source code and comparing it with f CO 2 from the other packages.So yes of course, we were aware of it, as clearly stated in the Discussion paper.Should we have hidden it?
By being open about this difference, the community now knows that calculated pCO 2 made with csys and published previously contains the same error.We are happy that the csys developers took this finding into account and have now provided a new version of csys that rectifies this confusion between pCO 2 and f CO 2 , while providing not one C4423 but both variables as output.In the revised manuscript, we will not hide this finding, but it will be be much less apparent.It will be confined to only 1 figure that discusses how old versions of packages diverge from the reference.Elsewhere, we will show the pCO 2 from the most recent version of csys, which has been corrected.
For users dealing with high-quality data/applications, csys has been updated to compute and output both fCO2 and pCO2.We have also added a user option to csys if one wishes to compare csys output to that of CO2SYS: http://www.soest.hawaii.edu/oceanography/faculty/zeebe_files/CO2_System_in_Seawater/csys.htmlThank you for this correction and the useful addition of the option to facilitate comparison with CO2SYS.
However, note that in order to resolve measurement differences of order 1 uatm, high laboratory standards and reference materials are required (see above and below).
We think there is no reason not to correct a systematic bias if it is known, even if it is relatively small.Our Discussion paper made it known, and the cys developers corrected it.
(3) The more daunting issues A realistic error assessment of CO2 system calculations needs to consider the uncertainties associated with input variables and the fundamental constants used in the calculations.For example, pCO2 may be calculated from measured TA and DIC (a "good" combination as opposed to e.g.calculating DIC from pH and pCO2).To estimate the error in calculated pCO2, one needs to consider (at least) errors in TA, DIC, pK1, and pK2.Assuming realistic uncertainties of +-3 umol/kg in TA and DIC and +-0.006/+-0.011 in pK1 and pK2 (e.g.Millero et al., 2006), the uncertainty in calculated pCO2 (DpCO2) may be estimated as: where the four terms of the sum are the individual squared uncertainties owing to uncertainties in TA, DIC, pK1, and pK2 at T = 25C, S = 35, P = 0, TA = 2400 umol/kg and DIC = 2080 umol/kg.This uncertainty is about 100 times larger than the difference between numerical packages (order 0.1 uatm) as discussed by Orr.
Does this propagation of uncertainty, which is intended for random errors, mean that as long as systematic errors in calculated pCO 2 are less than 11 µatm, then packages can be used interchangeably?Does it mean that we should have never mentioned the systematic difference in pCO 2 of more than 1 µatm from csys in the Discussion paper?Other fundamental issues for CO2 system calculations include an apparent large effect on pK2 at rising pCO2/DIC (Millero et al., 2002) and possible inconsistencies between parameters when over-determining the system (Hoppe et al., 2010).Regrettably, no suggestions for solutions of these more daunting issues are offered by Orr.Neither are the implications of the numerical package comparison discussed on the background of the fundamental uncertainties listed here.
The focus of our Discussion paper is on discerning systematic differences between packages that compute carbonate chemistry.But that comparison did not put us in a position to resolve all other open issues.Actually, we did mention the work by Hoppe et al. (2012) in our Discussion paper (see p. 5331, lines 8-12).We pointed out that we could not contribute to resolve their dilemma of large inconsistencies between calculated and measured pCO 2 , differences that are much worse than found by studies from marine chemists.
(4) Pressure corrections The significance of the discussion and the recommendations made regarding pressure corrections and scale conversions as described at length in the manuscript is difficult to comprehend (see Secs. 2.7,4.2.1,4.2.3,4.2.4,4.2.6,and 5).The key issue here is to recall that measurements of pressure effects on acid-base equilibria in seawater are sparse and/or uncertain and that estimated values from molal C4425 volume and compressibility have substantial uncertainties.For example, the estimated P-correction for KS involves steps using HCO3-(!) as a model for HSO4-; the measured and calculated value for KS in water seem to differ by some 10% at 1000 bar (Millero, 1983).For K1 and K2, the P-corrections from Millero (1983) are based on Culberson and Pytkowicz's (1968) data (short CP68), which were obtained in artificial seawater.CP68 found a difference in the K2 P-correction of 6% relative to previous work and also found significant differences relative to a few measurements made in natural seawater.Millero's (1983) calculated P-corrections for K1 and K2 differ by up to 3% and 8% from CP68's data at 1000 bar and 2 degC.Overall, the uncertainties for the P-corrections for K1 and K2 could also be well of order 10%.
These are all good points.In the Discussion paper, our focus was on comparing results from the seven packages that all used the same approach from Millero (1995) to make pressure adjustments.None of the packages offers another approach.In the revised manuscript, although our focus will remain the same, we will also discuss more general uncertainties in making pressure adjustments to the equilibrium constants.
In the manuscript, Orr suggests that constants should be first converted to the seawater scale (SWS), then P-corrections be applied on the SWS, and finally constants be convert back to say, the total scale.Now the scale conversion depends on KS and KF, which are themselves P-dependent, including large uncertainties (e.g.probably more than 10% for KS).Note that as a result of pressure change, the ratio of the scale conversion at P = 0 and 1000 bar may differ by perhaps 0.4%.
Our discussion paper already states that making pressure corrections without converting to the seawater scale, as done by csys, caused only small differences, e.g., 0.5% for K 1 and K 2 at 4000 db.We wish to clarify though that it is not us who first suggested that the pressure corrections should be made on the seawater scale.Rather it was Lewis and Wallace (1998), who based that assertion on an extensive study of the literature associated with the pressure adjustment proposed by Millero (1995).The same approach is used in all packages except csys.This does not mean that csys C4426 is wrong, but it is inconsistent with the others.Once again our focus is on precision among packages.
In summary, this would mean one should apply a correction of ∼0.4% to a value that is uncertain by perhaps 10%.To make it worse, the 0.4% correction itself is uncertain because the P-correction for the scale conversion is uncertain, owing to its pressuredependence through KS (order 10%?) and KF.The bottom line is that attempts to gain apparent accuracy by applying scale conversions to pressure corrections are completely lost in a sea of uncertainty that traces back to the original data and theoretical estimates of the pressure coefficients.
It is a small correction and uncertainties are larger.However, in terms of our goal to quantify precision among packages, it is the main reason why at depth, calculated variables in csys differ from those in other packages.We thank the csys developers for recently adding the option to allow users to choose to do the pressure correction as in the other packages or to stay with the original csys approach.
(5) Equilibrium constant for water (KW) Dickson (2007) converts Kw approximately from the SWS to the total scale by subtracting 0.015 from the constant term, which appears well-justified given the much larger fundamental uncertainties described above.We have followed this approach.Some may believe one could strive for more precision than the very laboratory that currently supplies the reference materials for the quality control of ocean CO2 measurements.I do not.
The simplified approach used in cys is certainly legitimate; it follows the approach recommended by Dickson et al. (2007).In our Discussion paper, we only say that the other packages use a variable offset approach, which is more sophisticated, but we do not say it is better.We thank the csys developers for recently introducing an option to allow the user to choose either approach.

C4427
The pressure coefficients a0, a1, and a2 for KW given in Millero (1995, Table 9) appear to be for water rather than seawater (Millero, 1983).They have been changed in csys.
We appreciate that the csys developers have now fixed these 3 pressure-adjustment coefficients for K W , errors which we identified in our Discussion paper.
(6) Equilibrium constant for hydrogen fluoride (KF) Orr devotes the better part of an entire section (4.2.4) discussing and analyzing KF in our code.However, KF was never used to calculate output because csys works on the total scale (the code lines for KF had been kept from a different version).Thus, KF has zero effect on csys output on the total scale.If a revised version of the ms will be invited, the discussion of KF in csys and the comparison of KF with other packages (Fig. 17) should be deleted.Subsection 4.2.4 in the Discussion paper contained only 2 paragraphs on K F .The first of those points out a common error in both seacarb and csys, an erroneous formulation for ionic strength.The second paragraph further details discrepancies in K F from csys.Since the ionic strength error was corrected in both packages and other discrepancies have also been remedied, our revised manuscript will be naturally much more concise in regards to K F .
The second sentence in the above comment was true for the version of csys that was compared in the Discussion paper.But K F is now used in the latest version of csys, released afterwards.With csys's new option to improve comparison with CO2SYS, K F is used to convert from the total scale to the seawater scale before pressure corrections of some equilibrium constants and then back again afterwards.We will modify discussion in the revised manuscript in order to account for these changes.
(7) Equilibrium constants for phosphoric acid Also, constants for phosphoric acid (discussed in ms Sec.4.2.6) are never used in csys to calculate any output (the code lines had been kept from a different version).

C4428
Thus, constants for phosphoric acid also have zero effect on csys output.If a revised version of the ms will be invited, the discussion of phosphoric acid constants in csys and the comparison with other packages (Fig. 18) needs to be deleted.
We already pointed out in the Discussion paper that the phosphoric acid constants that were computed by csys were not used.Since csys does compute K 1P , K 2P , and K 3P , and its source code is available and can be modified, we felt it appropriate to notify users of these differences.The latest version of csys still computes these 3 constants, but discrepancies have been reduced.Hence the revised manuscript will say less about the subject.
(8) Comments on csys description Section 3.2 and caption fig 7: "... csys, which does not allow pCO2 as an input variable" is inaccurate.csys allows pCO2 as input variable but only in combination with pH.
In the revised manuscript, we will correct the corresponding text in section 3.2 and in Fig. 7 to be consistent with Table 4, which does show that csys is able to use the pCO 2 -pH pair.Thanks for pointing out this inconsist description.Section 4.2.1:"... csys exhibits problems that can be traced back to its implementation of the Lueker et al. (2000) formulations for K1 and K2." is inaccurate.Differences may be due to pressure corrections but not due to the implementation of K1 and K2 based on Lueker et al. (2000).K1 and K2 agree within round-off error with CO2SYS at P = 0.
We thank the csys developers for pointing out this imprecision.We do say earlier in the same section that for K 1 and K 2 from Lueker et al. (2000), all packages agree at the surface.In the revised manuscript, we will modify the sentence mentioned above to be consistent with what was said previously.
Total boron: this variable is part of the input section in csys and can be readily changed by the user if Lee et al.'s (2010) value is preferred.
Certainly users are free to modify whatever they like given the source code.If the de-C4429 velopers offer this new option for total boron in a later release, we would then consider that csys offers it.
We have added a user option to csys if one wishes to compare csys output to that of CO2SYS: http://www.soest.hawaii.edu/oceanography/faculty/zeebe_files/CO2_System_in_Seawater/csys.htmlThis is particularly useful option.Much appreciated.