Comparison of ten packages that compute ocean carbonate chemistry

. Marine scientists often use two measured or modeled carbonate system variables to compute others. These carbonate chemistry calculations, based on well-known thermodynamic equilibria, are now available in a dozen public packages. Ten of those were compared using common input data and the set of equilibrium constants recommended for best practices. Current versions of all 10 packages agree within 0.2 µatm for p CO 2 , 0.0002 units for pH, and 0.1 µmol kg − 1 for CO 2 − 3 in terms of surface zonal-mean values. That represents more than a 10-fold improvement relative to outdated versions of the same packages. Differences between packages grow with depth for some computed variables but remain small. Discrepancies derive largely from differences in equilibrium constants. Analysis of the sensitivity of each computed


Introduction
Our ability to assess ocean carbon uptake and associated impacts from ocean acidification relies on an accurate representation of the marine carbonate system.Fortunately, the seawater carbonate system is well constrained, allowing any two of its variables to be used to calculate all others, given associated temperature T , salinity S, pressure P , and nutrient concentrations.For example, it is common to measure or simulate two conservative variables, dissolved inorganic carbon C T and total alkalinity A T , and then compute from associated thermodynamics the corresponding pH, partial pressure of carbon dioxide pCO 2 , concentrations of aqueous CO * 2 as well as carbonate CO 2− 3 and bicarbonate HCO − 3 ions, and the related Revelle factor and saturation states of aragonite A and calcite C .It is the CO 2 -driven changes in these variables that drive the biological impacts from ocean acidification (Gattuso and Hansson, 2011;Kroeker et al., 2013;Wittmann and Pörtner, 2013) and degrade most the ocean's capacity to absorb CO 2 (Sarmiento et al., 1995;Orr, 2011).
These equilibrium computations are made with numerous software packages, either those developed and used by individual scientists or, more commonly, those that have been made available publicly.The latter have become indispensable for many ocean scientists, whether they study marine chemistry or impacts of ocean acidification on marine biota.Yet how packages differ is seldom addressed.Lewis and Wallace (1998) documented differences in basic variables among three existing packages at a time when no such package was publicly available.Provided with the same input, computed output from the three packages differed by 21 µatm for pCO 2 , 0.16 units for pH, and 15 µmol kg −1 for CO 2− 3 as well as HCO − 3 .Packages used different pH scales, different formulations for some of the constants (K 1 , K 2 , K B , and K S ), and different definitions of total alkalinity, all apparently hard-coded.These differences prompted Lewis and Wallace (1998) to develop a publicly available package, CO2SYS, which provides many options to select from the available pH scales and constants while being based primarily on recommendations from Dickson and Goyet (1994).Since that time, other packages have also been developed and released publicly, yet to this day no study has been published that compares their results.One may assume, given continued efforts to establish and refine procedures for best practices (Dickson and Goyet, 1994;Dickson et al., 2007;Dickson, 2010), that differences among currently available packages are less than what was found 17 years ago.But even that poor level of agreement has not been established.A quantitative understanding of the accuracy and precision of these packages is needed to rigorously compare studies that aim to assess, e.g., air-sea CO 2 fluxes and thresholds associated with ocean acidification.
Ten publicly available software packages were included in this comparison.The first was CO2SYS, but that now exists in four different variants: the original program written in QBasic and running on DOS (Lewis and Wallace, 1998), two variants as Excel spreadsheets (Pierrot et al., 2006;Pelletier et al., 2007), and most recently a variant as MATLAB scripts (van Heuven et al., 2011).We will refer to these packages as CO2SYS-QBasic, CO2SYS-Excel-Pierrot, CO2SYS-Excel-Pelletier, and CO2SYS-MATLAB.Another package, csys, was also written in MATLAB, but it was released a decade earlier as a supplement to the book by Zeebe and Wolf-Gladrow (2001).The development of csys inspired seacarb, an R library (R Development Core Team, 2012) released 2 years later (Proye and Gattuso, 2003) and frequently improved with new revisions (Gattuso et al., 2015).About the same time as the release of the two Excel variants of CO2SYS, the swco2 package was also released with a similar spreadsheet interface but a distinct library of core routines written in Visual Basic (Hunter, 2007).Three years later, oceanographers saw the release of two new carbonate chemistry packages, CO2calc and ODV, both of which also exploit the core CO2SYS code.While CO2calc provided a new tool for Mac, PC, and iOS (Robbins et al., 2010), ODV provided carbonate chemistry calculations as an add-on to an already widely used visualization and analysis tool (Schlitzer, 2002).Parallel to those developments for the observational and experimental communities, the Ocean Carbon-Cycle Model Intercomparison Project (OCMIP) provided routines to compute surface pCO 2 and air-sea CO 2 fluxes from simulated A T and C T (Orr et al., 1999).Those were adapted to include the full suite of other carbonate system variables (Orr et al., 2005), then later improved and released publicly as the mocsy package (Orr and Epitalon, 2015).To assess the consistency of these packages, we compared results generated by running them with common sets of constants, pH scales, and input data.
Two other public packages, AquaEnv (Hofmann et al., 2010) and SolveSAPHE (Munhoven, 2013), were not included in this comparison.AquaEnv is particularly suited to aquatic chemical model generation in freshwater and estuaries; however, we found it to be designed for high-end users, e.g., finding no examples to quickly convert all of its results from its default free hydrogen ion scale to the total hydrogen ion scale, as recommend for best practices.SolveSAPHE defines the state of the art for the algorithm used to solve the pH-alkalinity equation because of its greater efficiency and stability.It always converges even under extreme conditions.Its solver routines have already been adopted by one of the packages (mocsy 2.0), but SolveSAPHE itself does not provide an adequate user interface for simple use, as needed for this comparison.
We limit this study to package comparison.For brevity, we avoid redocumenting the associated approaches and algorithms, which are now commonly used and for which abundant literature already exists (e.g., Dickson et al., 2007;Munhoven, 2013).Likewise, we do not address the debate raised by Hoppe et al. (2012) concerning poor agreement between measured pCO 2 and that computed from A T and C T , a disaccord found to be worse than in previous studies by marine chemists (e.g., Lueker et al., 2000).Nonetheless, we go beyond simply identifying differences between packages; we also seek to identify their causes.Our goal was to inspire subsequent package developments while facilitating evaluation and tightening agreement.Once packages are shown to provide essentially the same results, they can be legitimately chosen based on convenience, efficiency, functionality, and a user's programming experience.Some users may well prefer spreadsheet-based programs (CO2SYS-Excel-Pierrot, CO2SYS-Excel-Pelletier, CO2calc, and swco2).Others who use ODV for general oceanographic data analysis and visualization can easily compute carbonate system variables using its predefined derived-variable facility.Users with some programming experience may prefer packages that are available in languages that they are already familiar with.They can choose from CO2SYS-MATLAB and csys in MATLAB, seacarb in R, swco2 in Visual Basic, and mocsy in Fortran.Python programmers can use either mocsy or seacarb.

