# Comparison of ten packages that compute ocean carbonate chemistry

J. C. Orr, J.-M. Epitalon, Jean-Pierre Gattuso

<figure><img src="image_1.png" /></figure>

## To cite this version:

J. C. Orr, J.-M. Epitalon, Jean-Pierre Gattuso. Comparison of ten packages that compute ocean carbonate chemistry. Biogeosciences, 2015, 12 (5), pp.1483-1510. 10.5194/bg-12-1483-2015. hal-01207651

HAL Id: hal-01207651

https://hal.sorbonne-universite.fr/hal-01207651

Submitted on 1 Oct 2015

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_3.png" /></figure>

<figure><img src="image_4.png" /></figure>

# Comparison of ten packages that compute ocean carbonate chemistry

J. C. Orr<fnref n="1" />, J.-M. Epitalon<fnref n="2" />, and J.-P. Gattuso<fnref n="3" />,<fnref n="4" />

*Correspondence to:* J. C. Orr (james.orr@lsce.ipsl.fr)

Received: 24 March 2014 – Published in Biogeosciences Discuss.: 4 April 2014

Revised: 24 January 2015 – Accepted: 31 January 2015 – Published: 9 March 2015

**Abstract.** 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 pCO₂, 0.0002 units for pH, and 0.1 µmol kg⁻<fnref n="1" /> for CO₃<fnref n="2" />⁻ 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 variable to changes in each constant reveals the general dominance of K₁ and K₂ but also the comparable sensitivity to K_B for the A_T - C_T input pair. Best-practice formulations for K₁ and K₂ are implemented consistently among packages. Yet with more recent formulations designed to cover a wider range of salinity, packages disagree by up to 8 µatm in pCO₂, 0.006 units in pH, and 1 µmol kg⁻<fnref n="1" /> in CO₃<fnref n="2" />⁻ under typical surface conditions. They use different proposed sets of coefficients for these formulations, all of which are inconsistent. Users would do well to use up-to-date versions of packages and the constants recommended for best practices.

## 1 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*<sub>T</sub> and total alkalinity *A*<sub>T</sub>, and then compute from associated thermodynamics the corresponding pH, partial pressure of carbon dioxide *pCO*₂, concentrations of aqueous CO₃<fnref n="2" />⁻ as well as carbonate CO₃<fnref n="2" />⁻ and bicarbonate HCO₃⁻ ions, and the related Revelle factor and saturation states of aragonite Ω<sub>A</sub> and calcite Ω<sub>C</sub>. It is the CO₂-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₂ (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

<footnote n="1">¹LSCE/IPSL, Laboratoire des Sciences du Climat et de l'Environnement, CEA-CNRS-UVSQ, Gif-sur-Yvette, France</footnote>

<footnote n="2">²Geoscientific Programming Services, Fanjeaux, France</footnote>

<footnote n="3">³CNRS-INSU, Laboratoire d'Océanographie de Villefranche, Villefranche-sur-mer, France</footnote>