Packages and reference
All publicly available software packages (Table 1) take two ocean carbonate system variables as input and compute the others from basic thermodynamics.All 10 packages were first downloaded in November 2012.Our earliest findings led developers to update two packages (CO2calc and seacarb).
Another three packages (csys, mocsy, and CO2SYS-Excel-Pelletier) were updated following publication of our discussion paper (Orr et al., 2014).Different results from some older versions of these packages are briefly shown in one figure to illustrate the discrepancies associated with running software that is out of date.The remaining comparison refers only to the latest version of each package.
To compare packages, it was necessary to define a common reference.Although check values exist for most of the equilibrium constants (Dickson et al., 2007), none are available for computed variables.Hence we chose CO2SYS as a relative reference for three reasons: (1) it was the first publicly available package; (2) its core routines already serve as the base code for two other packages (CO2calc and ODV); and (3) its documentation and code reveal the intense effort that its developers have put into ferreting out the right coefficients from the literature and the most appropriate version of formulations for the constants.
Our choice of the reference had to be refined, however, because, as mentioned, CO2SYS comes in four variants: the original in QBasic running on DOS (Lewis and Wallace, 1998), two others that run with Excel (Pelletier et al., 2007;Pierrot et al., 2006), and finally MATLAB code (van Heuven et al., 2011).The original variant is still used by some, but it does not provide options to use formulations for K 1 and K 2 from Lueker et al. (2000), as recommended for best practices (Dickson et al., 2007).Thus, we reduced our choices for the reference to the other three variants of CO2SYS (MATLAB and both Excel variants), all of which provide options to use all constants recommended for best practices (with one minor exception).All three versions give nearly identical results (Fig. 1).But they differ significantly from the original version run with the closest substitutes for Lueker et al.'s K 1 and K 2 , namely the previous refits by Dickson and Millero (1987) of the same measured constants from Mehrbach et al. (1973) (DM87).With that older set of K 1 and K 2 , however, the original version produces results that match those from the Excel-Pierrot version when also run with DM87.Another requirement was efficiency, since one aspect of our comparison required use of a global-scale input data set with nearly 1 million records (Table 2).Thus we further narrowed our choice of the reference to the most efficient variants, CO2SYS-MATLAB and CO2SYS-Excel-Pelletier.Finally, we chose CO2SYS-MATLAB as the reference because it was the most efficient and because its source code was easiest for us to inspect, modify, and rerun for sensitivity tests.This arbitrary choice of the relative reference was necessary for our comparison, but it does not imply that the chosen reference is necessarily error free.

Features
The 10 packages compute results from the same thermodynamic equilibria, but software features differ, including available input pairs, pH scales, and constants.Package diversity covers all commonly used operating systems: all packages run on Windows, eight packages on Mac OSX, and five packages on Linux and Unix (Table 3).Source code is available in seven packages in standard programming languages (QBasic, Visual Basic, MATLAB, R, Fortran 95), thereby allowing code validation and improvements by users on all three operating systems mentioned above.The number of possible input pairs of carbonate system variables varies widely between packages, from 1 to 20 (Table 4).The mocsy package  The three most recent variants (MATLAB and both Excel versions) are run with constants recommended for best practices (BP).The QBasic variant does not offer the same K 1 and K 2 , so we used an earlier refit by Dickson and Millero (1987) of the same data (DM87) and compared it to one Excel version also with DM87.
treats only one input pair: A T -C T , the two carbonate system variables carried by models.The four CO2SYS variants and its two derivatives (CO2calc and ODV) allow the user to select from six commonly measured pairs (A T -C T , A T -pH, A T -pCO 2 , C T -pCO 2 , C T -pH, and pH-pCO 2 ); in addition, they allow equivalent pairs where f CO 2 replaces pCO 2 .The csys package provides 10 more input pairs by allowing pair members to include one or more of the three inorganic carbon species CO * 2 , HCO − 3 , and CO 2− 3 .Although the two former species can only be calculated, promising new techniques are being developed to measure the latter (Byrne and Yao, 2008;Martz et al., 2009;Easley et al., 2013).Yet despite csys's enhanced number of input pairs, it limits pCO 2 to be used as input only when combined with pH.The two remaining packages, seacarb and swco2, include the same 16 pairs as csys but also add four others, all including pCO 2 .
Computed variables are affected by the choice of the pH scale and the constants.All packages allow users to work on the total pH scale as recommended for best practices (Table 5) and as used for this comparison.The mocsy package provides only the total scale, while the others allow for conversion to the free scale.The others also allow users to work on the seawater scale except for csys.The 4 CO2SYS vari-ants as well as CO2calc and swco2 also offer the NBS scale.The choice of the pH scale affects the values of the constants for which H + is part of the equilibrium equation.For K 1 and K 2 , the CO2SYS variants and derivatives offer a large range of choices (Table 6).Yet most of those may now be considered out of date, having been replaced by more recent assessments, sometimes with some of the same data.All packages except CO2SYS-QBasic offer the K 1 and K 2 formulations from Lueker et al. (2000), as recommended for best practices.Six packages also offer the most recent formulations for K 1 and K 2 that have been proposed as more appropriate for low-salinity waters (Millero, 2010).The formulations for K 1 and K 2 from the two latter studies are used individually in this comparison to assess associated differences between packages.For the other constants, all packages provide the formulations recommended for best practices, except for K F , a difference shown later to have no consequence.Some packages also offer additional features.For example, CO2SYS variants and CO2calc allow users to compute variables at a temperature that differs from the in situ value.Some also distinguish different components of total alkalinity, including those from total B, P, and Si.
a Both variants: CO2SYS-Excel-Pierrot and CO2SYS-Excel-Pelletier. b CO2SYS, CO2calc, and ODV also allow input pairs containing f CO 2 instead of pCO 2 .c seacarb and swco2 include user-callable functions to convert between pCO 2 and f CO 2 .
sion of pH and constants between the free, total, and seawater scales; other packages make such conversions internally but do not provide user-callable functions.The seacarb package also offers functions to help design perturbation experiments to investigate effects of ocean acidification (Gattuso and Lavigne, 2009).Two packages, mocsy and seacarb, al-low users to account for pressure effects on subsurface f CO 2 and pCO 2 following Weiss (1974); other packages neglect these pressure effects.

Input data
To compare packages, we used two different kinds of input data.A first analysis compared variables computed in each package as a function of latitude and depth using as input the three-dimensional gridded data products for A T and C T from the Global Ocean Data Analysis Project (GLO-DAP; Key et al., 2004) combined with comparable products from the 2009 World Ocean Atlas (WOA2009) for T (Locarnini et al., 2010), S (Antonov et al., 2010), and concentrations of total dissolved inorganic phosphorus P T and total dissolved inorganic silicon Si T (Garcia et al., 2010).We will refer to this combined gridded input data as GLODAP-WOA2009.A second analysis focused on comparing pack-

Best-practices comparison
Comparisons were made using the total pH scale and constants recommended for best practices by Dickson et al. (2007).The equilibrium constant for the solubility of CO 2 in seawater K 0 is from Weiss (1974).The equilibrium constants K 1 and K 2 are from Lueker et al. (2000), who refit the constants determined by Mehrbach et al. (1973) to the total pH scale.The formulation for K B is from Dickson (1990b) and is also on the total pH scale.Formulations for K W , K 1P , K 2P , K 3P , and K Si are from Millero (1995), who provides equations for the seawater scale, and those are converted to the total scale.The formulation for K S is from Dickson (1990b) on the free scale (see above).The solubility products for aragonite K A and for calcite K C are from Mucci (1983).All these are equilibrium constants given in terms of concentrations, not activities.The only constant for which the formulation was not that recommended by Dickson et al. (2007) is K F , because that best-practices formulation (Perez and Fraga, 1987) is not offered by most CO2SYS variants, CO2calc, nor ODV.Instead, we used the K F formulation by Dickson and Riley (1979) on the total scale, which is offered by all packages and recommended by Dickson and Goyet (1994).Dickson et al. (2007) state that results from the two formulations are in reasonable agreement.Additionally, all packages used consistent formulations for total concentrations of boron (Uppström, 1974), sulfur (Morris andRiley, 1966), fluoride (Riley, 1965), and Ca 2+ (Riley and Tongudai, 1967), each proportional to salinity.Nine packages compute saturation states for aragonite A and calcite C from the product of the concentrations of Ca 2+ and CO 2− 3 divided by the corresponding solubility product, either for aragonite K A or calcite K C (Mucci, 1983), respectively.Only the csys package does not provide output for A and C .To simplify comparison, most figures plot results as absolute differences: computed values are shown after subtracting off corresponding results from the reference.

Sensitivity tests
The most extensive comparison was made with A T -C T as input, the only pair that is available in all packages.With that pair, packages were also compared in terms of how their computed variables were affected by nutrient concentrations, i.e., by varying P T and Si T across their observed ranges in the ocean.Additionally, the same pair was used to quantify effects of two important developments since the best practices were published in 2007.For the first, we quantified effects on computed variables of Lee's (2010) assessment that the total boron concentration in the ocean may be 4 % larger than considered previously (Uppström, 1974).For the second, we assessed impacts of using Millero's (2010) new K 1 and K 2 formulations, which are designed to cover a wider range of input S and T relative to the intended range for recommended constants (Lueker et al., 2000).

Constants
To better assess the most likely causes of differences in computed carbonate system variables, we also compared associated constants.For the four packages where source code was available and easily modified (CO2SYS-MATLAB, csys, seacarb, mocsy), we used existing routines or slightly modified versions to output all the constants for the same physical input data (T , S, and P ) that we used for computing variables.For CO2SYS-Excel-Pelletier, the constants are also available.For swco2, we retrieved its constants using its documented parameter numbers and its routine to extract anything with a parameter number.For packages where source code was not available, we computed constants from output variables when possible.With output from CO2calc and CO2SYS-Excel-Pierrot, we computed its K 0 , K 1 , K 2 , K B , K W , K A , and K C ; from ODV output, we computed its K 0 , K 1 , K 2 , K A , and K C .

Pressure corrections
Until recently, no public package accounted for pressure effects on K 0 (needed to convert CO * 2 to f CO 2 ) and the corresponding fugacity coefficient C f (needed to convert f CO 2 to pCO 2 ) as originally proposed (Weiss, 1974, Eqs. 5 and 9).Instead, the total pressure term in those equations was simply assigned to be that of the atmosphere only (1 atm).Hence, their computed subsurface f CO 2 and pCO 2 may be considered as being referenced to the surface.These pressure effects are accounted for, however, in the latest versions of two packages, mocsy 2.0 and seacarb 3.0.6.Both allow users to compute f CO 2 and pCO 2 in three ways: (1) the same "com-mon" approach that computes K 0 and C f with total pressure of 1 atm and in situ T , (2) the "potential" approach that likewise uses atmospheric pressure only but also uses potential temperature θ instead of in situ T , and (3) the "in situ" approach that uses the true total pressure (atmospheric + hydrostatic) and in situ T .Other packages offer only the first approach.
For the other equilibrium constants, all packages make pressure corrections following the approach of Millero (1995).That is, the effect of P on each equilibrium constant K i is given by the equation ln where the left-hand side contains the ratio between K i at depth (P in bars) and at the surface (P at 0 bars), R is the gas constant, T k is temperature in K, V i is the partial molal volume, and κ i is the change in compressibility.The latter two variables differ for each constant and were fitted empirically by Millero (1995) to be quadratic in temperature: where T c is temperature in • C. Some of these original coefficients (Millero, 1995, Table 9) contained typographical errors as identified in the code and documentation of CO2SYS-QBasic (Lewis and Wallace, 1998, Appendix).Nonetheless, these errors have persisted in some of the packages as well as in the literature (e.g., Millero, 2007).To help amend this situation, Table 7 lists these coefficients for each constant where known errors have been corrected.To determine the fidelity of packages to this array of coefficients, we studied available source code and evaluated patterns of discrepancies in results by carrying out sensitivity tests to decipher fingerprints characteristic of previous errors.
Although the same approach is used by all packages to make pressure adjustments (Eqs. 2 and 3), it is based on extremely limited data.Thus it may not be particularly accurate.For example for K 1 and K 2 , there are differences of 3 and 8% between adjusted values from Millero (1983) and data from Culberson and Pytkowicz (1968) for deep water at 2 • C at 10000 dbar.Although improving the accuracy of the pressure adjustments to equilibrium constants should be a high priority for future research, our aim here is to assess package precision.

What is significant?
If software packages with identical input cannot agree to within much less than the measurement precision of a computed variable (e.g., pCO 2 ), then their varied use would add substantially to the total uncertainty.To avoid this situation, it is necessary for these tools to have a numerical precision that is far superior to the measurement precision.By numerical precision, we mean their agreement, including all coding differences and errors as well as the usually much smaller numerical round-off error.Therefore, we arbitrarily define the cutoff level for numerical precision to be 10 times smaller than the best measurement uncertainty (Dickson, 2010, Table 1.5).A package that agrees with a given variable from the reference package within the numerical cutoff specified in Table 8 will be referred to here as having a negligible discrepancy relative to the reference; conversely, a package with a greater difference for a given variable will be considered to have a significant discrepancy.

Results
Because the CO2SYS variants agree so closely (Fig. 1), subsequent comparison usually shows results only for CO2SYS-MATLAB (our reference).Packages were compared in terms of how computed variables differed with latitude and depth (using global gridded data) and how individual physical variables and chemical choices affected results (using simplified data).Packages were also compared in terms of computational efficiency.

Global gridded data
In this section, all variables are computed from the GLODAP-WOA2009 gridded input data.With those data, we first compare two new approaches to the common approach of computing subsurface f CO 2 and pCO 2 , in two packages.Then we expand comparison to all packages and other variables.
The mocsy and seacarb packages offer the three approaches to compute f CO 2 and pCO 2 (Sect.2.7).Both packages agree within 0.008% (0.03 µatm at the surface) for each approach for each of the two variables (Tables 9 and  10).While the two packages always compare well throughout the water column, the three approaches diverge as depth increases, as detailed in our companion paper for one package (Orr and Epitalon, 2015, Figs. 1 and 3).Differences between common and potential pCO 2 reach 7 µatm at 5000 m, while differences between potential and in situ pCO 2 are much larger.The latter is 5% greater than the former at 100 m but 18 times larger at 5000 m.Differences between potential and in situ f CO 2 are smaller because they involve pressure corrections only to K 0 and not C f .Yet they still differ by more than a factor of 2 at 5000 m.Subsequent comparison of pCO 2 shows just the common approach, the only one offered by all packages.
More generally, surface zonal means from all packages agree within 0.2 µatm for pCO 2 , 0.006 µmol kg −1 for CO * 2 , 0.0002 units for pH, 0.1 µmol kg −1 for CO 2− 3 , 0.004 for A , and 0.1 for the Revelle factor (Fig. 2).Packages diverge as pressure increases, but agreement generally remain within a factor of 2 of that seen at the surface (Fig. 3).There are two exceptions: the disagreement in pH is 5 times larger  K HS and K NH 4 are the constants for dissociation of hydrogen sulfide and ammonium (Millero, 1995); other constants are defined in the text.at 5000 m, where csys is 0.001 larger than other packages, which agree within 0.0003; for the Revelle factor R f , packages agree within 0.02 throughout the water column except for seacarb, whose discrepancy grows to 0.2 at 5000 m.Although seacarb computes R f with an efficient analytical formula (Frankignoulle, 1994), that approach neglects effects of P T and Si T on total alkalinity, unlike the less efficient numerical approach used in other packages (Orr and Epitalon, 2015).Overall, discrepancies among packages are larger at depth, but they remain negligible (Table 8) except for pH and pCO 2 .Yet agreement was not always so close.For some perspective, the same CO2SYS-MATLAB reference was also compared to older versions of four packages: CO2calc (version 1.0.4revised on 18 June 2013), csys (version revised on 3 February 2010), seacarb (version 2.3.3 revised on 2 April 2010), and an early predecessor of mocsy developed by Orr et al. (2005) but not released publicly.Discrepan-cies relative to the same reference were once larger, e.g., more than 10 times as much for pCO 2 , pH, and CO 2− 3 (Fig. 4).With the mocsy precursor, there are significant discrepancies in pCO 2 reaching up to 1.5 µatm at the surface.Those grow with depth, e.g., reaching 4 µatm at 5000 m.At the same depth, there are discrepancies in CO 2− 3 reaching 0.5 µmol kg −1 and in pH up to 0.007.Subsurface discrepancies are mainly due to two common modeling approximations that were corrected in the first public release of mocsy (Orr and Epitalon, 2014).With CO2calc v1.0.4,surface discrepancies reach up to 2 µatm in pCO 2 , up to 1.3 µmol kg −1 in CO 2− 3 , and up to 0.007 in pH.Those discrepancies are associated with coding errors in the K 1 and K 2 formulations from Lueker et al. (2000), errors that were corrected in CO2calc version 1.2.0.With the previous version of csys, surface pCO 2 is about 1 µatm lower than the reference because that variable was mislabeled; it was actually f CO 2 .As for seacarb v2.3.3, there are no significant discrepancies.However, with an even earlier version of seacarb (v2.0.3 released in 2008, not shown), the only package that maintains public access to all previous versions, discrepancies at depth are much larger (e.g., −7 µmol kg −1 in CO 2− 3 and −0.165 in pH at 4000 m).Because earlier versions of packages often have much larger discrepancies, users would be wise to keep their carbonate system software up to date.
Previous analysis has illustrated how discrepancies vary spatially across the global ocean, but the realistic gridded input data sets that were exploited did not allow us to isolate how discrepancies vary with individual physical variables and chemical input options.We will now focus on those factors, individually, by exploiting simple artificial input data.(Locarnini et al., 2010;Antonov et al., 2010;Garcia et al., 2010).Curves are shown for each package and variable after subtracting off corresponding results for the CO2SYS-MATLAB reference.The csys package does not provide results for A and the Revelle factor.It also neglects nutrient alkalinity, but its curves were adjusted to include the effects of P T and Si T as computed by mocsy.
Figure 3. Global-mean vertical profiles of variables computed from the same gridded data products as in Fig. 2. For each software package, corresponding results from the reference (CO2SYS-MATLAB) have been subtracted.The csys curves are adjusted as in Fig. 2. In all comparisons, the csys results are computed with the option ocdflag = 1; its discrepancies would be larger with ocdflag = 0.  b Following Weiss (1974), f CO 2 = CO * 2 / K 0 exp (1 − P ) vCO 2 /RT , The exponential term vanishes when P is set to 1 atm (common and potential approaches) but is intended to represent total pressure (in situ approach).The potential approach uses θ in place of in situ T (common and in situ approaches) in the above equation and in the calculation of K 0 .a Area-weighted global means from same gridded input data as in Table 9. b Following Weiss (1974), pCO 2 = f CO 2 /C f , where C f = exp B + 2 x 2 2 δ 12 P /RT and P is the total pressure (atmospheric + hydrostatic) as adopted for the in situ approach; the other two approaches assume that P is only atmospheric pressure.The potential approach also uses θ in place of in situ T .