<footnote n="4">⁴Sorbonne Universités, UPMC Univ. Paris 06, Observatoire Océanologique, Villefranche-sur-mer, France</footnote>
was publicly available. Provided with the same input, computed output from the three packages differed by 21 µmatt for pCO₂, 0.16 units for pH, and 15 µmol kg⁻¹ for CO₃²⁻ as well as HCO₃⁻. Packages used different pH scales, different formulations for some of the constants (K₁, K₂, K<sub>B</sub>, and K<sub>S</sub>), 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₂ 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₂ and air-sea CO₂ fluxes from simulated A<sub>T</sub> and C<sub>T</sub> (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₂ and that computed from A<sub>T</sub> and C<sub>T</sub>, 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.

## 2 Methods

To compare all publicly available packages, our approach was to install them in one location, define common input data and constants, and use those with each package to generate a data archive for centralized analysis.
<table><caption>**Table 1.** Carbonate system software packages.</caption><thead><tr><th>Package</th><th>Language</th><th>Version</th><th>Reference</th></tr></thead><tbody><tr><td>CO2SYS<sup>a</sup></td><td>QBasic</td><td>1.05</td><td>Lewis and Wallace (1998)</td></tr><tr><td>CO2SYS<sup>b</sup></td><td>Excel</td><td>24</td><td>Pelletier et al. (2007)</td></tr><tr><td>CO2SYS<sup>a</sup></td><td>Excel</td><td>2.1</td><td>Pierrot et al. (2006)</td></tr><tr><td>CO2SYS<sup>a</sup></td><td>MATLAB</td><td>1.1</td><td>van Heuven et al. (2011)</td></tr><tr><td>CO2calc<sup>c</sup></td><td>Visual Basic</td><td>1.3.0</td><td>Robbins et al. (2010)</td></tr><tr><td>csys<sup>d</sup></td><td>MATLAB</td><td>04–2014</td><td>Zeebe and Wolf-Gladrow (2001)</td></tr><tr><td>ODV<sup>e</sup></td><td>C++</td><td>4.5.0</td><td>Schlitzer (2002)</td></tr><tr><td>mocsy<sup>f</sup></td><td>Fortran 95</td><td>2.0</td><td>Orr and Epitalon (2015)</td></tr><tr><td>seacarb<sup>g</sup></td><td>R</td><td>3.0.6</td><td>Gattuso et al. (2015)</td></tr><tr><td>swco2<sup>h</sup></td><td>Excel</td><td>2</td><td>Hunter (2007); Mosley et al. (2010)</td></tr><tr><td>swco2<sup>h</sup></td><td>Visual Basic</td><td>2</td><td>Hunter (2007)</td></tr></tbody></table>

a http://cdiac.ornl.gov/oceans/co2rprt.html

b http://www.ecy.wa.gov/programs/eap/models.html

c http://pubs.usgs.gov/of/2010/1280/

d http://www.soest.hawaii.edu

e http://odv.awi.de/

f http://ocmip5.ipsl.jussieu.fr/mocsy

g http://cran.r-project.org/package=seacarb

h http://neon.otago.ac.nz/research/mfc/people/keith_hunter/software/swco2/

## 2.1 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

<table><caption>**Table 2.** Computational time required to process the GLODAP-WOA2009 gridded data product<sup>a</sup>.</caption><thead><tr><th>Package</th><th>Total time</th><th>Run time</th><th>Write time</th></tr></thead><tbody><tr><td>swco2 (Excel)</td><td>897.1</td><td></td><td></td></tr><tr><td>CO2calc</td><td>91.2</td><td></td><td></td></tr><tr><td>ODV</td><td>73.3</td><td></td><td></td></tr><tr><td>seacarb (R)<sup>b</sup></td><td>58.3</td><td></td><td></td></tr><tr><td>swco2 (Visual Basic)</td><td>28.0</td><td></td><td></td></tr><tr><td>csys (MATLAB)<sup>c</sup></td><td>7.7</td><td></td><td></td></tr><tr><td>CO2SYS (Excel-Pelletier)<sup>d</sup></td><td>7.2</td><td></td><td></td></tr><tr><td>CO2SYS (MATLAB)<sup>c</sup></td><td>5.9</td><td>5.5</td><td>0.4</td></tr><tr><td>mocsy (Fortran 95)</td><td>0.5</td><td>0.2</td><td>0.3</td></tr></tbody></table>

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 CO2SYS-Excel-Pelletier run on a 2.2 GHz Intel Core 2 duo T7500 under Windows 7 and Excel 2007.

e All other packages run on a 2.0 GHz Intel Pentium dual-core T4200 under Linux or Windows Vista.

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.

## 2.2 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
<figure><img src="image_5.png" /><figcaption>**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_3^{2-}$ (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 by Dickson and Millero (1987) of the same data (DM87) and compared it to one Excel version also with DM87.</figcaption></figure>

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 $fCO_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_3^{2-}$. 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. The seacarb package provides explicit functions to the user to allow conver-
<table><caption>**Table 3.** Operating system and code details for each package.</caption><thead><tr><th rowspan="2">OS & details</th><th colspan="6">CO2SYS</th><th rowspan="2">mocsy</th></tr><tr><th>QBasic</th><th>Excel<sup>a</sup></th><th>Matlab</th><th>CO2calc</th><th>ODV</th><th>csys</th></tr></thead><tbody><tr><td>Linux/Unix</td><td></td><td></td><td>•</td><td></td><td>•</td><td>•</td><td>•</td></tr><tr><td>Windows</td><td>•</td><td>•</td><td>•</td><td>•</td><td>•</td><td>•</td><td>•</td></tr><tr><td>Mac OS</td><td></td><td>•</td><td>•</td><td>•</td><td>•</td><td>•</td><td>•</td></tr><tr><td>iOS</td><td></td><td></td><td></td><td>•</td><td></td><td></td><td></td></tr><tr><td>Public source code</td><td>•</td><td>•</td><td>•</td><td></td><td>•</td><td>•</td><td>•</td></tr><tr><td>User programmable</td><td></td><td></td><td>•</td><td></td><td>•</td><td>•</td><td>•</td></tr><tr><td>Software platform</td><td>E<sup>d</sup></td><td>M<sup>b</sup></td><td></td><td>M<sup>b</sup></td><td>R<sup>f, h</sup></td><td>E<sup>e</sup></td><td>F<sup>g, h</sup></td></tr></tbody></table>

<sup>a</sup> Both variants: CO2SYS-Excel-Pierrot and CO2SYS-Excel-Pelletier.

<sup>b</sup> Package runs under MATLAB (commercial software) or octave (free software).

<sup>c</sup> Spreadsheet interface is not code; core library is callable (Visual Basic) but not modifiable.

<sup>d</sup> Package runs under Excel.

<sup>e</sup> Package runs under Excel (commercial) or LibreOffice (free and open source).

<sup>f</sup> Package runs under R.

<sup>g</sup> Fortran 95 code.

<sup>h</sup> Also runs under Python.

<table><caption>**Table 4.** Available input pairs for each package.</caption><thead><tr><th rowspan="2">Pair</th><th colspan="6">CO2SYS<sup>b</sup></th><th rowspan="2">mocsy</th></tr><tr><th>QBasic</th><th>Excel<sup>a</sup></th><th>Matlab</th><th>CO2calc<sup>b</sup></th><th>ODV<sup>b</sup></th><th>csys</th></tr></thead><tbody><tr><td>A<sub>T</sub>-C<sub>T</sub></td><td>•</td><td>•</td><td>•</td><td>•</td><td>•</td><td>•</td><td>•</td></tr><tr><td>A<sub>T</sub>-pCO<sub>2</sub></td><td>•</td><td>•</td><td>•</td><td>•</td><td>•</td><td>•</td><td>•</td></tr><tr><td>A<sub>T</sub>-pH</td><td>•</td><td>•</td><td>•</td><td>•</td><td>•</td><td>•</td><td>•</td></tr><tr><td>A<sub>T</sub>-CO<sub>3</sub><sup>2-</sup></td><td></td><td></td><td></td><td></td><td>•</td><td>•</td><td>•</td></tr><tr><td>A<sub>T</sub>-CO<sub>2</sub><sup>b</sup></td><td></td><td></td><td></td><td></td><td>•</td><td>•</td><td>•</td></tr><tr><td>A<sub>T</sub>-HCO<sub>3</sub><sup>-</sup></td><td></td><td></td><td></td><td></td><td>•</td><td>•</td><td>•</td></tr><tr><td>C<sub>T</sub>-pCO<sub>2</sub></td><td>•</td><td>•</td><td>•</td><td>•</td><td>•</td><td>•</td><td>•</td></tr><tr><td>C<sub>T</sub>-pH</td><td>•</td><td>•</td><td>•</td><td>•</td><td>•</td><td>•</td><td>•</td></tr><tr><td>C<sub>T</sub>-CO<sub>3</sub><sup>2-</sup></td><td></td><td></td><td></td><td></td><td>•</td><td>•</td><td>•</td></tr><tr><td>C<sub>T</sub>-CO<sub>2</sub><sup>b</sup></td><td></td><td></td><td></td><td></td><td>•</td><td>•</td><td>•</td></tr><tr><td>C<sub>T</sub>-HCO<sub>3</sub><sup>-</sup></td><td></td><td></td><td></td><td></td><td>•</td><td>•</td><td>•</td></tr><tr><td>pCO<sub>2</sub>-pH</td><td>•</td><td>•</td><td>•</td><td>•</td><td>•</td><td>•</td><td>•</td></tr><tr><td>pCO<sub>2</sub>-CO<sub>3</sub><sup>2-</sup></td><td></td><td></td><td></td><td></td><td></td><td>•</td><td>•</td></tr><tr><td>pCO<sub>2</sub>-HCO<sub>3</sub><sup>-</sup></td><td></td><td></td><td></td><td></td><td></td><td>•</td><td>•</td></tr><tr><td>pH-CO<sub>3</sub><sup>2-</sup></td><td></td><td></td><td></td><td></td><td>•</td><td>•</td><td>•</td></tr><tr><td>pH-CO<sub>2</sub><sup>b</sup></td><td></td><td></td><td></td><td></td><td>•</td><td>•</td><td>•</td></tr><tr><td>pH-HCO<sub>3</sub><sup>-</sup></td><td></td><td></td><td></td><td></td><td>•</td><td>•</td><td>•</td></tr><tr><td>CO<sub>3</sub><sup>2-</sup>-CO<sub>2</sub><sup>b</sup></td><td></td><td></td><td></td><td></td><td>•</td><td>•</td><td>•</td></tr><tr><td>CO<sub>3</sub><sup>2-</sup>-HCO<sub>3</sub><sup>-</sup></td><td></td><td></td><td></td><td></td><td>•</td><td>•</td><td>•</td></tr><tr><td>CO<sub>2</sub><sup>b</sup>-HCO<sub>3</sub><sup>-</sup></td><td></td><td></td><td></td><td></td><td>•</td><td>•</td><td>•</td></tr></tbody></table>

<sup>a</sup> Both variants: CO2SYS-Excel-Pierrot and CO2SYS-Excel-Pelletier.

<sup>b</sup> CO2SYS, CO2calc, and ODV also allow input pairs containing fCO₂ instead of pCO₂.

<sup>c</sup> seacarb and swco2 include user-callable functions to convert between pCO₂ and fCO₂.

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 fCO₂ and pCO₂ following Weiss (1974); other packages neglect these pressure effects.
<table><caption>**Table 5.** Available pH scales for each package.</caption><thead><tr><th rowspan="2">pH scale</th><th colspan="7">CO2SYS</th></tr><tr><th>QBasic</th><th>Excel<sup>a</sup></th><th>Matlab</th><th>CO2calc</th><th>ODV</th><th>csys</th><th>seacarb<sup>c</sup></th><th>swco2<sup>c</sup></th><th>mocsy</th></tr></thead><tbody><tr><td>NBS</td><td>•</td><td>•</td><td>•</td><td>•</td><td></td><td></td><td></td><td>•</td><td></td></tr><tr><td>Free</td><td>•</td><td>•</td><td>•</td><td>•</td><td>•</td><td>•</td><td>•</td><td>•</td><td>•</td></tr><tr><td>Total</td><td>•</td><td>•</td><td>•</td><td>•</td><td>•</td><td>•</td><td>•</td><td>•</td><td>•</td></tr><tr><td>Seawater</td><td>•</td><td>•</td><td>•</td><td>•</td><td>•</td><td>•</td><td>•</td><td>•</td><td>•</td></tr><tr><td>Convert pH between scales</td><td></td><td></td><td>•</td><td></td><td></td><td></td><td>•</td><td>•</td><td></td></tr><tr><td>Convert K's between scales</td><td>•</td><td></td><td></td><td></td><td></td><td>•</td><td></td><td></td><td></td></tr></tbody></table>

<table><caption>**Table 6.** Available constants for each package.</caption><thead><tr><th rowspan="2">Constant</th><th colspan="7">CO2SYS</th></tr><tr><th>QBasic</th><th>Excel<sup>a</sup></th><th>Matlab</th><th>CO2calc</th><th>ODV</th><th>csys</th><th>seacarb</th><th>swco2</th><th>mocsy</th></tr></thead><tbody><tr><td>K<sub>1</sub> and K<sub>2</sub></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>(Lueker et al., 2000)</td><td></td><td>•</td><td>•</td><td>•</td><td>•</td><td>•</td><td>•</td><td>•</td><td>•</td></tr><tr><td>(Roy et al., 1993)</td><td>•</td><td>•</td><td>•</td><td>•</td><td>•</td><td>•</td><td>•</td><td>•</td><td>•</td></tr><tr><td>(Goyet and Poisson, 1989)</td><td>•</td><td>•</td><td>•</td><td>•</td><td></td><td></td><td></td><td>•</td><td></td></tr><tr><td>(Hansson, 1973a, b)<sup>b</sup></td><td>•</td><td>•</td><td>•</td><td>•</td><td></td><td></td><td></td><td>•</td><td></td></tr><tr><td>(Mehrbach et al., 1973)<sup>b</sup></td><td>•</td><td>•</td><td>•</td><td>•</td><td></td><td></td><td></td><td>•</td><td></td></tr><tr><td>(Millero, 1979)</td><td>•</td><td>•</td><td>•</td><td>•</td><td></td><td></td><td></td><td>•</td><td></td></tr><tr><td>(Mojica Prieto and Millero, 2002)</td><td></td><td>•</td><td>•</td><td></td><td></td><td></td><td></td><td>•</td><td></td></tr><tr><td>(Cai and Wang, 1998)</td><td></td><td>•</td><td>•</td><td></td><td></td><td></td><td></td><td>•</td><td></td></tr><tr><td>(Millero et al., 2006)</td><td></td><td>•</td><td>•</td><td>•</td><td></td><td></td><td></td><td>•</td><td></td></tr><tr><td>(Millero, 2010)</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td>•</td><td></td></tr><tr><td>K<sub>0</sub> (Weiss, 1974)</td><td>•</td><td>•</td><td>•</td><td>•</td><td>•</td><td>•</td><td>•</td><td>•</td><td>•</td></tr><tr><td>K<sub>B</sub> (Dickson, 1990b)</td><td>•</td><td>•</td><td>•</td><td>•</td><td>•</td><td>•</td><td>•</td><td>•</td><td>•</td></tr><tr><td>K<sub>F</sub> (Perez and Fraga, 1987)</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td>•</td><td>•</td></tr><tr><td>K<sub>F</sub> (Dickson and Riley, 1979)</td><td>•</td><td>•</td><td>•</td><td>•</td><td>•</td><td>•</td><td>•</td><td>•</td><td>•</td></tr><tr><td>K<sub>W</sub> (Millero, 1995)</td><td>•</td><td>•</td><td>•</td><td>•</td><td>•</td><td>•</td><td>•</td><td>•</td><td>•</td></tr><tr><td>K<sub>S</sub> (Dickson, 1990a)</td><td>•</td><td>•</td><td>•</td><td>•</td><td>•</td><td>•</td><td>•</td><td>•</td><td>•</td></tr><tr><td>K<sub>S</sub> (Khoo et al., 1977)</td><td>•</td><td>•</td><td>•</td><td>•</td><td>•</td><td>•</td><td>•</td><td>•</td><td>•</td></tr><tr><td>K<sub>1P</sub>, K<sub>2P</sub>, K<sub>3P</sub> (Millero, 1995)</td><td>•</td><td>•</td><td>•</td><td>•</td><td>•</td><td>•</td><td>•</td><td>•</td><td>•</td></tr><tr><td>K<sub>Si</sub> (Millero, 1995)</td><td>•</td><td>•</td><td>•</td><td>•</td><td>•</td><td>•</td><td>•</td><td>•</td><td>•</td></tr><tr><td>K<sub>A</sub> (Mucci, 1983)</td><td>•</td><td>•</td><td>•</td><td>•</td><td>•</td><td>•</td><td>•</td><td>•</td><td>•</td></tr><tr><td>K<sub>C</sub> (Mucci, 1983)</td><td>•</td><td>•</td><td>•</td><td>•</td><td>•</td><td>•</td><td>•</td><td>•</td><td>•</td></tr></tbody></table>

## 2.3 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$ (Lo-carnini 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-

<footnote>a Both variants: CO2SYS-Excel-Pierrot and CO2SYS-Excel-Pelletier.</footnote>

<footnote>b All packages convert pH and K's between scales, internally.</footnote>

<footnote>c Some packages have user-callable routines to make these conversions between scales.</footnote>

<footnote>a Both variants: CO2SYS-Excel-Pierrot and CO2SYS-Excel-Pelletier.</footnote>

<footnote>b Refit by Dickson and Millero (1987).</footnote>
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₂, C_T-pCO₂, 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. More precisely, all other carbonate system variables were first calculated with one package (seacarb) from A_T = 2300 µmol kg⁻¹ and pCO₂ = 400 µatm at global average surface conditions (T = 18 °C, S = 35, and P = 0 db). Then two surface data sets were produced for each pair (and each package) by varying T and S, individually, and recalculating all other carbonate system variables from the fixed input pair. For the first, T was varied from -2 to 50 °C, while for the second S was varied from 0 to 50. In both cases, pressure was held at 0 db. To assess how packages differ below the surface, we used the same approach, varying pressure between 0 and 10 000 db and maintaining S = 35. But pressure corrections are highly sensitive to temperature (Sect. 2.7), so for each package we made two data sets: (1) holding T = 2 °C (typical of the deep open ocean) and (2) holding T = 13 °C (typical of deep waters of the Mediterranean Sea).

## 2.4 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₂ in seawater K₀ is from Weiss (1974). The equilibrium constants K₁ and K₂ 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 and Riley, 1966), fluoride (Riley, 1965), and Ca²⁺ (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²⁺ and CO₃²⁻ 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.

## 2.5 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₁ and K₂ 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).

## 2.6 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₀, K₁, K₂, K_B, K_W, K_A, and K_C; from ODV output, we computed its K₀, K₁, K₂, K_A, and K_C.

## 2.7 Pressure corrections

Until recently, no public package accounted for pressure effects on K₀ (needed to convert CO₂ to fCO₂) and the corresponding fugacity coefficient C_f (needed to convert fCO₂ to pCO₂) 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 fCO₂ and pCO₂ 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 fCO₂ and pCO₂ 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 $\theta$ 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(K_i^P/K_i^0) = -(\Delta V_i/RT_k)P + (0.5\Delta\kappa_i/RT_k)P^2, \quad (1) $$

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, $\Delta V_i$ is the partial molal volume, and $\Delta \kappa_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:

$$ \Delta V_i = a_0 + a_1 T_c + a_2 T_c^2, \quad (2) $$

$$ \Delta \kappa_i = b_0 + b_1 T_c + b_2 T_c^2, \quad (3) $$

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.

## 2.8 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.

## 3 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.

### 3.1 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 $fCO_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 $fCO_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 $fCO_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_3^{2-}$, 0.004 for $\Omega_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
<table><caption>**Table 7.** Coefficients used in Eqs. (2) and (3) to correct for effect of pressure on equilibrium constants.</caption><thead><tr><th>K</th><th>a<sub>0</sub></th><th>a<sub>1</sub></th><th>a<sub>2</sub></th><th>b<sub>0</sub></th><th>b<sub>1</sub></th><th>b<sub>2</sub></th></tr></thead><tbody><tr><td>K<sub>1</sub></td><td>-25.50</td><td>0.1271</td><td>0</td><td>-0.00308</td><td>-0.877 × 10<sup>-4</sup></td><td>0</td></tr><tr><td>K<sub>2</sub></td><td>-15.82</td><td>-0.0219</td><td>0</td><td>0.00113</td><td>-1.475 × 10<sup>-4</sup></td><td>0</td></tr><tr><td>K<sub>B</sub></td><td>-29.48</td><td>0.1622</td><td>-0.002608</td><td>-0.00284</td><td>0</td><td>0</td></tr><tr><td>K<sub>W</sub></td><td>-20.02</td><td>0.1119</td><td>-0.001409</td><td>-0.00513</td><td>-0.794 × 10<sup>-4</sup></td><td>0</td></tr><tr><td>K<sub>S</sub></td><td>-18.03</td><td>0.0466</td><td>-0.000316</td><td>-0.00453</td><td>-0.900 × 10<sup>-4</sup></td><td>0</td></tr><tr><td>K<sub>F</sub></td><td>-9.78</td><td>-0.0090</td><td>-0.000942</td><td>-0.00391</td><td>-0.540 × 10<sup>-4</sup></td><td>0</td></tr><tr><td>K<sub>C</sub></td><td>-48.76</td><td>0.5304</td><td>0</td><td>-0.01176</td><td>-3.692 × 10<sup>-4</sup></td><td>0</td></tr><tr><td>K<sub>A</sub></td><td>-45.96</td><td>0.5304</td><td>0</td><td>-0.01176</td><td>-3.692 × 10<sup>-4</sup></td><td>0</td></tr><tr><td>K<sub>1P</sub></td><td>-14.51</td><td>0.1211</td><td>-0.000321</td><td>-0.00267</td><td>-0.427 × 10<sup>-4</sup></td><td>0</td></tr><tr><td>K<sub>2P</sub></td><td>-23.12</td><td>0.1758</td><td>-0.002647</td><td>-0.00515</td><td>-0.900 × 10<sup>-4</sup></td><td>0</td></tr><tr><td>K<sub>3P</sub></td><td>-26.57</td><td>0.2020</td><td>-0.003042</td><td>-0.00408</td><td>-0.714 × 10<sup>-4</sup></td><td>0</td></tr><tr><td>K<sub>HS</sub></td><td>-14.80</td><td>0.0020</td><td>-0.000400</td><td>0.00289</td><td>-0.540 × 10<sup>-4</sup></td><td>0</td></tr><tr><td>K<sub>NH<sub>4</sub></sub></td><td>-26.43</td><td>0.0889</td><td>-0.000905</td><td>-0.00503</td><td>-0.814 × 10<sup>-4</sup></td><td>0</td></tr><tr><td>K<sub>Si</sub></td><td>-29.48</td><td>0.1622</td><td>-0.002608</td><td>-0.00284</td><td>0</td><td>0</td></tr></tbody></table>

$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.

<table><caption>**Table 8.** Desired measurement and numerical uncertainties.</caption><thead><tr><th rowspan="2">Variable</th><th colspan="2">Uncertainties</th><th rowspan="2">Units</th></tr><tr><th>Measurement</th><th>Numerical</th></tr></thead><tbody><tr><td>A<sub>T</sub></td><td>1</td><td>0.1</td><td>µmol kg<sup>-1</sup></td></tr><tr><td>C<sub>T</sub></td><td>1</td><td>0.1</td><td>µmol kg<sup>-1</sup></td></tr><tr><td>pCO<sub>2</sub></td><td>1</td><td>0.1</td><td>µatm</td></tr><tr><td>CO<sub>3</sub><sup>2-</sup></td><td>1</td><td>0.1</td><td>µmol kg<sup>-1</sup></td></tr><tr><td>pH</td><td>0.003</td><td>0.0003</td><td></td></tr><tr><td>pK<sub>0</sub></td><td>0.002</td><td>0.0002</td><td></td></tr><tr><td>pK<sub>1</sub></td><td>0.01</td><td>0.001</td><td></td></tr><tr><td>pK<sub>2</sub></td><td>0.02</td><td>0.002</td><td></td></tr><tr><td>pK<sub>i</sub> (other)</td><td>0.01</td><td>0.001</td><td></td></tr></tbody></table>

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.4 revised 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. Discrepancies

relative to the same reference were once larger, e.g., more than 10 times as much for $pCO_2$, pH, and $CO_3^{2-}$ (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_3^{2-}$ 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_3^{2-}$, 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 $fCO_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_3^{2-}$ 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.
<figure><img src="image_6.png" /><figcaption>**Figure 2.** Global zonal-mean surface values for variables computed from gridded data products for $A_T$ and $C_T$ from GLODAP (Key et al., 2004) combined with $T$, $S$, and nutrients from the 2009 World Ocean Atlas (WOA2009) (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 $\Omega_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.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**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.</figcaption></figure>
<table><caption>**Table 9.** Oceanic fCO₂ᵃ from three approachesᵇ in two packages.</caption><thead><tr><th rowspan="2">Depth (m)</th><th colspan="2">Common</th><th colspan="2">Potential</th><th colspan="2">In situ</th></tr><tr><th>mocsy</th><th>seacarb</th><th>mocsy</th><th>seacarb</th><th>mocsy</th><th>seacarb</th></tr></thead><tbody><tr><td>0</td><td>331.99</td><td>331.97</td><td>331.99</td><td>331.97</td><td>331.99</td><td>331.97</td></tr><tr><td>10</td><td>332.17</td><td>332.15</td><td>332.16</td><td>332.13</td><td>332.62</td><td>332.59</td></tr><tr><td>50</td><td>348.50</td><td>348.48</td><td>348.42</td><td>348.40</td><td>350.85</td><td>350.83</td></tr><tr><td>100</td><td>399.82</td><td>399.79</td><td>399.65</td><td>399.62</td><td>405.28</td><td>405.25</td></tr><tr><td>500</td><td>628.72</td><td>628.68</td><td>627.67</td><td>627.62</td><td>674.20</td><td>674.15</td></tr><tr><td>1000</td><td>671.15</td><td>671.10</td><td>669.24</td><td>669.19</td><td>773.06</td><td>773.01</td></tr><tr><td>2000</td><td>551.12</td><td>551.08</td><td>548.18</td><td>548.14</td><td>733.22</td><td>733.17</td></tr><tr><td>5000</td><td>438.05</td><td>438.02</td><td>430.74</td><td>430.71</td><td>900.90</td><td>900.83</td></tr></tbody></table>

ᵃ Area-weighted global means computed from the GLODAP-WOA2009 gridded data set.

ᵇ Following Weiss (1974), $f\text{CO}_2 = [\text{CO}_2]/(K_0 \exp[(1-P)\bar{v}\text{CO}_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₀.

<table><caption>**Table 10.** Oceanic pCO₂ᵃ from three approachesᵇ in two packages.</caption><thead><tr><th rowspan="2">Depth (m)</th><th colspan="2">Common</th><th colspan="2">Potential</th><th colspan="2">In situ</th></tr><tr><th>mocsy</th><th>seacarb</th><th>mocsy</th><th>seacarb</th><th>mocsy</th><th>seacarb</th></tr></thead><tbody><tr><td>0</td><td>333.15</td><td>333.12</td><td>333.15</td><td>333.12</td><td>333.1</td><td>333.1</td></tr><tr><td>10</td><td>333.33</td><td>333.30</td><td>333.31</td><td>333.28</td><td>334.9</td><td>334.9</td></tr><tr><td>50</td><td>349.73</td><td>349.70</td><td>349.65</td><td>349.62</td><td>358.3</td><td>358.3</td></tr><tr><td>100</td><td>401.27</td><td>401.24</td><td>401.09</td><td>401.06</td><td>421.6</td><td>421.6</td></tr><tr><td>500</td><td>631.24</td><td>631.20</td><td>630.19</td><td>630.14</td><td>826.0</td><td>825.9</td></tr><tr><td>1000</td><td>673.95</td><td>673.90</td><td>672.03</td><td>671.98</td><td>1175.3</td><td>1175.2</td></tr><tr><td>2000</td><td>553.48</td><td>553.44</td><td>550.53</td><td>550.49</td><td>1729.7</td><td>1729.6</td></tr><tr><td>5000</td><td>439.95</td><td>439.92</td><td>432.62</td><td>432.58</td><td>7976.1</td><td>7975.6</td></tr></tbody></table>

ᵃ Area-weighted global means from same gridded input data as in Table 9.

ᵇ Following Weiss (1974), $p\text{CO}_2 = f\text{CO}_2/C_f$, where $C_f = \exp[(B + 2x_2^2\delta_{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.

## 3.2 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-p\text{CO}_2$ and $C_T-p\text{CO}_2$, excluded the csys package, which does offer pCO₂ as an input variable but only when paired with pH.

### 3.2.1 $A_T-C_T$

With the $A_T-C_T$ pair, packages agree within 0.2 µatm in pCO₂, 0.05 µmol kg⁻¹ in CO₂⁻, 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₂ from ODV, but those re-

main quite small (less than twice our arbitrary numerical cut-off 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°C but only negligible discrepancies at 2°C. At 5000 db, its discrepancies at 13 °C reach -2 µatm for pCO₂, +1 µmol kg⁻¹ for CO₂⁻, and +0.002 for pH. Discrepancies in other packages remain negligible even at 13°C.
<figure><img src="image_8.png" /></figure>

**Figure 4.** Global zonal-mean surface values (top) and global-mean vertical profiles (bottom) from outdated versions of packages for *pCO₂* (left), *CO₃²⁻* (middle), and pH (right) as computed from GLODAP *A*<sub>T</sub> and *C*<sub>T</sub> as in Figs. 3 and 2. The four older versions include CO₂calc (v1.0.4), csys (from 3 February 2010), seacarb (v2.3.3), and mocsy (non public predecessor from Orr et al., 2005). As before, results are shown after subtracting off corresponding results from the same CO₂SYS-MATLAB reference.

### 3.2.2 A<sub>T</sub>-pH

With the *A*<sub>T</sub>-pH input pair (Fig. 6), surface discrepancies between packages remain negligible for all variables. All packages agree within 0.02 µatm in *pCO₂*, 0.02 µmol kg<sup>-1</sup> in *CO₃²⁻*, and 0.08 µmol kg<sup>-1</sup> in *C*<sub>T</sub> 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.

### 3.2.3 A<sub>T</sub>-pCO₂

With the *A*<sub>T</sub>-*pCO₂* input pair (Fig. 7), surface discrepancies are always negligible. The five packages differ by less than 0.05 µmol kg<sup>-1</sup> in *C*<sub>T</sub> and 0.015 µmol kg<sup>-1</sup> in *CO₃²⁻*. Likewise for pH, packages generally agree within 0.0001; only CO₂calc exhibits larger variability (within ±0.0004), but those variations are randomly distributed with a mean near zero, a consequence of CO₂calc'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 +1 µmol kg<sup>-1</sup> in *C*<sub>T</sub>, +0.1 µmol kg<sup>-1</sup> in *CO₃²⁻*, and +0.0002 in pH at 5000 db.

### 3.2.4 C<sub>T</sub>-pH

With the *C*<sub>T</sub>-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₂*, 0.007 µmol kg<sup>-1</sup> for *CO₃²⁻*, and 0.1 µmol kg<sup>-1</sup> for *A*<sub>T</sub> at surface pressure. With *C*<sub>T</sub>-pH, unlike with previously analyzed pairs, the swco2 package's pressure corrections do not induce substantial discrepancies in computed subsurface *pCO₂* and *CO₃²⁻*, even at 13 °C. Yet swco2 does have significant discrepancies in computed subsurface *A*<sub>T</sub> (e.g., 1 µmol kg<sup>-1</sup> at 4000 db); conversely, with the pressure correction at *T* = 2 °C, swco2's *A*<sub>T</sub> 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.

### 3.2.5 C<sub>T</sub>-pCO₂

With the *C*<sub>T</sub>-*pCO₂* input pair (Fig. 9), all of the five packages have negligible surface discrepancies for computed *A*<sub>T</sub> (≤ 0.1 µmol kg<sup>-1</sup>), *CO₃²⁻* (≤ 0.01 µmol kg<sup>-1</sup>), and pH (≤ 0.003 units). Out of the five packages offering both the *C*<sub>T</sub>-*pCO₂* and the *C*<sub>T</sub>-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*<sub>T</sub> grow linearly with depth, reaching 1 µmol kg<sup>-1</sup> at 4000 db, similar to those seen with the *C*<sub>T</sub>-pH input pair
<figure><img src="image_9.png" /><figcaption>**Figure 5.** Variables computed from $A_T$ and $C_T$ for each package minus corresponding results from CO2SYS-MATLAB. The computed $p\text{CO}_2$ (top), $\text{CO}_3^{2-}$ (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.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Figure 6.** Variables computed from $A_T$ and pH with each package minus corresponding results from CO2SYS-MATLAB. Shown are computed $p\text{CO}_2$ (top), $\text{CO}_3^{2-}$ (middle), and $C_T$ (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).</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>**Figure 7.** Variables computed from $A_T$ and $pCO_2$ with each package minus corresponding results from CO2SYS-MATLAB. Shown are computed $C_T$ (top), $CO_3^{2-}$ (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. Packages not included are mocsy, which allows only the $A_T-C_T$ pair, and csys, which does not allow $pCO_2$ as an input variable.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**Figure 8.** Variables computed from $C_T$ and pH with each package minus corresponding results from CO2SYS-MATLAB. Shown are computed $pCO_2$ (top), $CO_3^{2-}$ (middle), and $A_T$ (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 software package.</figcaption></figure>
(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.

## 3.3 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).

### 3.3.1 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⁻¹ in $CO_3^{2-}$, 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_3^{2-}$ is reduced by about 2 µmol kg⁻¹ in the same waters (Orr and Epitalon, 2015).

### 3.3.2 Total boron

Relative to the standard formulation for total boron (Oppströ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⁻¹ for $CO_3^{2-}$, 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_3^{2-}$ by 2.1 µmol kg⁻¹, 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⁻¹ for $CO_3^{2-}$, respectively). Conversely, with the $C_T$-pH and $C_T-pCO_2$ pairs, changes are negligible for all computed carbonate system variables except total alkalinity.

### 3.3.3 $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⁻¹ in $CO_3^{2-}$, and +0.006 units in pH. Discrepancies are also found for seacarb, reaching up to -20 µatm in $pCO_2$, ±0.2 µmol kg⁻¹ in $CO_3^{2-}$, 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⁻¹ in $CO_3^{2-}$, and -0.001 units in pH, changes that are notably less than its surface discrepancies. Although
<figure><img src="image_13.png" /><figcaption>**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_3^{2-}$ (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.</figcaption></figure>

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.

## 3.4 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.

## 4 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.

### 4.1 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 $\partial y/y : \partial 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.185$ and $A_T = 2300 \mu mol kg^{-1}$). Then, we increased each input variable by 1% ($\partial x_i/x_i = 0.01$), individually, and recalculated out-
<figure><img src="image_14.png" /><figcaption>**Figure 10.** 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_3^{2-}$ (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.</figcaption></figure>

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 $\partial y/y$.

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_3^{2-}$ as well as $\Omega_A$ and $\Omega_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_3^{2-}$ 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 $fCO_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_3^{2-}$, $\Omega_A$, and $\Omega_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 $\Omega_A$ and $\Omega_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

$$ \frac{u(y)}{y} = \sqrt{\sum_{i=1}^{n} \left( \frac{\partial y/y}{\partial x_i/x_i} \right)^2 \left( \frac{u(x_i)}{x_i} \right)^2}, \quad (4) $$
<figure><img src="image_15.png" /><figcaption>**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_3^{2-}$ (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.</figcaption></figure>

<table><caption>**Table 11.** Ratio of relative change<sup>a,b</sup> between output and input variables<sup>c,d</sup> $(\partial y/\partial x) / (\partial x_i/x_i)$.</caption><thead><tr><th rowspan="2">Output</th><th colspan="10">Input</th></tr><tr><th>A<sub>T</sub></th><th>C<sub>T</sub></th><th>K<sub>0</sub></th><th>K<sub>1</sub></th><th>K<sub>2</sub></th><th>K<sub>B</sub></th><th>K<sub>W</sub></th><th>K<sub>S</sub></th><th>K<sub>F</sub></th><th>K<sub>A</sub></th><th>K<sub>C</sub></th></tr></thead><tbody><tr><td>CO<sub>2</sub><sup>*</sup></td><td>-9.5</td><td>11.7</td><td></td><td>-0.92</td><td>0.63</td><td>0.28</td><td>0.02</td><td></td><td></td><td></td><td></td></tr><tr><td>HCO<sub>3</sub><sup>-</sup></td><td>-0.7</td><td>1.7</td><td></td><td>0.01</td><td>-0.03</td><td>0.02</td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>CO<sub>3</sub><sup>2-</sup></td><td>8.8</td><td>-7.4</td><td></td><td>-0.05</td><td>0.30</td><td>-0.24</td><td>-0.01</td><td></td><td></td><td></td><td></td></tr><tr><td>H<sup>+</sup></td><td>-8.8</td><td>9.8</td><td></td><td>0.06</td><td>0.66</td><td>0.26</td><td>0.01</td><td></td><td></td><td></td><td></td></tr><tr><td>pCO<sub>2</sub></td><td>-9.5</td><td>11.7</td><td>-0.99</td><td>-0.92</td><td>0.63</td><td>0.28</td><td>0.02</td><td></td><td></td><td></td><td></td></tr><tr><td>fCO<sub>2</sub></td><td>-9.5</td><td>11.7</td><td>-0.99</td><td>-0.92</td><td>0.63</td><td>0.28</td><td>0.02</td><td></td><td></td><td></td><td></td></tr><tr><td>Ω<sub>A</sub></td><td>8.8</td><td>-7.4</td><td></td><td>-0.05</td><td>0.30</td><td>-0.24</td><td>-0.01</td><td></td><td>-0.99</td><td></td><td></td></tr><tr><td>Ω<sub>C</sub></td><td>8.8</td><td>-7.4</td><td></td><td>-0.05</td><td>0.30</td><td>-0.24</td><td>-0.01</td><td></td><td></td><td>-0.99</td><td></td></tr></tbody></table>

<footnote><sup>a</sup> Percent change in output ($\partial y/\partial x$) computed from a 1% change in input ($\partial x/x$).</footnote>

<footnote><sup>b</sup> Missing values indicate changes of less than 0.001 %.</footnote>

<footnote><sup>c</sup> Ratios were computed under the standard conditions described in Sect. 4.1.</footnote>

<footnote><sup>d</sup> Ratios are sensitive to solution composition.</footnote>
<figure><img src="image_16.png" /><figcaption>**Figure 12.** 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 $p\text{CO}_2$ (top), $\text{CO}_3^{2-}$ (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).</figcaption></figure>

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 ($\partial y/y$) / ($\partial 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.

## 4.2 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.

### 4.2.1 $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.

### 4.2.2 Alternative $K_1$ and $K_2$ for low salinities

Section 3.3.3 detailed 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_1$ and $K_2$ relative to those computed by
<table><caption>**Table 12.** Check values vs. computed equilibrium constants (pKs) at T = 25°C, S = 35, and P = 0.</caption><thead><tr><th></th><th>pK<sub>0</sub></th><th>pK<sub>1</sub></th><th>pK<sub>2</sub></th><th>pK<sub>B</sub></th><th>pK<sub>W</sub></th><th>pK<sub>S</sub><sup>a</sup></th><th>pK<sub>F</sub><sup>a</sup></th></tr></thead><tbody><tr><td>Check value</td><td>1.5468</td><td>5.8472</td><td>8.966</td><td>8.5975</td><td>13.217</td><td>0.999</td><td>2.627</td></tr><tr><td>swco2</td><td>1.54681</td><td>5.84715</td><td>8.96595</td><td>8.59747</td><td>13.2204</td><td>0.9987</td><td>2.6261</td></tr><tr><td>mocsy</td><td>1.54681</td><td>5.84715</td><td>8.96595</td><td>8.59747</td><td>13.2204</td><td>0.9987</td><td>2.6261</td></tr><tr><td>CO2SYS-MATLAB</td><td>1.54681</td><td>5.84715</td><td>8.96595</td><td>8.59747</td><td>13.2204</td><td>0.9987</td><td>2.6261</td></tr><tr><td>CO2SYS-Excel-Pelletier</td><td>1.54681</td><td>5.84715</td><td>8.96595</td><td>8.59747</td><td>13.2204</td><td>0.9987</td><td>2.6261</td></tr><tr><td>csys</td><td>1.54681</td><td>5.84715</td><td>8.96595</td><td>8.59747</td><td>13.2204</td><td>0.9987</td><td>2.6261</td></tr><tr><td>seacarb</td><td>1.54681</td><td>5.84715</td><td>8.96595</td><td>8.59747</td><td>13.2204</td><td>0.9987</td><td>2.6261</td></tr><tr><td>CO2calc</td><td>1.5468</td><td>5.847<sup>b</sup></td><td>8.966<sup>b</sup></td><td>8.597</td><td>13.220</td><td></td><td></td></tr><tr><td>ODV</td><td>1.54680</td><td>5.84715</td><td>8.96595</td><td></td><td></td><td></td><td></td></tr></tbody></table>

<table><caption>**Table 12.** Continued.</caption><thead><tr><th></th><th>pK<sub>C</sub></th><th>pK<sub>A</sub></th><th>pK<sub>1P</sub></th><th>pK<sub>2P</sub></th><th>pK<sub>3P</sub></th><th>pK<sub>Si</sub></th></tr></thead><tbody><tr><td>Check value</td><td></td><td></td><td>1.61</td><td>5.962</td><td>8.79</td><td>9.39</td></tr><tr><td>swco2</td><td>6.3693</td><td>6.1883</td><td>1.6215</td><td>5.9714</td><td>8.799</td><td>9.393</td></tr><tr><td>mocsy</td><td>6.3693</td><td>6.1883</td><td>1.6150</td><td>5.9649</td><td>8.792</td><td>9.387</td></tr><tr><td>CO2SYS-MATLAB</td><td>6.3693</td><td>6.1883</td><td>1.6150</td><td>5.9649</td><td>8.793</td><td>9.387</td></tr><tr><td>CO2SYS-Excel-Pelletier</td><td>6.3693</td><td>6.1883</td><td>1.6150</td><td>5.9649</td><td>8.792</td><td>9.387</td></tr><tr><td>csys</td><td>6.3693</td><td>6.1883</td><td>1.6150</td><td>5.9649</td><td>8.792</td><td></td></tr><tr><td>seacarb</td><td>6.3693</td><td>6.1883</td><td>1.6150</td><td>5.9649</td><td>8.792</td><td>9.387</td></tr><tr><td>CO2calc</td><td>6.369</td><td>6.188</td><td></td><td></td><td></td><td></td></tr><tr><td>ODV</td><td>6.3693</td><td>6.1883</td><td></td><td></td><td></td><td></td></tr></tbody></table>

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₁ and 590 measurements of K₂ on the seawater pH scale to basic equations of the form

$$pK_i - pK_i^0 = A_i + B_i/T_k + C_i \ln(T_k), \quad (5)$$

where pKᵢ⁰ was from his previous fit of the same form for pure water (Millero et al., 2006); Tₖ is the absolute temperature; and Aᵢ, Bᵢ, and Cᵢ are functions of salinity:

$$ \begin{align} A_i &= \alpha_0 S^{0.5} + \alpha_1 S + \alpha_2 S^2, \notag \\ B_i &= \alpha_3 S^{0.5} + \alpha_4 S, \tag{6} \\ C_i &= \alpha_5 S^{0.5}. \notag \end{align} $$

Millero (2010) provides this set of the six αᵢ coefficients for each of K₁ and K₂ on the original seawater pH scale but also on the free and total scales, i.e., by making 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₁ and K₂; 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₁ and for K₂, e.g., once the results from the seawater set are converted back to the total scale. Yet they do not. Indeed when K₁ and K₂ 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ₔ. That is, Millero converted pK₁ and pK₂ on the measured

<footnote>a Free scale (all other pKs are on the total scale).</footnote>

<footnote>b All CO2calc constants shown here are computed from other output variables; the pK₁ and pK₂ calculated directly by CO2calc, provided by its developers, agree with the reference to beyond six decimal places (L. Robbins, personal communication, 2014).</footnote>
<figure><img src="image_17.png" /><figcaption>**Figure 13.** 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$ from Weiss (1974) and $K_1$ and $K_2$ from Lueker 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 (Table 12). Line signatures are as in Fig. 5.</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>**Figure 14.** The $pK_1$ (top) and $pK_2$ (bottom) computed in each of four packages minus corresponding values for CO2SYS-MATLAB, as in Fig. 13 except that formulations are from Millero (2010), not Lueker et al. (2000). The filled black circles indicate the check values (Table 12). 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.</figcaption></figure>
<figure><img src="image_19.png" /><figcaption>**Figure 15.** Relative differences between four formulations for $pK_1$ (top) as well as $pK_2$ (bottom), i.e., between values computed from four different sets of coefficients from Millero (2010) over ranges of $T$ (left) and $S$ (right). Two sets of coefficients are from Millero's (2010) Tables 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).</figcaption></figure>

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 $\alpha_5$ coefficient 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.

<table><caption>**Table 13.** Coefficients from Millero (2010) for formulations of $K_1$ and $K_2$ (seawater scale).</caption><thead><tr><th></th><th>K<sub>1</sub></th><th>K<sub>2</sub></th></tr></thead><tbody><tr><td>α<sub>0</sub></td><td>13.4038</td><td>21.3728</td></tr><tr><td>α<sub>1</sub></td><td>0.03206</td><td>0.1218</td></tr><tr><td>α<sub>2</sub></td><td>-5.242 × 10<sup>-5</sup></td><td>-3.688 × 10<sup>-4</sup></td></tr><tr><td>α<sub>3</sub></td><td>-530.659</td><td>-788.289</td></tr><tr><td>α<sub>4</sub></td><td>-5.8210</td><td>-19.189</td></tr><tr><td>α<sub>5</sub></td><td>-2.0664</td><td>-3.374*</td></tr></tbody></table>

<figure><img src="image_20.png" /><figcaption>**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.</figcaption></figure>

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

<footnote>* Value is -3.3738 in Millero's spreadsheet (F. J. Millero, personal communication, 2013).</footnote>
<figure><img src="image_21.png" /><figcaption>**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 (Table 12).</figcaption></figure>

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 $\alpha_1$ for $K_1$ for the total scale would need to be extended from 4 to 5 decimal places, while the $\alpha_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.

### 4.2.3 $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. This same error was only identified and corrected in other packages (e.g., in csys and seacarb) in 2010, whereas the current version of swco2 (v2) dates from 2007. Because the $a_2$ coefficient is multiplied by the square of $T_c$ (Eq. 2), the swco2 discrepancies associated with the pressure correction of $K_B$ are about
<figure><img src="image_22.png" /><figcaption>**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).</figcaption></figure>

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.

### 4.2.4 $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 $\Delta 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 $\Delta 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.

### 4.2.5 $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 $\Delta pK_A \sim 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_3^{2-}$). Pressure-correction discrepancies are negligible in all packages but one, swco2. At 5000 db, the discrepancy for 13 °C water reaches $\Delta pK_A = 0.015$, thereby biasing $K_A$ to be 3.4% too low and $\Omega_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.
<figure><img src="image_23.png" /><figcaption>**Figure 19.** The $pK_{1P}$ (top), $pK_{2P}$ (middle), and $pK_{3P}$ (bottom) computed from each package minus corresponding values from CO2SYS-MATLAB. Formulations are from Millero (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 (Table 12).</figcaption></figure>

### 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 KSi are lower when carried out at 2°C rather than at 13°C. The total discrepancies in ΔpKSi 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.

## 5 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 pCO2 and CO3²⁻ 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 pCO2 and CO3²⁻ are computed from the A_T – C_T input pair. Those discrepancies appear to derive from a sign error in a pressure-correction coefficient for KB. Similar sign errors in pressure-correction 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₁ and K₂ (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₁ and K₂ 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 pCO2. One package differs because it uses the set of published coefficients to compute the salinity dependence of K₁ and K₂ 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 coefficients 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₁ and K₂ 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₁ and K₂ from Millero (2010) until discrepancies are resolved. For now, users are also advised to avoid the new total boron-to-salinity ratio (Lee et al., 2010) and favor the “best-practices” ratio (Uppström, 1974), which was used to compute K₁ and K₂ 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 pCO2, 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.**

*Acknowledgements.* This comparison was conceived by the SOLAS-IMBER Ocean Acidification Working Group (SIOA). Design of the comparison, data production, and preliminary analysis was originally supported by the International Ocean Carbon Coordination Project (IOCCP) with funds from the U.S. National Science Foundation (grant OCE-1243377) to the Scientific Committee on Oceanic Research (SCOR). Follow-up support came from the Ocean Acidification International Coordination Centre (OA-ICC), a project of the Peaceful Uses Initiative of the International Atomic Energy Agency (IAEA). Final analysis and the publication of this article were supported by the French ANR project MACROES (MACRoscope for Oceanic Earth System ANR-09-CEP-003) and two EU FP7 projects, CARBOCHANGE (grant 264879) and MedSeA (grant 265103). We thank F. J. Millero for providing the spreadsheet associated with his 2010 publication, with which we were able decipher causes of some differences in implementations of his formulations for K₁ and K₂. We are indebted to A. Dickson and J. Ólafsson, who provided reviews that helped to improve this study. We also appreciate the many helpful comments from E. Lewis, D. Wallace, L. Robbins, R. Zeebe, and D. Wolf-Gladrow.

Edited by: L. Cotrim da Cunha
## References

Antonov, J. I., Seidov, D., Boyer, T. P., Locarnini, R. A., Mishonov, A. V., Garcia, H. E., Baranova, O. K., Zweng, M. M., and Johnson, D. R.: World Ocean Atlas 2009, Volume 2: Salinity, edited by: Levitus, S., Atlas NESDIS 69, NOAA, US Government Printing Office, Washington DC, 184 pp., 2010.

Byrne, R. H. and Yao, W.: Procedures for measurement of carbonate ion concentrations in seawater by direct spectrophotometric observations of Pb (II) complexation, Mar. Chem., 112, 128–135, 2008.

Cai, W.-J. and Wang, Y.: The chemistry, fluxes, and sources of carbon dioxide in the estuarine waters of the Satilla and Altamaha Rivers, Georgia, Limnol. Oceanogr., 43, 657–668, 1998.

Culberson, C. and Pytkowicz, R. M.: Effect of pressure on carbonic acid, boric acid and the pH in seawater, Limnol. Oceanogr., 13, 403–417, 1968.

Dickson, A. G.: Standard potential of the reaction: $AgCl(s) + 1/2H_2(g) = Ag(s) + HCl(aq)$, and the standard acidity constant of the ion $HSO_4^-$ in synthetic seawater from 273.15 to 318.15 K, Chemical Thermodynamics, 22, 113–127, 1990a.

Dickson, A. G.: Thermodynamics of the dissociation of boric acid in synthetic seawater from 273.15 to 318.15 K, Deep-Sea Res., 37, 755–766, 1990b.

Dickson, A. G.: The carbon dioxide system in seawater: equilibrium chemistry and measurements, in: Guide to Best Practices for Ocean Acidification Research and Data Reporting, edited by: Riebesell, U., Fabry, V. J., Hansson, L., and Gattuso, J.-P., Publications Office of the European Union, 17–40, 2010.

Dickson, A. G. and Goyet, C.: Handbook of methods for the analysis of the various parameters of the carbon dioxide system in seawater; version 2, ORNL/CDIAC-74, ORNL, 1994.

Dickson, A. G. and Millero, F.: A comparison of the equilibrium constants for the dissociation of carbonic acid in seawater media, Deep-Sea Res., 34, 1733–1743, 1987.

Dickson, A. G. and Riley, J. P.: The effect of analytical error on the evaluation of the components of the aquatic carbon-dioxide system, Mar. Chem., 6, 77–85, 1978.

Dickson, A. G. and Riley, J. P.: The estimation of acid dissociation constants in seawater media from potentiometric titrations with strong base, I. The ionic product of water – $K_w$, Mar. Chem., 7, 89–99, 1979.

Dickson, A. G., Sabine, C. L., and Christian, J. R.: Guide to best practices for ocean $CO_2$ measurements, PICES Special Publication 3, 191 pp., 2007.

Easley, R. A., Patsavas, M. C., Byrne, R. H., Liu, X., Feely, R. A., and Mathis, J. T.: Spectrophotometric Measurement of Calcium Carbonate Saturation States in Seawater, Environ. Sci. Technol., 47, 1468–1477, 2013.

Frankignoulle, M.: A complete set of buffer factors for acid/base $CO_2$ system in seawater, J. Marine Syst., 5, 111–118, 1994.

Garcia, H. E., Locarnini, R. A., Boyer, T. P., Antonov, J. I., Zweng, M. M., Baranova, O. K., and Johnson, D. R.: World Ocean Atlas 2009, Volume 4: Nutrients (Phosphate, Nitrate, Silicate), edited by: Levitus, S., Atlas NESDIS 71, NOAA, US Government Printing Office, Washington DC, 398 pp., 2010.

Gattuso, J.-P. and Hansson, L. (Eds.): Ocean Acidification, Oxford Univ. Press, Oxford, UK, 408 pp., 2011.

Gattuso, J.-P. and Lavigne, H.: Technical Note: Approaches and software tools to investigate the impact of ocean acidification, Biogeosciences, 6, 2121–2133, doi:10.5194/bg-6-2121-2009, 2009.

Gattuso, J.-P., Epitalon, J.-M. and Lavigne, H.: seacarb: seawater carbonate chemistry with R package version 3.0.6, The Comprehensive R Archive Network, http://CRAN.R-project.org/package/seacarb (last access: 24 February 2015), 2015.

Goyet, C. and Poisson, A.: New determination of carbonic acid dissociation constants in seawater as a function of temperature and salinity, Deep-Sea Res., 36, 1635–1654, 1989.

Hansson, I.: A new set of acidity constants for carbonic acid and boric acid in sea water, Deep-Sea Res., 20, 461–478, 1973a.

Hansson, I.: The determination of dissociation constants of carbonic acid in synthetic sea water in the salinity range of 20–40‰ and temperature range of 5–30 °C, Acta Chemica Scandanavia, 27, 931–944, 1973b.

Hofmann, A. F., Soetaert, K., Middelburg, J. J., and Meysman, F. J.: AquaEnv: An Aquatic Acid-Base Modelling Environment in: R. Aquatic Geochemistry, 16, 507–546, 2010.

Hoppe, C. J. M., Langer, G., Rokitta, S. D., Wolf-Gladrow, D. A., and Rost, B.: Implications of observed inconsistencies in carbonate chemistry measurements for ocean acidification studies, Biogeosciences, 9, 2401–2405, doi:10.5194/bg-9-2401-2012, 2012.

Hunter, K. A.: XLCO2 – Seawater $CO_2$ Equilibrium Calculations Using Excel Version 2, University of Otago, New Zealand, available at: http://neon.otago.ac.nz/research/mfc/people/keith_hunter/software/swco2/ (last access: 24 February 2015), 2007.

Key, R. M., Sabine, C. L., Lee, K., Wanninkhof, R., Bullister, J., Feely, R. A., Millero, F. J., Mordy, C., and Peng, T.-H.: A Global Ocean carbon climatology: results from Global Data Analysis Project (GLODAP), Global Biogeochem. Cy., 18, GB4031, doi:10.1029/2004GB002247, 2004.

Khoo, K. H., Ramette, R. W., Culberson, C. H., and Bates, R. G.: Determination of hydrogen ion concentrations in seawater from 5 to 40 °C: standard potentials at salinities from 20 to 45‰, Anal. Chem., 49, 29–34, 1977.

Kroeker, K. J., Kordas, R. L., Crim, R., Hendriks, I. E., Ramajo, L., Singh, G. S., Duarte, C. M., and Gattuso, J.-P.: Impacts of ocean acidification on marine organisms: quantifying sensitivities and interaction with warming, Glob. Change Biol., 19, 156–159, 2013.

Lee, K., Kim, T.-W., Byrne, R. H., Millero, F. J., Feely, R. A., and Liu, Y.-M.: The universal ratio of boron to chlorinity for the North Pacific and North Atlantic oceans, Geochim. Cosmochim. Ac., 74, 1801–1811, doi:10.1016/j.gca.2009.12.027, 2010.

Lewis, E. and Wallace, D. W. R.: Program Developed for $CO_2$ System Calculations, ORNL/CDIAC-105, Carbon Dioxide Inf. Anal. Cent., Oak Ridge Natl. Lab., Oak Ridge, Tenn., 38 pp., 1998.

Locarnini, R., Mishonov, A., Antonov, J., Boyer, T., Garcia, H., Baranova, O., Zweng, M., and Johnson, D.: World Ocean Atlas 2009, Volume 1: Temperature, edited by: Levitus, S., Atlas NESDIS 68, NOAA, US Government Printing Office, Washington DC, 184 pp., 2010.

Lueker, T. J., Dickson, A. G., and Keeling, C. D.: Ocean $pCO_2$ calculated from dissolved inorganic carbon, alkalinity, and equations for $K_1$ and $K_2$: validation based on laboratory measurements of $CO_2$ in gas and seawater at equilibrium, Mar. Chem., 70, 105–119, 2000.
Martz, T. R., Jannasch, H. W., and Johnson, K. S.: Determination of carbonate ion concentration and inner sphere carbonate ion pairs in seawater by ultraviolet spectrophotometric titration, Mar. Chem., 115, 145–154, 2009.

Mehrbach, C., Culberson, C. H., Hawley, J. E., and Pytkowicz, R. M.: Measurement of the apparent dissociation constants of carbonic acid in seawater at atmospheric pressure, Limnol. Oceanogr., 18, 897–907, 1973.

Millero, F. J.: The thermodynamics of the carbonate system in seawater, Geochim. Cosmochim. Ac., 43, 1651–1661, 1979.

Millero, F. J.: Influence of pressure on chemical processes in the sea, in: Chemical Oceanography, edited by: Riley, J. P. and Chester, R., 1–88, Academic Press, New York, 1983.

Millero, F. J.: Thermodynamics of the carbon dioxide system in the oceans, Geochim. Cosmochim. Ac., 59, 661–677, 1995.

Millero, F. J.: The marine inorganic carbon cycle, Chem. Rev., 107, 308–341, 2007.

Millero, F. J.: Carbonate constants for estuarine waters, Mar. Freshwater Res., 61, 139–142, doi:10.1071/MF09254, 2010.

Millero, F. J., Graham, T. B., Huang, F., Bustos-Serrano, H., and Pierrot, D.: Dissociation constants of carbonic acid in seawater, Mar. Chem., 100, 80–94, doi:10.1016/j.marchem.2005.12.001, 2006.

Mojica Prieto, F. J. and Millero, F. J.: The values of $pK_1 + pK_2$ for the dissociation of carbonic acid in seawater, Geochim. Cosmochim. Ac., 66, 2529–2540, 2002.

Morris, A. and Riley, J.: The bromide/chlorinity and sulphate/chlorinity ratio in sea water, Deep-Sea Res., 13, 699–705, 1966.

Mosley, L. M., Peake, B. M., and Hunter, K. A.: Modelling of pH and inorganic carbon speciation in estuaries using the composition of the river and seawater end members, Environ. Modell. Softw., 25, 1658–1663, 2010.

Mucci, A.: The solubility of calcite and aragonite in seawater at various salinities, temperatures, and one atmosphere total pressure, Am. J. Sc., 283, 780–799, 1983.

Munhoven, G.: Mathematics of the total alkalinity-pH equation – pathway to robust and universal solution algorithms: the SolveSAPHE package v1.0.1, Geosci. Model Dev., 6, 1367–1388, doi:10.5194/gmd-6-1367-2013, 2013.

Orr, J. C.: Recent and future changes in ocean carbonate chemistry, in: Ocean Acidification, edited by: Gattuso, J.-P., and Hansson, L., chap. 3, 41–66, Oxford Univ. Press, 2011.

Orr, J. C. and Epitalon, J.-M.: Improved routines to model the ocean carbonate system: mocsy 1.0, Geosci. Model Dev. Discuss., 7, 2877–2902, doi:10.5194/gmdd-7-2877-2014, 2014.

Orr, J. C. and Epitalon, J.-M.: Improved routines to model the ocean carbonate system: mocsy 2.0, Geosci. Model Dev., 8, 485–499, doi:10.5194/gmd-8-485-2015, 2015.

Orr, J. C., Najjar, R., Sabine, C. L., and Joos, F.: Abiotic-HOWTO. Internal OCMIP Report, LSCE/CEA Saclay, Gif-sur-Yvette, France, 25 pp., 1999.

Orr, J. C., Fabry, V. J., Aumont, O., Bopp, L., Doney, S. C., Feely, R. A., Gnanadesikan, A., Gruber, N., Ishida, A., Joos, F., Key, R. M., Lindsay, K., Maier-Reimer, E., Matear, R., Monfray, P., Mouchet, A., Najjar, R. G., Plattner, G.-K., Rodgers, K. B., Sabine, C. L., Sarmiento, J. L., Schlitzer, R., Slater, R. D., Totterdell, I. J., Weirig, M.-F., Yamanaka, Y., and Yool, A.: Anthropogenic ocean acidification over the twenty-first century and its impact on calcifying organisms, Nature, 437, 681–686, doi:10.1038/nature04095, 2005.

Orr, J. C., Epitalon, J.-M., and Gattuso, J.-P.: Comparison of seven packages that compute ocean carbonate chemistry, Biogeosciences Discuss., 11, 5327–5397, doi:10.5194/bgd-11-5327-2014, 2014.

Pelletier, G., Lewis, E., and Wallace, D.: CO2SY S.XLS: A calculator for the CO₂ system in seawater for Microsoft Excel/VBA, Wash. State Dept. of Ecology/Brookhaven Nat. Lab., Olympia, WA/Upton, NY, USA, 2007.

Perez, F. F. and Fraga, F.: Association constant of fluoride and hydrogen ions in seawater, Mar. Chem., 21, 161–168, 1987.

Pierrot, D., Lewis, E., and Wallace, D. W. R.: MS Excel Program Developed for CO₂ System Calculations, Tech. rep., Carbon Dioxide Inf. Anal. Cent., Oak Ridge Natl. Lab., US DOE, Oak Ridge, Tenn., 2006.

Proye, A. and Gattuso, J.-P.: Seacarb, an R package to calculate parameters of the seawater carbonate system, 2003.

R Development Core Team: R: A Language and Environment for Statistical Computing, R Foundation for Statistical Computing, Vienna, Austria, ISBN 3-900051-07-0, 2012.

Riley, J. P.: The occurrence of anomalously high fluoride concentrations in the North Atlantic, Deep-Sea Res., 12, 219–220, 1965.

Riley, J. and Tongudai, M.: The major cation/chlorinity ratios in sea water, Chem. Geol., 2, 263–269, 1967.

Robbins, L. L., Hansen, M. E., Kleypas, J. A., and Meylan, S. C.: CO2calc: a user-friendly carbon calculator for Windows, Mac OS X, and iOS (iPhone), USGS Open-File Report, 2010–1280, 17, 2010.

Roy, R. N., Roy, L. N., Vogel, K. M., Porter-Moore, C., Pearson, T., Good, C. E., Millero, F. J., and Campbell, D. M.: The dissociation constants of carbonic acid in seawater at salinities 5 to 45 and temperatures 0 to 45 °C, Mar. Chem., 44, 249–267, 1993.

Sarmiento, J. L., Le Quéré, C., and Pacala, S.: Limiting future atmospheric carbon dioxide, Global Biogeochem. Cy., 9, 121–137, 1995.

Schlitzer, R.: Interactive analysis and visualization of geoscience data with Ocean Data View, Comput. Geosci., 28, 1211–1218, 2002.

Uppström, L. R.: The Boron/Chlorinity Ratio of Deep-Sea Water from the Pacific Ocean, in: Deep-Sea Research and Oceanographic Abstracts, 21, 161–162, Elsevier, 1974.

van Heuven, S., Pierrot, D., Rae, J. W. B., Lewis, E., and Wallace, D. W. R.: MATLAB program developed for CO₂ system calculations, ORNL/CDIAC-105b, Carbon Dioxide Inf. Anal. Cent., Oak Ridge Natl. Lab., US DOE, Oak Ridge, Tenn., 2011.

Weiss, R. F.: Carbon dioxide in water and seawater: the solubility of a non-ideal gas, Mar. Chem., 2, 203–215, 1974.

Wittmann, A. C. and Pörtner, H.-O.: Sensitivities of extant animal taxa to ocean acidification, Nature Climate Change, 3, 995–1001, 2013.

Zeebe, R. E. and Wolf-Gladrow, D. A.: CO₂ in Seawater: Equilibrium, Kinetics, Isotopes, Elsevier, Amsterdam, 346 pp., 2001.