Physical factors
Packages were compared with five common input pairs with the same simple data sets where T , S, and P were varied individually.All packages were compared with the A T -C T pair.Comparison with the four other pairs excluded the mocsy package, which is designed to use only A T -C T .Comparison with two of the pairs, A T -pCO 2 and C T -pCO 2 , excluded the csys package, which does offer pCO 2 as an input variable but only when paired with pH.

A T -C T
With the A T -C T pair, packages agree within 0.2 µatm in pCO 2 , 0.05 µmol kg −1 in CO 2− 3 , and 0.0004 in pH across the observed ranges of ocean T and S at surface pressure (Fig. 5).Surface discrepancies are significant only for one variable from one package, pCO 2 from ODV, but those re-main quite small (less than twice our arbitrary numerical cutoff of 0.1 µatm).Away from the surface, in the open ocean with its cold deep waters at around 2 • C, pressure corrections in all packages do not add significantly to the discrepancies seen at the surface.Yet some deep waters can be warmer, for instance around 13 • C in the Mediterranean Sea.At that temperature, inconsistencies would be more apparent if they were due to errors in coefficients of pressure corrections, which are quadratic functions of temperature (Eqs. 2 and 3).One package, swco2, does indeed exhibit substantial discrepancies with deep water at 13

A T -pH
With the A T -pH input pair (Fig. 6), surface discrepancies between packages remain negligible for all variables.All packages agree within 0.02 µatm in pCO 2 , 0.02 µmol kg −1 in CO 2− 3 , and 0.08 µmol kg −1 in C T across ranges of observed T and S. Below the surface, the swco2 package's subsurface discrepancies remain negligible with the pressure correction at 2 • C, but for water at 13 • C they start to become significant below 4000 db.For the other packages, pressure corrections lead to negligible discrepancies for all variables.

A T -pCO 2
With the A T -pCO 2 input pair (Fig. 7), surface discrepancies are always negligible.The five packages differ by less than 0.05 µmol kg −1 in C T and 0.015 µmol kg −1 in CO 2− 3 .Likewise for pH, packages generally agree within 0.0001; only CO2calc exhibits larger variability (within ±0.0004), but those variations are randomly distributed with a mean near zero, a consequence of CO2calc's limited output precision of only 3 decimal places for pH.The pressure correction when performed at 2 • C does not add significant discrepancies, unlike that performed at 13 • C, for which discrepancies in swco2 grow linearly with pressure, e.g., reaching 3 , and +0.0002 in pH at 5000 db.

C T -pH
With the C T -pH input pair (Fig. 8), there is similar agreement for all packages across ranges of surface T and S. Packages agree within 0.0015 µatm for pCO 2 , 0.007 µmol kg −1 for CO 2− 3 , and 0.1 µmol kg −1 for A T at surface pressure.With C T -pH, unlike with previously analyzed pairs, the swco2 package's pressure corrections do not induce substantial discrepancies in computed subsurface pCO 2 and CO −2 3 , even at 13 • C. Yet swco2 does have significant discrepancies in computed subsurface A T (e.g., 1 µmol kg −1 at 4000 db); conversely, with the pressure correction at T = 2 • C, swco2's A T discrepancies are negligible, consistent with previous patterns.In contrast, there is little temperature sensitivity associated with the slight yet always negligible subsurface discrepancies from ODV.

C T -pCO 2
With the C T -pCO 2 input pair (Fig. 9), all of the five packages have negligible surface discrepancies for computed A T (≤ 0.1 µmol kg −1 ), CO 2− 3 (≤ 0.01 µmol kg −1 ), and pH (≤ 0.003 units).Out of the five packages offering both the C T -pCO 2 and the C T -pH input pairs (excluding csys and mocsy), only swco2 develops significant subsurface discrepancies and only at 13 • C for one variable, in both cases.At that temperature, the swco2 package's discrepancies in computed A T grow linearly with depth, reaching 1 µmol kg −1 at 4000 db, similar to those seen with the  (Fig. 8).As before with the low-temperature correction, discrepancies in swco2's A T remain negligible.
Considering results from the five input pairs together, we can now make several general comments.For all intents and purposes, surface discrepancies remain negligible.The only exception is pCO 2 computed by ODV with the A T -C T input pair, but its discrepancies still remains less than one-fifth of the best measurement precision.Subsurface discrepancies are not significantly worse than those at the surface, i.e., for common deep waters at 2 • C. Conversely, with deep waters at 13 • C, characteristic of the Mediterranean Sea, one package (swco2) does exhibit significant subsurface discrepancies.Yet even under those extreme conditions, swco2 discrepancies above 1000 m remain less than the best measurement precision.They concern either computed A T or other variables computed when A T is a member of the input pair.

Chemical factors
In Sect.3.2, we compared differences among packages while varying physical input for different input pairs.Here we assess differences due to chemical factors, namely accounting for alkalinity from silicic and phosphoric acids (nutrient alkalinity) and opting for potentially important developments since publication of the best-practices guide (Dickson et al., 2007).

Nutrients
Both P T and Si T contribute to the total alkalinity.Thus they affect computed carbonate alkalinity A C when their concentrations are significant and one member of the input pair is A T .One of the packages, csys, neglects this nutrient alkalinity, assuming P T and Si T concentrations are always zero.All other packages account for nutrient alkalinity.Two of those exhibit discrepancies relative to CO2SYS-MATLAB that become significant as nutrient concentrations are increased to the maxima observed in the ocean (Fig. 10).Discrepancies for swco2 grow linearly with nutrient concentrations, reaching −0.2 µatm in pCO 2 , +0.07 µmol kg −1 in CO 2− 3 , and +0.0002 units in pH.These discrepancies are largely associated with Si T ; those from P T are more than 10 times smaller.For CO2calc, discrepancies in pH seem to reach up to nearly 0.001, but those are due to the precision in CO2calc's pH output (given to only three decimal places).
These differences between packages are at least 70 times smaller than the actual changes in computed variables attributable to alkalinity from P T and Si T .With the A T -C T input pair, this nutrient alkalinity increases computed pCO 2 by 6 µatm for average surface waters in the Southern Ocean and by 12 µatm for average deep waters (below 2000 m); simultaneously, CO 2− 3 is reduced by about 2 µmol kg −1 in the same waters (Orr and Epitalon, 2015).

Total boron
Relative to the standard formulation for total boron (Uppström, 1974), the new formulation (Lee et al., 2010) represents about a 3% increase of borate alkalinity throughout the ocean.Hence we first assessed whether or not packages gave consistent responses when changing from the standard to the new formulation.For the six packages that allow for the new formulation (CO2SYS-MATLAB, both CO2SYS-Excel variants, CO2calc, mocsy, and seacarb), computed changes agree within 0.15 µatm for pCO 2 , 0.02 µmol kg −1 for CO 2− 3 , and 0.00006 for pH, i.e., with the A T -C T pair across observed ranges of T , S, and P (Fig. 11).With other input pairs, agreement is closer still, but the comparison is limited to fewer packages (mocsy treats only A T -C T ).Much larger are the actual changes themselves.With the A T -C T pair, given global average surface conditions (T = 18 • C, S = 35), changing from the standard to the new formulation for total boron increases pCO 2 by 5.7 µatm, decreases CO 2− 3 by 2.1 µmol kg −1 , and decreases pH by 0.0056 units.Changes are generally smaller with the A T -pH and A T -pCO 2 pairs (e.g., −0.3 and −0.4 µmol kg −1 for CO 2− 3 , respectively).Conversely, C T and C T -pCO 2 pairs, changes are negligible for all computed carbonate system variables except total alkalinity.

K 1 and K 2
The formulations for K 1 and K 2 from Lueker et al. (2000) are recommended for best practices (Dickson et al., 2007), but they are intended to be restricted to waters with S between 19 and 43 and T between 2 and 35 • C. For waters with physical conditions outside of those ranges, there are no recommended K 1 and K 2 formulations.However, formulations exist, such as those from the most recent reassessment by Millero (2010), which are applicable over wider ranges of S (1-50) and T (0-50 • C).
With an analysis analogous to that shown in Fig. 5 (Sect.3.2.1),we replaced formulations for K 1 and K 2 from Lueker et al. (2000) with those from Millero (2010) to assess the consistency of computed variables in the six packages that include this newer option (CO2SYS-MATLAB, both CO2SYS-Excel variants, CO2calc, seacarb, and mocsy).With the A T -C T pair, four out of six packages agree closely at surface pressure across ranges of T and S (Fig. 12).Conversely, CO2calc differs from the CO2SYS-MATLAB reference by up to −12 µatm in pCO 2 , −1.2 µmol kg −1 in CO 2− 3 , and +0.006 units in pH.Discrepancies are also found for seacarb, reaching up to −20 µatm in pCO 2 , ±0.2 µmol kg −1 in CO 2− 3 , and −0.0025 units in pH.However, seacarb's discrepancies at the surface are inconsistent with its negligible subsurface discrepancies.Pressure corrections alter CO2calc's discrepancies by less than +1 µatm in pCO 2 , −0.05 µmol kg −1 in CO 2− 3 , and −0.001 units in pH, changes that are notably less than its surface discrepancies.Although fewer packages offer the Millero (2010) formulations for K 1 and K 2 , the resulting differences between packages reach levels that are orders of magnitude larger than with the Lueker et al. (2000) formulations.

Computational efficiency
Besides the accuracy and precision of the different packages that compute carbonate system variables, some users with large data sets may be concerned with computational efficiency.To assess differences in computation time between packages, we chose to use the global gridded data set described in Sect.3.1.With nearly 1 million ocean grid points, the computational time needed to compute all carbonate system variables varies by more than a factor of 1800 between packages (Table 2).The slowest package (swco2-Excel) required more than half a day, while the fastest (mocsy) needed 30 s. Packages based on spreadsheets are generally slower than those run by directly calling routines with programming languages (Fortran 95, MATLAB, Visual Basic).The latter are usually coded so that the equilibrium calculations are made one time for each set of input data, whereas spreadsheets often repeat the same set of calculations for each computed variable (each cell).An exception is CO2SYS-Excel-Pelletier, whose execution time rivals that of CO2SYS-MATLAB.Fortunately, even with the slowest of the packages shown in Table 2, the computational time is triv-ial for most observational analysis efforts, because the number of samples is much smaller.Hence developers of most packages have not concerned themselves with computational efficiency.Nonetheless, for very large data sets and for models, which may have millions of grid cells to treat every time step, computational efficiency is critical.

Discussion
To diagnose why computed variables differ between packages, we computed their sensitivities to each constant, assessed errors in individual constants, and used both to assign causes.

Sensitivity to individual constants
A computed variable y is affected by errors in all input variables (including constants as well as members of the input pair), each denoted here as x i (for i = 1, 2, . ..n).Thus, we calculated the sensitivity ratio as the relative change of y to the relative change in each x i , namely ∂y/y : ∂x i /x i .These sensitivity ratios were determined numerically in three successive steps.First, we calculated variables with seacarb under our standard conditions (S = 35, T = 18 • C, P = 0 db, and zero nutrients, along with C T = 2058.185and A T = 2300 µmol kg −1 ).Then, we increased each input variable by 1 % (∂x i /x i = 0.01), individually, and recalculated out- put variables with seacarb for each perturbation.Finally, we took the difference between the first and second computations to obtain the proportional change in the computed variable ∂y/y.

Biogeosciences
Table 11 shows these sensitivity ratios for each variable and constant.Sensitivities of computed variables to input A T , C T , K 0 , and K 1 are much like those from Dickson and Riley (1978), who used the same approach but with different software, input data, and values for constants.Yet our sensitivities to K 2 differ.Although only ∼ 30 % smaller for computed CO * 2 , pCO 2 , and H + , they are 7 times larger for CO 2− 3 as well as A and C .The former moderate reductions occur because the K 2 from Lueker et al. (2000) at standard conditions is 30 % greater than Dickson and Riley's (pK 2 = 9.115).Conversely, our sevenfold-greater sensitivity of CO 2− 3 to K 2 is mysterious.It does not appear to come from our numerical derivative, which we have verified by computing analytical solutions by hand in other pairs when feasible.However, part of the difference could derive from slightly different approaches: Dickson and Riley (1978) used A C as an input variable, while we used A T .
With the A T -C T input pair, the sensitivities to K 1 and K 2 generally dominate, as expected from Dickson and Riley (1978).There is a similarly large sensitivity to K 0 when computing pCO 2 and f CO 2 .The sensitivities to other constants have not been discussed previously.We find a large sensitivity to K B , even surpassing that to K 1 for computed H + , HCO − 3 , CO 2− 3 , A , and C .The sensitivity to K W is also significant but remains 14 to 26 times smaller than that for K B .The sensitivity to other constants remain small (absolute values less than 0.001), except for the solubility products K A and K C , which are inversely proportional to A and C (large negative sensitivities of −1).Sensitivities differ with other input pairs as shown by Dickson and Riley (1978).
These sensitivities are fundamental to the classic propagation of relative errors that has already been applied to the carbonate system (Dickson and Riley, 1978;Dickson et al., 2007).In that, the uncertainty u of computed variable y is expressed as Table 11.Ratio of relative change a,b between output and input variables c, d (∂y/y) / (∂x i /x i ).Variables computed from A T and C T with each package minus corresponding results from CO2SYS-MATLAB, as in Fig. 5 but with K 1 and K 2 from Millero (2010) instead of from Lueker et al. (2000).Shown are computed pCO 2 (top), CO 2− 3 (middle), and pH (bottom) across ranges of T (column 1), S (column 2), and P when T = 2 • C (column 3) and when T = 13 • C (column 4).
where the left-hand side is the relative error in y, a function of the right-hand side's individual relative errors of each input variable and constant (u(x i )/x i ) multiplied by the square of the associated sensitivity term (∂y/y) / (∂x i /x i ) (Table 11).Hence to assess the relative importance of each input variable, we need not only its sensitivity but also its uncertainty.For the case where relative errors for each of the constants are assumed to be similar (Table 8), it is largely the sensitivity term which determines the relative contribution of each constant to the overall error.Yet numerical errors in computed constants are neither identical nor entirely negligible in all packages.

Errors in equilibrium constants
In order to identify sources of error, equilibrium constants were plotted in the same manner as computed variables.By characterizing errors graphically, we were able to use patterns in discrepancies to help isolate problems and eventually identify causes, particularly in packages where source code was available.For packages without source code, we attempted to reproduce discrepancy patterns by making temporary modifications to another package where source code was available.

K 0 , K 1 , and K 2 : best practices
The constants that have the greatest potential to cause the differences in computed variables seen in Fig. 5 are K 0 , K 1 , and K 2 (Table 11), simply because of their prominence in the fundamental equilibria.For K 0 , all packages have negligible discrepancies relative to CO2SYS-MATLAB across ranges of T and S, and all packages agree with the check value from Dickson et al. (2007) to its fourth and final decimal place (Table 12).For K 1 and K 2 with the Lueker et al. (2000) formulation, the story is the same (Fig. 13).None of these three best-practices constants can be responsible for significant discrepancies in computed variables.

Alternative K 1 and K 2 for low salinities
Section 3.3.3detailed the large discrepancies among packages, in terms of computed variables, that were generated simply by replacing the Lueker et al. (2000) formulations for K 1 and K 2 with those from Millero (2010), i.e., the most recent alternative for low-salinity waters and the one based on the greatest number of measurements.Hence the implementation of the Millero (2010) formulations must be done inconsistently between the different packages.Indeed, six packages (three CO2SYS variants and three others) offer that new option, but three of them show significant surface discrepancies in K  a Free scale (all other pKs are on the total scale).
b All CO2calc constants shown here are computed from other output variables; the pK 1 and pK 2 calculated directly by CO2calc, provided by its developers, agree with the reference to beyond six decimal places (L.Robbins, personal communication, 2014).CO2SYS-MATLAB (Fig. 14).Patterns of those discrepancies are qualitatively consistent with patterns of computed variables found with the same option (Fig. 12).A comparison of the source code in three packages revealed that their implemented formulations are strictly identical; however, the sets of coefficients differ.More precisely, Millero fit 551 measurements of K 1 and 590 measurements of K 2 on the seawater pH scale to basic equations of the form where pK 0 i was from his previous fit of the same form for pure water (Millero et al., 2006); T k is the absolute temperature; and A i , B i , and C i are functions of salinity: Millero (2010) provides this set of the six α i coefficients for each of K 1 and K 2 on the original seawater pH scale but also on the free and total scales, analogous fits after converting the measured constants to each of those other scales (Millero, 2010, Eqs. 9 to 12).
With the Millero (2010) formulation, the four packages with internally consistent results are CO2SYS-MATLAB, CO2SYS-Excel-Pierrot, CO2SYS-Excel-Pelletier, and mocsy (Fig. 14).All use the coefficients for his formulation on the seawater scale.In seacarb, however, it is more complicated.At the surface, seacarb uses Millero's set of coefficients on the total scale to compute K 1 and K 2 ; conversely, below the surface seacarb uses Millero's coefficients for the seawater scale and then converts the resulting constants to the total scale after making the pressure correction (as is appropriate).Although the two approaches should yield equivalent results, seacarb's constants computed for the surface differ from those calculated by the reference.Conversely, below the surface, seacarb agrees with the reference.Because seacarb's formulations and coefficients are strictly identical to those published by Millero (2010), which we have verified closely, the surface discrepancies imply an inconsistency between Millero's sets of coefficients for the total and seawater scales.Both sets should yield the same results for K 1 and for K 2 , e.g., once the results from the seawater set are converted back to the total scale.Yet they do not.Indeed when K 1 and K 2 were computed separately from Millero's published sets of coefficients, we found similar patterns of surface discrepancies as already seen between seacarb and CO2SYS-MATLAB (Fig. 15).Patterns matched exactly when we also accounted for the differences in K F .That is, Millero converted pK 1 and pK 2 on the measured  (2000).Constants were computed across the same ranges of T (column 1), S (column 2), and P when T = 2 • C (column 3) and when T = 13 • C (column 4) as used in previous figures.Constants from CO2calc and ODV were not directly available but are estimated from variables computed with the A T -C T pair.The filled black circles indicate the check values (Table 12).Line signatures are as in Fig. 5.  2 and 3, i.e., for the total scale (T , orange dotted line) and the seawater scale (SWS, blue dotted line).The two other sets, also on the T and SWS scales (solid lines), have greater precision, coming from the spreadsheet used for calculations in the same publication (F.J. Millero, personal communication 2013).For consistent comparison, both SWS curves were converted to the total scale using the standard approach (Millero, 2010, Eqs. 11 and 12) and K F from Perez and Fraga (1987).Curves are shown after subtracting values from the preferred formulation (spreadsheet coefficients for the SWS scale converted to the T scale).

Biogeosciences
seawater scale to the total scale using K F from Perez and Fraga (1987), whereas for this study all packages use K F from Dickson and Riley (1979).Thus seacarb shows surface discrepancies (relative to CO2SYS-MATLAB) primarily because it uses Millero's coefficients on the total scale, which are inconsistent with those on the seawater scale, and secondly because the K F used by Millero is inconsistent with that used in this study.
The second package that differs substantially from CO2SYS-MATLAB is CO2calc, but only for K 2 .Although CO2calc's main code was taken from CO2SYS, the CO2calc developers included the Millero (2010) K 1 and K 2 formulations themselves, before they were available in CO2SYS.Lacking the CO2calc source code, we studied discrepancy patterns and made sensitivity tests.These differences appear to come from a different number of significant figures in one of the coefficients.More precisely, Fig. 15 compares the constants computed from the published set of coefficients (Millero, 2010, Tables 2 and 3) to those computed with the unpublished yet more precise set of coefficients used by Millero, i.e., his spreadsheet for the same publication (F.J. Millero, personal communication 2013).Only the α 5 coeffi- Figure 16.The ratio of K 1 /K 2 for each set of coefficients from Millero (2010) over the same ratio for the reference.We refer to that ratio of ratios as the K 1 /K 2 quotient.The reference is chosen arbitrarily as Millero's set of coefficients on SWS scale from his spreadsheet.Curve colors and line patterns are as in Fig. 15.cient differs, being given with one additional decimal place in the spreadsheet (Table 13).The difference between constants computed with the published coefficients and those computed with unpublished coefficients, i.e., the spreadsheet coefficients on the seawater scale, match the pattern and magnitude of the differences between CO2calc and CO2SYS-MATLAB (compare Figs. 15 and Fig. 14).Hence CO2calc developers appear to have used the set of seawater-scale coefficients from Millero's spreadsheet.
To further confirm that the different sets of coefficients yield fundamentally different results, we made two additional comparisons.First, we compared their K 1 /K 2 ratios (Fig. 16).By definition, that ratio should be independent of the pH scale, yet with different sets of coefficients it differs by up to 1.5% over the observed range of T and S. Secondly, we made a sensitivity test with CO2SYS-Excel-Pelletier's two options for the Millero (2010) formulation (two sets of coefficients), both of which became available after our discussion paper was published.With that package, users choose to use either the published seawater-scale coefficients (as in the CO2SYS-MATLAB) or the unpublished seawater-scale coefficients from the spreadsheet (as used by CO2calc).With the former, CO2SYS-Excel-Pelletier agrees with the reference; with the latter, results are like those for CO2calc.Therefore, all three simple analyses indicate that it  12).

Biogeosciences
is the different sets of coefficients that cause substantial differences in computed variables.
It is tempting to conclude that all package developers should simply adopt the spreadsheet's more precise set of coefficients on the seawater scale, given their greater precision and their consistency with original results (Millero, 2010, Fig. 3).However, even that set of unpublished coefficients, which is more precise than those published, may have two few significant figures To test that concern, we exploited the same spreadsheet's coefficients for the total scale, which are given to more decimal places than its seawater-scale coefficients.By incrementally reducing the number of decimal places in each unpublished total-scale coefficient (one at a time), we determined the point at which calculated variables change significantly.Thus we found that the published α 1 for K 1 for the total scale would need to be extended from 4 to 5 decimal places, while the α 5 for K 1 as well as that for K 2 would both need to be extended from 4 to 6 decimal places before results match those computed with the unpublished spreadsheet's total-scale coefficients.It follows that the published seawater-scale coefficients should be extended likewise, because they have similar magnitudes and the same number of significant figures as the published total-scale coefficients.The much larger disagreement among packages found when changing from the Lueker et al. (2000) formulations for K 1 and K 2 to those from Millero (2010) emphasizes the danger of applying conclusions from one comparison to cases with different sets of constants.

K B and K W : principal non-carbonate alkalinity constants
Previous comparison revealed significant discrepancies in subsurface variables computed from the swco2 package (Fig. 5-9), which are not due to K 1 and K 2 (Fig. 13).These discrepancies occur only when A T is a member of the input pair or when A T is computed, suggesting that they stem from the need to correct from total to carbonate alkalinity.The largest factor to correct for is borate alkalinity; hence it is also the most likely cause.Indeed, comparison of K B computed by the different packages does reveal discrepancies for swco2 (Fig. 17).Furthermore, swco2's divergence from CO2SYS-MATLAB increases linearly with depth, consistent with discrepancies in computed variables.Although we do not have access to the code for swco2, this discrepancy is consistent with a sign error in its a 2 pressure-correction coefficient for K B , as verified by sensitivity tests in seacarb.Figure 18.The pK F (top), pK S (middle), and pK A (bottom) computed from each package minus corresponding values from CO2SYS-MATLAB.The formulation for K F is from Dickson and Riley (1979) as recommended by Dickson and Goyet (1994), with all packages on the free scale.The formulation for K S is from Dickson (1990a) and on the free scale, while that for K A is from Mucci (1983) with no scale, as recommended for best practices.The filled black circles indicate check values (Table 12).
2 orders of magnitude smaller and become negligible when carried out at 2 • C rather than 13 • C.
Regarding K W , none of the packages diverges significantly from CO2SYS-MATLAB.Although CO2calc's K W appears to oscillate about the reference, the extremes of those oscillations remain negligible.Moreover, our estimates of CO2calc's equilibrium constants are not precise.Having to compute them from output variables with limited precision, we were unable to estimate CO2calc's equilibrium constants typically beyond the third place after the decimal.

K F and K S : constants to change pH scales
For K F , all packages agree with the check value to its third and final digit after the decimal.The K F from the swco2 package does diverge from the others as salinity increases, but it still agrees with CO2SYS-MATLAB within pK F = 0.00006, and its discrepancies do not change with temperature or pressure (Fig. 18).
For K S , all packages agree with the check value, again to its third and final digit after the decimal.Beyond that, packages agree even more tightly, with the largest divergence reaching pK S = 0.0002 for swco2.There are visible differences for K S computed from the different packages that merit further investigation (e.g., positive excursions under surface conditions for swco2), but they remain quite small.Given the negligible consequences, we leave their resolution to future work.

K A and K C : solubility products
For K A , the situation is similar at the surface but not at depth (Fig. 18).Under surface conditions, no packages have significant discrepancies.Although CO2calc appears to have discrepancies larger than other packages (average pK A ∼ 0.0002), they remain quite small; moreover, they may be exaggerated because we had to calculate them from computed variables with limited output precision (2 decimal places for CO 2− 3 ).Pressure-correction discrepancies are negligible in all packages but one, swco2.At 5000 db, the discrepancy for 13 • C water reaches pK A = 0.015, thereby biasing K A to be 3.4 % too low and A to be 3.4 % too high.However, for more typical deep waters at 2 • C, those discrepancies are reduced by a factor of 7. The form of the swco2 discrepancy curve, quadratic with pressure, suggests an error in the b 1 pressure-correction coefficient.Without access to the swco2 source code, we made sensitivity tests with seacarb.Inversing the sign of its b 1 coefficient reproduced the form and magnitude of the swco2 discrepancies.Thus it appears that swco2's b 1 coefficient for K A needs to be changed to +0.0003692.Other swco2 pressure-correction coefficients for K A appear to be correct.12).
4.2.6 K 1P , K 2P , and K 3P : constants for phosphoric and silicic acids Constants for phosphoric and silicic acids enter into the calculations only when nutrient concentrations are significant and A T is a member of the input pair.Under those conditions, A T must be corrected for nutrient alkalinity to provide an accurate estimate of A C , as needed to compute other variables.For K 1P , K 2P , and K 3P , three packages (seacarb, mocsy, and CO2SYS-Excel-Pelletier) agree with results from CO2SYS-MATLAB across ranges of T , S, and P (Fig. 19).Three other packages (CO2SYS-Excel-Pierrot, CO2calc, and ODV) do not provide these constants as output, nor could they be calculated from available variables.Yet they exhibit negligible discrepancies for computed variables as a function of nutrient concentrations (Fig. 10), suggesting that their discrepancies in associated equilibrium constants must also be negligible.
The remaining package, swco2, differs significantly from CO2SYS-MATLAB, with a constant shift of 0.006 for each of K 1P , K 2P , and K 3P across ranges of T and S. Sensitivity tests with seacarb suggest that this constant shift is caused by swco2 making the necessary conversion from the seawater pH scale to the total pH scale, but doing so twice.The original formulations of these constants are from Millero (1995) and are on the seawater scale.For their conversion to the total scale, the best-practices approach simply subtracts 0.015 (Dickson and Goyet, 1994;Dickson et al., 2007, Chap. 5, footnote 5), i.e., a constant correction.Conversely, the more rigorous approach to convert between those two pH scales (e.g., Millero, 2010, Eq. 6) results in an offset that varies with the hydrogen fluoride concentration [HF].For example, with K F from Perez and Fraga (1987) the offset ranges from 0 to 0.032 across observed ocean T and S; with K F from Dickson and Riley (1979), it ranges from 0 to 0.024.Our tests suggest that all packages make the variable correction but only swco2 does not first remove the 0.015 offset (equivalent to a 0.006 shift in pK).
For K Si , only the swco2 package reveals any discrepancies relative to CO2SYS-MATLAB (Fig. 17).Out of the other packages, three agree with K Si from CO2SYS-MATLAB (seacarb, mocsy, and CO2SYS-Excel-Pelletier), three do not provide K Si as output (CO2calc, ODV, and CO2SYS-Excel-Pierrot), and one does not compute K Si (csys).Discrepancies for K Si in swco2 under surface conditions are identical to those for its K 1P , K 2P , and K 3P .The constant positive excursion of 0.006 appears due to the same cause, correcting the equilibrium constant from the seawater to the total scale two times.Below the surface, swco2's discrepancy for K Si grows linearly with pressure, just as does its discrepancy for K B (Sect.4.2.3).And the cause appears to be identical, a sign error in the a 2 pressure correction coefficient, based on our sensitivity tests in seacarb.Thus we recommend that swco2's a 2 coefficient should be checked and changed if necessary to −0.002608.Just as for K B , discrepancies in swco2's pres- sure correction for K Si are lower when carried out at 2 • C rather than at 13 • C. The total discrepancies in pK Si are 0.007 and 0.015, respectively (at 4000 db).That implies that pressure-correction discrepancies are about 10 times larger at 13 • C than at 2 • C, after removing the 0.006 constant offset at surface conditions.

Conclusions
To assess the consistency of carbonate chemistry software packages, we have compared computed variables from 10 publicly available distributions, identifying discrepancies and causes.This comparison has led to improved agreement.Since our discussion paper was published (Orr et al., 2014), there has been a fivefold reduction in discrepancies in pCO 2 and CO 2− 3 when all packages use the set of constants recommended for best practices (Dickson et al., 2007).The small discrepancies that currently exist remain insignificant even after pressure adjustments are made for the cold waters that pervade the deep ocean.Only in warm deep waters, such as found in the Mediterranean Sea, are there significant discrepancies and only for one package (swco2), e.g., when pCO 2 and CO 2− 3 are computed from the A T -C T input pair.Those discrepancies appear to derive from a sign error in a pressurecorrection coefficient for K B .Similar sign errors in pressurecorrection coefficients for swco2's K A and K C appear to cause its A and C to be underestimated by 3 % at 4000 db when at 13 • C but only 0.5 % when at 2 • C.
The choice of equilibrium constants affects package agreement.Despite the excellent agreement found when packages use the best-practices set of equilibrium constants, their accuracy under extreme conditions is questionable.Best-practices formulations for K 1 and K 2 (Lueker et al., 2000) are based on measurements that did not include conditions such as found in estuaries (S < 19) nor in the ocean's coldest waters (T < 2 • C), which comprise 11% of its global surface area and 42% of its global volume, based on an annual climatology (Locarnini et al., 2010).Thus we also compared packages changing only the formulations for K 1 and K 2 to those that consider low-salinity and low-temperature waters (Millero, 2010).Out of the six packages where that newer option is available, three agreed with the reference, while three others differed, e.g., by up to 7 µatm in surface pCO 2 .One package differs because it uses the set of published coefficients to compute the salinity dependence of K 1 and K 2 on the total scale; conversely, the reference uses another set to compute those constants on the seawater scale, later converting to the total scale.Their disagreement after conversion indicates a fundamental inconsistency between the two published sets of coefficients.Two other packages differ because one coefficient, also for the seawater scale, has an additional significant figure taken from an unpublished spreadsheet.Other published may lack up to 2 significant figures based on our tests with the same spreadsheet's coefficients on the total scale, which have still greater precision.These discrepancies emphasize the fundamental need for new measurements of K 1 and K 2 at low salinities, low temperatures, and high pressures.
To limit future inconsistencies, we offer several practical recommendations.Users are encouraged to use up-to-date software, to use the set of constants recommended for best practices (Dickson et al., 2007), and to avoid the K 1 and K 2 from Millero (2010) until discrepancies are resolved.For now, users are also advised to avoid the new total boron-tosalinity ratio (Lee et al., 2010) and favor the "best-practices" ratio (Uppström, 1974), which was used to compute K 1 and K 2 from laboratory measurements (Mehrbach et al., 1973, Eq. 8).For reproducibility, users should cite not only the package name but also its version number and all chosen equilibrium constants as well as any other selected package options.Developers would facilitate future comparison and reproducibility by providing computed constants as output, releasing their source code, and offering access to older versions, ideally through a public revision control system.Developers are also encouraged to consider providing results for potential and in situ pCO 2 , which differ greatly in deep waters.Although our focus has been on public packages, it is just as necessary to validate privately developed code.To facilitate such validation, we provide a supplement containing a subset of the data produced for this comparison.
The Supplement related to this article is available online at doi:10.5194/bg-12-1483-2015-supplement.

Figure 1 .
Figure 1.Differences ( ) between the variants of CO2SYS relative to the reference MATLAB code for variables computed from A T and C T .Differences are shown across ranges of T (left), S (center), and P (right) for pCO 2 (top), CO 2−3 (middle), and pH (bottom).The three most recent variants (MATLAB and both Excel versions) are run with constants recommended for best practices (BP).The QBasic variant does not offer the same K 1 and K 2 , so we used an earlier refit byDickson and Millero (1987) of the same data (DM87) and compared it to one Excel version also with DM87.

Figure 5 .
Figure 5. Variables computed from A T and C T for each package minus corresponding results from CO2SYS-MATLAB.The computed pCO 2 (top), CO 2−3 (middle), and pH (bottom) are shown across ranges of T (column 1), S (column 2), and P when T = 2 • C (column 3) and when T = 13 • C (column 4).For each range, there is one curve per package and per variable.

Figure 9 .
Figure 9. Variables computed from C T and pCO 2 with each package minus corresponding results from CO2SYS-MATLAB.Shown are computed A T (top), CO 2− 3 (middle), and pH (bottom) across ranges of T (column 1), S (column 2), and P when T = 2 • C (column 3) and when T = 13 • C (column 4) for each package.

Figure 10 .
Figure10.Effect of nutrients on variables computed from A T and C T with each package minus results for CO2SYS-MATLAB.Shown are effects on computed pCO 2 (top), CO 2− 3 (middle), and pH (bottom) across the observed oceanic ranges of P T (right), Si T (center), and their combined effect (left) for each package.Results from the two CO2SYS Excel variants are not shown, but both agree with the reference.Results for csys are not included as it assumes that nutrient concentrations are always zero.

Figure 11 .
Figure 11.Effect of increased total boron on variables computed from A T and C T with each package minus results for CO2SYS-MATLAB.Shown are the effects of the increased boron (Lee et al. (2010) minus Uppström (1974)) on computed pCO 2 (top), CO 2− 3 (middle), and pH (bottom) across ranges of T (left), S (center), and P (right).The Lee et al. (2010) formulation is included in six packages: CO2SYS-MATLAB, CO2calc, seacarb, mocsy, and both CO2SYS Excel variants.The latter two are not shown but agree with the reference.The recommendation to use the Uppström (1974) formulation by Dickson et al. (2007) came before the Lee et al. (2010) study.
1 and K 2 relative to those computed by www

Figure 13 .
Figure13.The pK 0 (top), pK 1 (middle), and pK 2 (bottom) computed with each package minus corresponding values for CO2SYS-MATLAB.Formulations are those recommended for best practices, namely K 0 fromWeiss (1974) and K 1 and K 2 fromLueker et al. (2000).Constants were computed across the same ranges of T (column 1), S (column 2), and P when T = 2 • C (column 3) and when T = 13 • C (column 4) as used in previous figures.Constants from CO2calc and ODV were not directly available but are estimated from variables computed with the A T -C T pair.The filled black circles indicate the check values (Table12).Line signatures are as in Fig.5.

Figure 17 .
Figure 17.The pK B (top), pK W (middle), and pK Si (bottom) computed from each package minus corresponding values from CO2SYS-MATLAB.Formulations are from Millero (1995), as recommended for best practices.The filled black circles indicate the check values (Table12).

Figure 19 .
Figure19.The pK 1P (top), pK 2P (middle), and pK 3P (bottom) computed from each package minus corresponding values from CO2SYS-MATLAB.Formulations are fromMillero (1995), as recommended for best practices.Constants are computed across the same ranges of T (column 1), S (column 2), and P when T = 2 • C (column 3) and when T = 13 • C (column 4) as used in previous figures.The filled black circles indicate the check values (Table12).

Table 1 .
Carbonate system software packages.

Table 2 .
Computational time required to process the GLODAP-WOA2009 gridded data product a .
a Time required in minutes to treat 958 557 ocean grid points.b Time does not include calculation of the Revelle factor.c Time for code run under octave; it may run faster under MATLAB.d

Table 3 .
Operating system and code details for each package.

Table 4 .
Available input pairs for each package.

Table 5 .
Available pH scales for each package.

Table 6 .
Available constants for each package.

Biogeosciences, 12, 1483-1510, 2015 www.biogeosciences.net/12/1483/2015/ ages
while separating the effects of physical input variables (T , S, and P ) on computed variables.For that, we started with five commonly used input pairs:A T -C T , A T -pH, A T -pCO 2 , C T -pCO 2 ,and C T -pH.Then for each pair, we computed the other carbonate system variables over ranges of T , S, and P , assuming zero nutrient concentrations.

Table 7 .
Coefficients used in Eqs.(2) and (3) to correct for effect of pressure on equilibrium constants.

Table 8 .
Desired measurement and numerical uncertainties.

Table 9 .
Oceanic f CO 2 a from three approaches b in two packages.
a Area-weighted global means computed from the GLODAP-WOA2009 gridded data set.

Table 10 .
Oceanic pCO 2 a from three approaches b in two packages.
Lueker et al. (2000)K 2 (bottom) computed in each of four packages minus corresponding values for CO2SYS-MATLAB, as in Fig.13except that formulations are fromMillero (2010), notLueker et al. (2000).The filled black circles indicate the check values (Table12).Six packages include this option.Only CO2SYS-MATLAB, CO2calc, seacarb, and mocsy are shown.Both CO2SYS Excel variants agree with the reference, but when CO2SYS-Excel-Pelletier is switched from the published to the unpublished seawater-scale coefficients, it resembles CO2calc.