Articles | Volume 15, issue 1
Development and technical paper
04 Jan 2022
Development and technical paper |  | 04 Jan 2022

PyCO2SYS v1.8: marine carbonate system calculations in Python

Matthew P. Humphreys, Ernie R. Lewis, Jonathan D. Sharp, and Denis Pierrot

Oceanic dissolved inorganic carbon (TC) is the largest pool of carbon that substantially interacts with the atmosphere on human timescales. Oceanic TC is increasing through uptake of anthropogenic carbon dioxide (CO2), and seawater pH is decreasing as a consequence. Both the exchange of CO2 between the ocean and atmosphere and the pH response are governed by a set of parameters that interact through chemical equilibria, collectively known as the marine carbonate system. To investigate these processes, at least two of the marine carbonate system's parameters are typically measured – most commonly, two from TC, total alkalinity (AT), pH, and seawater CO2 fugacity (fCO2; or its partial pressure, pCO2, or its dry-air mole fraction, xCO2) – from which the remaining parameters can be calculated and the equilibrium state of seawater solved. Several software tools exist to carry out these calculations, but no fully functional and rigorously validated tool written in Python, a popular scientific programming language, was previously available. Here, we present PyCO2SYS, a Python package intended to fill this capability gap. We describe the elements of PyCO2SYS that have been inherited from the existing CO2SYS family of software and explain subsequent adjustments and improvements. For example, PyCO2SYS uses automatic differentiation to solve the marine carbonate system and calculate chemical buffer factors, ensuring that the effect of every modelled solute and reaction is accurately included in all its results. We validate PyCO2SYS with internal consistency tests and comparisons against other software, showing that PyCO2SYS produces results that are either virtually identical or different for known reasons, with the differences negligible for all practical purposes. We discuss insights that guided the development of PyCO2SYS: for example, the fact that the marine carbonate system cannot be unambiguously solved from certain pairs of parameters. Finally, we consider potential future developments to PyCO2SYS and discuss the outlook for this and other software for solving the marine carbonate system. The code for PyCO2SYS is distributed via GitHub (, last access: 23 December 2021) under the GNU General Public License v3, archived on Zenodo (Humphreys et al.2021), and documented online (, last access: 23 December 2021).

1 Introduction

The ocean absorbs about a quarter of the anthropogenic carbon dioxide (CO2) currently emitted each year (Friedlingstein et al.2020). This absorption is a double-edged sword. Removing CO2 from the atmosphere reduces the impact of these emissions on Earth's climate. However, CO2 uptake causes seawater pH and calcium carbonate mineral saturation states (Ω) to decline through a process termed ocean acidification, which has adverse effects on some marine species and ecosystems (Doney et al.2009).

Exchange of CO2 between the atmosphere and ocean, and the biogeochemical consequences of this process, are governed by a series of equilibrium chemical reactions and parameters collectively known as the marine carbonate system (Millero2000). The core parameters are the substance contents of aqueous CO2, the bicarbonate and carbonate ions formed by its hydration and dissociation (HCO3- and CO32-), and the sum of these three components (dissolved inorganic carbon, TC); total alkalinity (ATDickson1981); the fugacity, partial pressure, or dry-air mole fraction of CO2 in seawater (fCO2, pCO2, or xCO2Weiss1974); and pH (Dickson et al.2015). If any valid pair of these parameters is known, plus auxiliary data including temperature, pressure, salinity, and nutrient contents, then all the other parameters can be calculated (Park1969; Zeebe and Wolf-Gladrow2001).

Many research questions require solving the marine carbonate system from some measured or modelled pair of its parameters. Several software tools have been developed for this purpose such that most scientific software environments and programming languages have a widely accepted marine carbonate system solver (Orr et al.2015). However, there is not yet an established and fully functional tool for the popular scientific programming language Python, although partial solutions exist (e.g. Branson2018). Here, we present PyCO2SYS, a Python package designed to fill this capability gap and provide a robust platform for future developments in calculating marine chemical speciation. Being free, open-source, and working across all major operating systems, a Python package is a highly accessible, desirable, and useful tool.

As its name suggests, PyCO2SYS originates from the existing CO2SYS family of software. The original CO2SYS program for MS-DOS (Lewis and Wallace1998) has been further developed and “translated”, with implementations now available for Microsoft Excel (Pierrot et al.2006; Orr et al.2018; Pierrot et al.2021) and MATLAB/GNU Octave (van Heuven et al.2011; Xu et al.2017; Orr et al.2018; Sharp and Byrne2019; Sharp et al.2020). PyCO2SYS was created as an as-close-as-possible translation of CO2SYS-MATLAB v2.0.5 (Orr et al.2018), but we have since made several additional developments to it. Many of these developments involved reshaping the internal code into a more Pythonic style. These changes did not affect the calculations and so are not discussed further. Other developments added new functionality or made minor differences to the calculated results; these are documented and justified here.

As the original CO2SYS software is so well-established in the research field, we provide a relatively brief summary of the components of PyCO2SYS that are identical to CO2SYS-MATLAB in Sect. 2, before describing the areas where PyCO2SYS differs in more detail in Sect. 3. Equations that were inherited from CO2SYS-MATLAB or taken from the literature are generally reported in appendices rather than being reproduced in these sections. We go on to validate PyCO2SYS in Sect. 4 by examining its internal consistency and by comparing its calculations with another CO2SYS implementation. In Sect. 5, we discuss some nuances of solving the marine carbonate system that were explored during development and compare its computational speed with CO2SYS-MATLAB, before concluding with our perspectives on the outlook for PyCO2SYS and other related software.

2 Methods inherited from CO2SYS

The components of PyCO2SYS that have been inherited directly from CO2SYS-MATLAB v2.0.5 (Orr et al.2018), with only the minimal changes needed to translate to Python plus aesthetic code restructuring, are described in this section.

2.1 Units and pH scales

The abundances of all solutes and total alkalinity provided as arguments to PyCO2SYS or returned from it as results are in units of micromoles per kilogram (µmol kg−1), with kilograms (kg) being of the total solution. This means that they are neither concentrations nor molarity values, which are both per unit volume rather than mass, nor are they molality values, which are per kilogram of H2O. Although sometimes referred to as molinity, the correct term is substance content (IUPAC1997), which we abbreviate to content.

Temperature is in degrees Celsius (C) and salinity is practical salinity, which is dimensionless (Millero et al.2008).

Pressure is in decibars (dbar) and represents the hydrostatic pressure exerted by the overlying water column, consistent with typical oceanographic conductivity–temperature–depth (CTD) measurement reporting. Atmospheric pressure is not included, so pressure is effectively zero in the laboratory and at the sea surface.

The pH can be provided on the free, total, seawater, and/or NBS (now NIST) scale, with [H+] as a substance content, as noted above, and thus in moles per kilogram of solution (mol kg solution−1) (Appendix AZeebe and Wolf-Gladrow2001; Velo et al.2010). In the results, pH is returned on all four of these scales.

2.2 Parameterisations and constants

A notable feature of all CO2SYS software is the variety of different parameterisation options to calculate the various equilibrium constants and some components' total contents from salinity, temperature, and pressure. Which parameterisations the user selects can appreciably alter the results, so these choices should always be explicitly reported.

Roy et al. (1993)Goyet and Poisson (1989)Dickson and Millero (1987)Dickson and Millero (1987)Dickson and Millero (1987)Mehrbach et al. (1973)Mehrbach et al. (1973)Millero (1979)Cai and Wang (1998)Lueker et al. (2000)Mojica Prieto and Millero (2002)Millero et al. (2002)Millero et al. (2006)Millero (2010)Waters and Millero (2013)Sulpis et al. (2020)Schockman and Byrne (2021)

Table 1Parameterisations of the dissociation constants of carbonic acid available in PyCO2SYS and corresponding implicit settings (Table 2).

a Refit of Hansson (1973a, b) data. b Refit of Mehrbach et al. (1973) data. c Refit of Hansson (1973a, b) and Mehrbach et al. (1973) data. d Constants for zero-salinity freshwater. e Including the corrections of Waters et al. (2014).

Download Print Version | Download XLSX

Uppström (1974)Culkin (1965)Lee et al. (2010)Riley and Tongudai (1967)Culkin (1965)Riley (1965)Riley (1965)Morris and Riley (1966)Morris and Riley (1966)Millero (1995)Takahashi et al. (1982)Millero (1983)Millero (1995)Millero (1979)Millero (1979)Millero (1995)Millero (1995)Millero (1983)Dickson (1990b)Li et al. (1969)Millero (1979)Edmond and Gieskes (1970)Yao and Millero (1995)Kester and Pytkowicz (1967)Millero (1983)Millero (1983)Yao and Millero (1995)Sillén et al. (1964)Millero (1995)Millero (1995)Millero (1983)Ingle (1975)Ingle (1975)Takahashi et al. (1982)Millero (1983)Ingle et al. (1973)Ingle (1975)Takahashi et al. (1982)Weiss (1974)Weiss (1974)

Table 2Parameterisations that vary depending on the case of the selected carbonic acid constants (Table 1). P: pressure.

a Depending on user input. b In GEOSECS-Takahashi, phosphate is not included in the definition of total alkalinity; in GEOSECS-Peng, phosphate is included, though the contribution of each species to alkalinity is determined incorrectly based on charge rather than a zero level of protons at pK 4.5. c Includes all dissociation constants for this system: KP1*, KP2*, and KP3* (Appendix B). d Copies the pressure correction for boric acid. e A constant value of 1 is used in this case, i.e. pCO2=fCO2.

Download Print Version | Download XLSX

Some of these options also influence other, seemingly unrelated, parameters of other chemical systems. This is not widely appreciated because this happens internally, hidden within the code. The most influential choice is for the carbonic acid dissociation constants, K1* and K2*, for which there are 17 different options in PyCO2SYS (Table 1). We organise these options into three groups based on their effect on the “hidden” internal parameterisations (Table 2).

  1. Standard case. These are all identical, aside from their varying carbonic acid constants.

  2. GEOSECS cases: GEOSECS-Takahashi and GEOSECS-Peng. GEOSECS-Peng treats phosphate differently with respect to its contribution to alkalinity, and this difference is reported in the results as the “Peng correction”; see Lewis and Wallace (1998) for a more detailed explanation.

  3. Freshwater case. Salinity and other total salt contents (ammonia, borate, calcium, fluoride, phosphate, silicate, sulfate, and sulfide) are set to zero, irrespective of the user inputs.

Khoo et al. (1977)Dickson (1990a)Waters and Millero (2013)Millero (1995)Dickson and Riley (1979)Perez and Fraga (1987)Millero (1995)Clegg and Whitfield (1995)Millero (1995)Yao and Millero (1995)Millero (1995)Takahashi et al. (1982)Peng et al. (1987)Weiss and Price (1980)Weiss (1974)

Table 3Parameterisations that (except where noted) are not influenced by the case of the selected carbonic acid constants (Table 1).

a As selected by the user. b Including the corrections of Waters et al. (2014). c This option was written into the code for CO2SYS-MATLAB v2.0.5 and other versions, but commented out and therefore not directly usable. It is available in CO2SYS-MATLAB v3.2.0.

Download Print Version | Download XLSX

Other internal settings are consistent across all cases (Table 3). These three cases have been present since the original CO2SYS for MS-DOS (Lewis and Wallace1998). That program included only options 1–8 for the carbonic acid dissociation constants (Table 1), with the others being published subsequent to its release. All subsequently added carbonic acid options follow the standard case. While it is beyond the scope of this paper to judge the relative merits of the different options, in general we recommend that one of the standard cases be used unless there is a specific reason for doing otherwise.

In addition to the carbonic acid equilibria, the user has multiple parameterisation options for each of the following: (i) the ratio between total borate and salinity, (ii) the bisulfate dissociation constant KSO4*, and (iii) the hydrogen fluoride dissociation constant KHF* (Tables 2 and 3). However, note that for (i), the user's choice is not respected in the GEOSECS cases, and neither (ii) nor (iii) is included at all in the freshwater case (Table 2). It should also be noted that choices (ii) and (iii) affect pH scale conversions, including of equilibrium constants, which can have a small (but practically negligible) effect on the results.

Equilibrium constants in PyCO2SYS are all stoichiometric rather than thermodynamic and thus denoted with K*. This means that they represent the equilibrium balance of solute substance contents, not of their chemical activities. They are evaluated as follows:

  1. calculated on the pH scale reported in the literature as a function of temperature and salinity at zero in-water pressure;

  2. converted to the seawater pH scale (Appendix A);

  3. corrected to the in situ pressure;

  4. converted to the pH scale indicated by the user's input (Appendix A).

There are some exceptions to the evaluation steps listed above. First, the pH scale conversions (steps 2 and 4) are not applied to KSO4*, KHF*, Ksp*(calcite), Ksp*(aragonite), or K0*. For KSO4* and KHF*, this is because these constants always remain on the free pH scale. The other constants in this list are for equilibria that do not directly involve H+ and are therefore independent of the pH scale. Second, no pressure correction (step 3) is applied to the CO2 solubility factor K0* (Weiss1974). This value and calculations of fCO2, pCO2, and xCO2 are thus valid only for the surface ocean (Sect. 5.3).

In PyCO2SYS, users can also specify their own values for any or all of the equilibrium constants or total salt contents. Any values specified in this way are used as-is throughout PyCO2SYS: no pH scale or pressure corrections are applied, so it is left to the user to ensure that the values are provided on the appropriate pH scale and at the relevant temperature and pressure.

2.3 Input and output conditions

A useful feature of all CO2SYS software that can nonetheless cause confusion is calculations at “input” and “output” conditions; “conditions” refers to temperature and pressure. There is an unhelpful overlap of nomenclature, with input and output used firstly in a programming context to refer to arguments that are passed into functions and returned from them as results, and secondly in a measurement context in which they refer to the temperatures and pressures under which the known parameter pairs are provided and at which results are to be calculated. For clarity, we therefore use the terms “arguments” and “results” in the programming context, while input and output always refer to the measurement context. Thus, we provide values at both input and output conditions as arguments to PyCO2SYS and we receive calculations at both input and output conditions as results from the program.

Input and output conditions are used when measurements were conducted at a different temperature and/or pressure from what the sample would experience in situ or to evaluate the effect of changing these conditions on the solution chemistry. All core carbonate system parameters except for AT and TC are temperature- and pressure-sensitive, so the values of other measured arguments and calculated results may differ between the input and output conditions. For example, measurements might be conducted in a laboratory at 25 C on samples that were collected from several kilometres' depth in the ocean at sub-zero temperatures. In this case, we would provide the measurement conditions (i.e. temperature and pressure in the laboratory) as input arguments and the environmental conditions (i.e. temperature and pressure in the ocean) as output arguments. The corresponding output-condition results from PyCO2SYS then represent the true state of the sample in situ in its environment. The input-condition results are of less environmental interest but may be useful for quality-control purposes.

If calculations are conducted using only in situ values, for example from model output or with the temperature and pressure corrections already applied, then output-condition arguments need not be supplied. Results are then calculated only under the input conditions for computational efficiency.

2.4 Solving the marine carbonate system

We refer to the parameters from which PyCO2SYS can solve the marine carbonate system as the “core” marine carbonate system parameters. These are AT, TC, pH (on any scale), pCO2, fCO2, xCO2, [CO2(aq)], [HCO3-], and [CO32-]. Any pair of these can be provided, except for two of pCO2, fCO2, xCO2, and [CO2(aq)], which would not be valid because these are all directly proportional to each other at a given temperature, salinity, and atmospheric pressure.

Figure 1Overview of the process by which PyCO2SYS and other CO2SYS implementations solve the marine carbonate system (MCS) and calculate other results. Arguments provided by the user are shown as open symbols on a yellow background, while calculations and results use filled symbols. Components under input conditions are shown in light blue, those under output conditions are in red towards the right, and components that are independent of input/output conditions are in dark blue. Any pair of the parameters in the “MCS arguments” box at the top left can be provided, noting that only one of [CO2(aq)], pCO2, fCO2, or xCO2 may be included in a pair. Coupled with user-provided nutrients, total salts calculated from salinity (Totals), and stoichiometric dissociation constants calculated from salinity and input temperature and pressure (K* values), all core MCS parameters are determined (Input MCS results) from the known pair (Appendix C). Other results (e.g. carbonate mineral saturation states, buffer factors) are then calculated from the results under input conditions (Others). If the user provides output-condition temperature and/or pressure values, then the dissociation constants are recalculated under these new conditions, the core MCS is solved again (Output MCS results) from these updated constants (K* values out), the original “totals”, and the now-known AT and TC, which are independent of temperature and pressure. Finally, other results are calculated again from the output-condition results (Others out).


To calculate its results (Fig. 1), PyCO2SYS first determines the unknown core parameters from whichever pair is provided by the user under the input conditions (Appendix C). The parameter pairs that require an iterative solver to find pH (i.e. AT plus TC or one of its components) are solved using a scheme that has been updated from previous versions of CO2SYS (Sect. 3.1). The AT and TC provided or determined under the input conditions are then used to solve the core marine carbonate system again under the output conditions if these have been provided. This is possible because both AT and TC are unaffected by temperature and pressure changes.

Other properties of interest are subsequently calculated from whichever core parameters are most convenient under both input and (if provided) output conditions. These properties include all the individual components of alkalinity (Appendix B), calcite, and aragonite saturation states (Appendix D), as well as various chemical buffer factors (Sect. 3.3.4).

3 New developments in PyCO2SYS

3.1 Solving the alkalinity–pH equation

3.1.1 Automatic differentiation

Solving the alkalinity–pH equation is a critical component of marine carbonate system modelling. Like other implementations of CO2SYS, PyCO2SYS uses the Newton–Raphson method. The general equation is

(1) pH n + 1 = pH n - Δ A T ( pH n , v ) Δ A T ( pH n , v ) ,

where AT=dAT/dpH and

(2) Δ A T ( pH n , v ) = A T ( pH n , v ) - A T ( known ) ,

in which v is any of TC, fCO2, [HCO3-], or [CO32-]. AT(pHn,v) is determined as described in Sects. C2.1 (when v is TC) and C2.3C2.5 (when v is one of fCO2, [CO32-], or [HCO3-]).

Unlike other implementations of CO2SYS, the equations that determine the relative abundances of different chemical species as functions of pH and their total contents (Appendix B) appear only once in PyCO2SYS in what we term the “main chemical speciation function”. While this approach does not alter the calculated results, it does make the software more robust by reducing the opportunity for typographical errors when similar equations are repeated across the code.

The derivative term in Eq. (1) is evaluated from the main chemical speciation function using automatic differentiation, as implemented by the Python package Autograd (Maclaurin2016). Distinct from numerical or symbolic differentiation, the automatic approach breaks down the code to be differentiated into a sequence of individual arithmetic operations (addition, subtraction, etc.) and simple functions (logarithms, exponentials, etc.), then combines the derivatives of these components together using the chain rule. The overall differentials to arbitrary order of complicated functions can thus be evaluated efficiently and are accurate to the computer's precision.

Through our approach, the effect of every component of alkalinity in the main chemical speciation equation is included in the derivative term in Eq. (1). In contrast, some other implementations of CO2SYS use simplified expressions that only include the contributions of carbonate, borate, and water to the total alkalinity. Under typical open-ocean conditions, this makes little practical difference because the simplified equations include the most important components of the seawater solution. However, including every modelled component does make the solver more robust for more unusual solution compositions.

Automatic differentiation is also used to evaluate chemical buffer factors, again ensuring that the influence of every modelled equilibrium system is accurately included. The calculated buffer factors are described in more detail in Sect. 3.3.4.

A further advantage of the automatic differentiation approach is that if the main chemical speciation function is modified in the future, for example to include additional components of alkalinity, then these changes are automatically incorporated into all the alkalinity–pH solvers without needing to modify the various solver functions. In short, our approach ensures that PyCO2SYS calculations will remain internally consistent and reflect the influence of every solute and equilibrium modelled in the main chemical speciation function, even if this function is modified in the course of future development (Sect. 5.5).

3.1.2 Vectorised arguments and solver jumps

PyCO2SYS adjusts how to determine when the alkalinity–pH solver should stop solving for vectorised arguments. In CO2SYS-MATLAB v2.0.5, the solvers continue to iterate and update all values until the change in every element of the array satisfies the ΔpH tolerance threshold (10−4 in CO2SYS-MATLAB, 10−8 in PyCO2SYS). This means that a given set of arguments could return slightly different results depending on what data appear in the other, supposedly independent, elements of the argument arrays. Although negligible for all practical purposes, these differences are detectable in code validation exercises. In PyCO2SYS (and in CO2SYS-MATLAB v3.2.0; Sharp et al.2020) this process has been changed such that each element stops being updated once it has reached the tolerance threshold, independent of the other elements.

The maximum solver jump – which constrains the greatest change in pH possible between solver iterations, thus helping to prevent overshoot – is implemented slightly differently in PyCO2SYS than in other CO2SYS programs. In CO2SYS-MATLAB, any ΔpH values with a magnitude greater than 1 are halved. In PyCO2SYS, the same applies, but any ΔpH values with a magnitude still greater than 1 after halving are decreased to 1 (while preserving the sign). This has a negligible effect on calculations, but it is sometimes detectable in intercomparisons.

3.1.3 pH scale conversions

PyCO2SYS fixes a simplification in earlier CO2SYS implementations regarding how pH scales are converted within the main chemical speciation function. This simplification is noted in the programmer's comments in the relevant CO2SYS-MATLAB functions, carried through from the original MS-DOS implementation (Lewis and Wallace1998): “Though it is coded for H on the total pH scale, for the pH values occurring in seawater (pH > 6) it will be equally valid on any pH scale (H terms negligible) as long as the K Constants are on that scale.”

In short, pH and the equilibrium constants are provided to these functions on the same pH scale as each other – except for KSO4* and KF*, which are always on the free scale (Sect. 2.2). Calculations of all alkalinity components except [HSO4-] and [HF] have therefore always been correct. However, because KSO4* and KF* are always on the free scale, pH must be converted to this scale in order to determine the contributions of [HSO4-] and [HF] to total alkalinity. Other versions of CO2SYS prior to CO2SYS-MATLAB v3.2.0 (Sharp et al.2020) and CO2SYS-Excel v3 (Pierrot et al.2021) assume that the user-selected pH scale is total and thus apply the total-to-free scale conversion (Appendix A) regardless of what the user-selected pH scale actually is.

This simplification makes a negligible difference to calculations at typical seawater pH (because [HSO4-] and [HF] are each on the order of 10−10µmol kg−1 relative to AT on the order of 2000 µmol kg−1) and then only when the user-selected pH scale is not total. But, as implied in the original programmer's note, it can have a noticeable adverse effect under other conditions, such as the low pH values encountered during the acidimetric titrations of seawater used to measure AT. In PyCO2SYS, CO2SYS-MATLAB v3.2.0, and CO2SYS-Excel v3, the correct conversion factor is used based on the user-selected pH scale.

3.2 Initial pH estimates

Like most iterative solvers, the Newton–Raphson method (Sect. 3.1) requires an initial pH value that is near the true value in order to prevent overshoot and guarantee convergence to a root. Previous versions of CO2SYS used 8 as the initial pH estimate in every case. This works well for typical open-ocean seawater but may be less appropriate in niche environments or when modelling acidimetric titrations. Munhoven (2013) found a better initial pH estimate for solving from known AT and TC by considering only the contributions of carbonate and borate species to AT, simplifying the AT equation.

(3) A CB = [ HCO 3 - ] + 2 [ CO 3 2 - ] + [ B ( OH ) 4 - ]

Following Munhoven (2013) and as also implemented elsewhere (e.g. Orr and Epitalon2015), PyCO2SYS and CO2SYS-MATLAB v3.2.0 also take this approach (Appendix F). Furthermore, we have extended it to apply to the pH solvers that use one of the components of TC as the second known parameter, as follows. We note that these extensions are equivalent to those described and discussed in greater detail by Munhoven (2021), although they were added to the PyCO2SYS code in its v1.3.0 release (, Humphreys et al.2020​​​​​​​) before the publication of that study.

3.2.1 Solving from AT and fCO2

For clarity in the equations in this section, we abbreviate [CO2(aq)] as s and [H+] as h. As noted in Appendix C1.2, the approach described here is also used for known parameter pairs of AT plus any of pCO2, xCO2, or [CO2(aq)].

First, fCO2 is converted to s using Eq. (C5). Carbonate–borate alkalinity (ACB) as a function of s and h is

(4) A CB ( h , s ) = K 1 * s ( h + 2 K 2 * ) h 2 + K B * T B h + K B * .

This can be rearranged into a third-order polynomial in h:

(5) P s ( h , s ) = h 3 + h 2 g 2 ( s ) + h g 1 ( s ) + g 0 ( s ) = 0 ,

with the following.


Following a scheme equivalent to Munhoven (2013), the initial h value is determined by

(9) h 0 ( s ) = 10 - 3 for  A T 0 h min + - P s ( h min ) g 2 2 - 3 g 1 for  A T > 0 .

Negative ACB is impossible because both terms in Eq. (4) are always positive, so the equations given above cannot be applied if AT is indeed negative (e.g. after the alkalinity end point in an acidimetric titration). The default h0 of 10−3 mol kg−1, corresponding to a pH of 3, is therefore used in this case. Otherwise, hmin in Eq. (9) is found following Munhoven (2013).

(10) h min = - g 2 + g 2 2 - 3 g 1 / 3 for  g 2 < 0 - g 1 / g 2 + g 2 2 - 3 g 1 for  g 2 0

When AT is positive, the square-rooted term g22-3g1 is always greater than zero, and thus hmin has a real value. However, there is an additional constraint: ACB cannot be greater than 2TC+TB (Munhoven2013). If AT is actually greater than this limit, then we use a default h0 of 10−7 mol kg−1 instead (pH 7).

3.2.2 Solving from AT and [HCO3-]

For clarity in the equations in this section, we abbreviate [HCO3-] as b and [H+] as h.

Carbonate–borate alkalinity as a function of b is

(11) A CB ( h , b ) = b + 2 K 2 * b h + K B * T B h + K B * .

This can be rearranged into a second-order polynomial in h:

(12) P b ( h , b ) = h 2 g 2 ( b ) + h g 1 ( b ) + g 1 ( b ) = 0 ,

with the following.


The initial h value is estimated following

(16) h 0 ( b ) = - g 1 - g 1 2 - 4 g 0 g 2 2 g 2 for  b < A T 10 - 3 for  b A T .

When b<AT, the square-rooted term g12-4g0g2 is always positive, and thus h0(b) has a real value. Otherwise, b can only be greater than AT if the negative components of AT such as [H+] are dominant, as happens at low pH. The default initial pH estimate used by PyCO2SYS in that case is therefore 3.

3.2.3 Solving from AT and [CO32-]

For clarity in the equations in this section, we abbreviate [CO32-] here as c and [H+] as h.

Carbonate–borate alkalinity as a function of c is

(17) A CB ( h , c ) = c h K 2 * + 2 c + K B * T B h + K B * .

This can be rearranged into a second-order polynomial in h:

(18) P c ( h , c ) = h 2 g 2 ( c ) + h g 1 ( c ) + g 0 ( c ) = 0 ,

with the following.


The initial h value is estimated following

(22) h 0 ( c ) = - g 1 + g 1 2 - 4 g 0 g 2 2 g 2 for  A T > 2 c + T B 10 - 3 for  A T 2 c + T B .

When 2c+TB<AT, the square-rooted term g12-4g0g2 is always positive, and thus h0(c) has a real value. Otherwise, 2c+TB can only be greater than AT if the negative components of AT such as [H+] are dominant, as happens at low pH. The default initial pH estimate used by PyCO2SYS in that case is therefore 3.

3.3 New calculations, components, and constants

3.3.1 Additional alkalinity components

The contributions of ammonia and bisulfide to alkalinity (Cai et al.2017; Xu et al.2017) plus the ability to solve from carbonate and/or bicarbonate ion content have been added in collaboration with Sharp et al. (2020) to ensure consistency between PyCO2SYS and CO2SYS-MATLAB v3.2.0. However, the GEOSECS alkalinity definition did not account for these species, so if using one of the GEOSECS options for the carbonic acid constants (Table 1) then users should be sure to set their total contents to zero for GEOSECS-compatible results. If values are provided, then they will be included in the alkalinity equation just as for the non-GEOSECS cases.

The total substance contents and stoichiometric dissociation constants for up to two additional acid–base systems that contribute to total alkalinity can be provided as arguments to PyCO2SYS and are part of its speciation model. The effects of these extra components are automatically incorporated into all PyCO2SYS calculations, including the iterative pH solvers (Sect. 3.1), buffer factors (Sect. 3.3.4), and uncertainty propagation (Sect. 3.6). These extra components are modelled following Sharp and Byrne (2020), as described in Appendix B11. No corrections of any sort (e.g. for pressure or pH scale; Sect. 2.2) are made to the dissociation constants for these user-defined additional components within PyCO2SYS; the user must ensure that they are already suitable for the conditions being analysed and on the user-indicated pH scale.

3.3.2 Gas constant

Previous versions of CO2SYS used an old value for the universal gas constant (R) of 8.31451 J mol−1 K−1. PyCO2SYS uses the 2018 CODATA-recommended value by default instead (i.e. 8.314462618 J mol−1 K−1), consistent with CO2SYS-MATLAB v3.2.0 and CO2SYS-Excel v3. This has a minor effect on conversions between pCO2, fCO2, and xCO2 (less than 10−4 %), as well as on the pressure corrections for the equilibrium constants (less than 10−3 % at 5000 dbar). It is detectable in comparisons with other versions of CO2SYS, but it is of no practical consequence.

3.3.3 Substrate : inhibitor ratio

Like CO2SYS-MATLAB v3.2.0, PyCO2SYS calculates the “substrate : inhibitor ratio” of Bach (2015), which quantifies the balance between the availability of a substrate for calcification (i.e. HCO3-) and the inhibition of calcification by H+ (Eq. D2).

3.3.4 Buffer factors

A buffer factor quantifies the sensitivity of a certain marine carbonate system parameter to a change in another parameter. Best known is the Revelle factor, which is the ratio of the fractional change in pCO2 corresponding to a fractional change in TC at constant AT (Revelle and Suess1957). Frankignoulle (1994) derived a broader set of buffer factors for the marine carbonate system, quantifying the responses of several different parameters to changes in TC and AT; these were later rediscovered by Egleston et al. (2010) and further extended by Hagens and Middelburg (2016). PyCO2SYS calculates the buffer factors of Egleston et al. (2010) and uses the nomenclature of that paper.

Closely related to these buffer factors, Frankignoulle et al. (1994) introduced the factor ψ, which quantifies the change in TC required to return to the original seawater pCO2 after the action of calcification (which reduces AT and TC in a 2 : 1 ratio) or CaCO3 dissolution (the reverse). Humphreys et al. (2018) introduced the “isocapnic quotient” (Q), which is the ratio of AT to TC change that does not affect seawater pCO2, thus generalising the concept of ψ for application to all biogeochemical processes that affect AT and TC (denoted ϕ). PyCO2SYS calculates both ψ and Q, the latter of which can be used to calculate ϕ for any biogeochemical process (Humphreys et al.2018).

PyCO2SYS offers two independent ways to evaluate the various buffer factors of the marine carbonate system: with explicit equations and by automatic differentiation. The latter is used by default.

The “explicit” approach follows equations reported in the literature (Frankignoulle et al.1994; Egleston et al.2010; Humphreys et al.2018), noting that the typographical errors in Egleston et al. (2010) identified in several studies (e.g. Orr2011; Álvarez et al.2014; Richier et al.2018; Orr et al.2018) have been corrected. In general, these equations do not include the effect of species beyond the carbonate, borate, and water contributions to total alkalinity, except that the buffer factors of Egleston et al. (2010) were extended to include phosphate and silicate effects by Orr et al. (2018).

The “automatic” approach uses automatic differentiation to find the derivative necessary to evaluate each buffer factor. The appropriate derivatives are taken from the functions that calculate a third carbonate system parameter from a known pair (Appendix C). All species modelled in the main chemical speciation function are therefore included, including any extra alkalinity components (Sect. 3.3.1), and typographical errors from the literature cannot influence these calculations. The details of the derivatives used are provided in Appendix E.

Of the buffer factors, only the Revelle factor was included in previous versions of CO2SYS. It was evaluated using finite-central-difference derivatives, which is replicated as the explicit option in PyCO2SYS (with the corrections described in Appendix G). However, as for all other buffer factors, the Revelle factor calculation uses automatic differentiation by default. To calculate the Revelle factor using a mathematical approach equivalent to the explicit calculation of the other buffer factors, one could calculate γTC of Egleston et al. (2010) (see Appendix E1) with the explicit approach and then use Eq. (E7).

3.3.5 Atmospheric pressure

For conversions between pCO2, fCO2, and xCO2, atmospheric pressure is assumed to be 1 atm by default, and it remains fixed at this value in CO2SYS-MATLAB and CO2SYS-Excel. However, in PyCO2SYS, the user can also specify a value other than 1 atm if necessary. Different values can be provided for input and output conditions (Sect. 2.3).

Atmospheric pressure can have a non-negligible effect on calculations in some regions: for example, over much of the Southern Ocean, atmospheric pressure is typically 3 % lower than the global mean, corresponding to a 10 µatm reduction in pCO2 and fCO2 relative to the values calculated at 1 atm (Orr et al.2017).

This optional argument is only intended for modelling the effects of variations in atmospheric pressure on samples from the surface ocean or in the laboratory. It is not suitable for determining interior ocean pCO2, fCO2, and xCO2 values that are corrected for the pressure of the overlying water column. This separate issue is discussed further in Sect. 5.3.

3.4 No-solve modes

As well as solving from a pair of parameters, PyCO2SYS can be run with one or no marine carbonate system parameter arguments.

If no parameters are provided, then PyCO2SYS returns all the equilibrium constants and total salt contents that are calculated from temperature, pressure, and salinity (Sect. 2.2) without actually using these to do any further computations.

If one parameter is provided, then the results that can be computed with that parameter alone are returned. This applies to pH, pCO2, fCO2, xCO2, and [CO2(aq)], as follows.

The pH can be converted between the different scales without knowledge of a second carbonate system parameter. Therefore, if pH alone is provided to PyCO2SYS, it is converted to every pH scale under the input conditions (Appendix A). Conversion to a different temperature and/or pressure does require solving the carbonate system (Fig. 1), so output-condition values are not calculated.

Seawater pCO2, fCO2, xCO2, and [CO2(aq)] can also be interconverted without knowledge of a second carbonate system parameter (Appendix C1.2). Therefore, if any of these parameters alone are provided to PyCO2SYS, all the others are calculated under the input conditions. If an output-condition temperature is provided, then pCO2 is also adjusted to the new temperature following Takahashi et al. (2009), and all others in this set of parameters are calculated under output conditions from the new pCO2 value.

3.5 Multidimensional arguments

All arguments to PyCO2SYS, including settings, can be multidimensional. A combination of scalar and multidimensional arguments can be provided, with the latter formatted as NumPy ndarrays (Harris et al.2020). Results that depend only on scalar arguments are themselves scalar, while results depending on multidimensional inputs are “broadcasted” into consistently shaped arrays (Fig. 2). The code is optimised to efficiently compute across multidimensional arrays following the approach of CO2SYS-MATLAB since its v1.1 (van Heuven et al.2011). However, all multidimensional arrays in CO2SYS-MATLAB are flattened into one-dimensional vectors and returned in the results in that same format.

Figure 2Schematic representation of broadcasting array shapes with NumPy in PyCO2SYS. (a) Two of the arguments to PyCO2SYS are provided as arrays, each containing 11 different values for TC and temperature. Other arguments could be similarly shaped vectors or single scalar values. (b) If the array arguments were all provided as one-dimensional rows, then the calculated results (e.g. aragonite saturation state) would also be one-dimensional rows. Each element of the results array corresponds to the element in the same position in each argument array. For scalar arguments, the same value is used across each result array. (c) If the array arguments are provided as a mixture of rows and columns, then the results are calculated on a broadcasted grid including every combination of the arguments' elements. The same principle applies to arguments and results of arbitrarily higher dimensionality.


3.6 Uncertainty propagation

Propagating the uncertainty in an argument through to a result requires knowing the derivative of the result with respect to the argument. Uncertainty propagation is available for a subset of the arguments in the original MS-DOS CO2SYS (Lewis and Wallace1998) and was added to the Excel and MATLAB implementations more recently (Orr et al.2018). However, while much of the code to solve the marine carbonate system in PyCO2SYS has been directly inherited from CO2SYS-MATLAB, its implementation of uncertainty propagation differs.

PyCO2SYS evaluates the derivatives using a finite-forward-difference approach. We use finite differences rather than automatic differentiation here because the latter, while possible, is computationally inefficient to apply over the entire PyCO2SYS program. We use forward- rather than central-difference derivatives because the former can be safely evaluated at zero for variables for which negative values are impossible (e.g. salinity). The derivative of a result r with respect to an argument a is thus calculated.

(23) r ( a ) a r ( a + Δ a ) - r ( a ) Δ a

The value of Δa is fixed for each argument (Appendix H). Different values for different arguments are necessary because some arguments can differ by over 20 orders of magnitude from others. If Δa is too large, then the derivative may be inaccurate because the equations governing the marine carbonate system are non-linear, but if Δa is too small, then the derivative may be inaccurate due to the limitations of solver tolerance and computer precision. We therefore tested a range of Δa values for each variable under typical open-ocean conditions and selected an appropriate value between these extremes (e.g. Fig. 3). The full list of Δa values is provided in Table H1.

Figure 3An example figure used to select a suitable Δa value for uncertainty propagation, in this case for ΔK1*. The pH was calculated from AT and TC, and then the value of K1* was incremented by the amounts shown on the horizontal axis; the vertical axis is for the corresponding gradient calculated from the pH response, shown by the red curve. The perpendicular blue lines show the Δa value selected in this case (i.e. 10−12), which falls within the flat section towards the centre of the figure. To the left of this (i.e. at higher Δa), the upwards curvature of the red line is due to non-linearity, while the erratic deviations to the right (i.e. at lower Δa) are due to solver tolerance and computer precision limitations.


PyCO2SYS can conveniently obtain derivatives of all its results with respect to all of its arguments and also with respect to all parameters that are normally calculated internally from temperature, pressure, and/or salinity, such as equilibrium constants and total salt contents.

The derivatives are calculated by a function that wraps the entire PyCO2SYS program, rather than by adding extra internal variables that keep track of the effects of differences in to the arguments, as has been implemented elsewhere (e.g. Orr et al.2018). The PyCO2SYS approach means that if the main program is producing valid results, then the derivatives can also be considered reliable without needing to verify some separate calculation mechanism.

To determine the overall uncertainty in each result, the uncertainty components from different arguments are combined using

(24) σ 2 ( r ) = i r a i 2 σ 2 ( a i ) ,

where σ is the uncertainty as a standard deviation (thus, σ2 is a variance). However, Eq. (24) is only valid if the uncertainties in all arguments are independent from each other. Propagation of co-varying uncertainties can still be carried out with PyCO2SYS because, as noted above, the derivative of any result with respect to any argument can be calculated. The user can therefore assemble the Jacobian matrix of partial derivatives needed to propagate any arbitrary set of co-varying argument uncertainties through to any result (JCGM2008).

4 Validation

There are no “certified” results of marine carbonate system calculations against which software like PyCO2SYS can be validated. But we can test its internal consistency, and we can compare its results with the calculations of other programs and values reported in the literature.

PyCO2SYS is developed and hosted on GitHub (, with releases archived on Zenodo (Humphreys et al.2021). Every validation test described in this section is built into PyCO2SYS's test suite; therefore, these tests are executed automatically by GitHub's continuous integration service every time the code is updated. Were any test to fail, an email report would be sent to us, the developers, and the failure displayed publicly in a badge on the GitHub repository's public web page (Fig. 4). Updates to PyCO2SYS are made in a developmental branch of the repository, and the tests must all pass before these changes may be incorporated into the main branch and publicly released in a new version. All validation tests described below were run with PyCO2SYS v1.8.0 (, Humphreys et al.2021), but these protocols should ensure that the quantitative statements made here will hold true as the code continues to be developed.

Figure 4The status badge for the validation tests, which are publicly visible at PyCO2SYS's GitHub repository (, last access: 23 December 2021), when the current version of the code (a) passes every test or (b) fails any test.


For all versions of PyCO2SYS up to v1.8.0, the test suite runs on Python v3.7, 3.8, and 3.9. Other versions of Python may also work but are untested.

4.1 Internal consistency

4.1.1 Round-robin test

In a “round-robin” test, we first determine all of the core carbonate system parameters from one pair and then solve the system again using every possible pair of determined parameters. Under typical seawater conditions, we find the same results for every parameter pair to within better than the tolerance of the iterative pH solvers (i.e. 10−8 in pH). The maximum absolute difference in each parameter across all possible input pair combinations is acceptably small (Table 4).

Table 4Results of an example round-robin test with PyCO2SYS with default parameterisation options. Other conditions: salinity 33, temperature 22 C, pressure 1234 dbar, total silicate 10 µmol kg−1, total phosphate 1 µmol kg−1, total ammonia 2 µmol kg−1, total sulfide 3 µmol kg−1. The pH solver tolerance in PyCO2SYS is 10−8 in terms of pH.

Download Print Version | Download XLSX

4.1.2 Buffer factors

If we include only the solution components that appear in the explicit equations for the buffer factors (i.e. zero nutrients and total salts, except for TB) then we can compare these results with the automatic values (Sect. 3.3.4). Under a range of typical seawater conditions, we find that the differences between these two calculation approaches are totally negligible: on the order of 10−12 % for the Egleston et al. (2010) buffers, 10−9 % for ψ and Q, and 10−7 % for the Revelle factor. The Revelle factor is less well-matched because its explicit value is computed using a finite-difference scheme (for consistency with CO2SYS-MATLAB), which is inherently less accurate than using a direct equation.

Typically, one would not set the total salt contents to zero when computing buffer factors with the default automatic approach. As a consequence, differences between the explicit and automatic buffer factors may be larger than described above but still practically negligible: keeping nutrients at zero but using TSO4 and TF calculated from a salinity of 35, we find that the automatic buffer factors change such that their differences with the corresponding explicit buffer factors increase to the order of 0.01 %.

4.1.3 Uncertainty propagation simulations

The propagation of independent uncertainties using forward-difference derivatives (Sect. 3.6) is tested by comparison with Monte Carlo simulations for all equilibrium constants and all known parameter pair combinations. In every case, the uncertainty determined from the simulations (n=104) as a standard deviation is either within 3 % of the directly calculated value if the latter is non-zero or negligibly small if it is zero (absolute value less than 10−10). The 3 % cutoff is relatively high because of the relatively small number of simulations; the cutoff can be reduced if a greater number of simulations is used, but then the computation time for the test suite becomes impractically long.

4.2 Comparison with other CO2SYS software

We used CO2SYS-MATLAB v2.0.5 (Orr et al.2018) as the main alternative software to compare our results with. PyCO2SYS was originally created as an as-close-as-possible Python translation of this particular version, so any differences in the results should be both understood and intentional. Its predecessor, CO2SYS-MATLAB v1.1 (van Heuven et al.2011), was included in the software intercomparison study of Orr et al. (2015). Indeed, it was selected as the reference software to test the others against. CO2SYS-MATLAB v2.0.5 differs from v1.1 only in that it contains one additional parameterisation for the carbonic acid dissociation constants plus some extra internal variables associated with uncertainty propagation. Comparing PyCO2SYS with CO2SYS-MATLAB v2.0.5 therefore also shows PyCO2SYS's performance and reliability in the context of the wider set of software tested and discussed by Orr et al. (2015).

However, these CO2SYS-MATLAB versions do not permit solving with either carbonate or bicarbonate ion content as a known parameter, nor do they include ammonia or sulfide speciation. They also lack the parameterisations of Sulpis et al. (2020) and Schockman and Byrne (2021) for the carbonic acid dissociation constants (options 16 and 17 in Table 1) and the parameterisation of Waters and Millero (2013) for bisulfate dissociation (Table 3). We therefore also tested PyCO2SYS against CO2SYS-MATLAB v3.2.0 (Sharp et al.2020), which does include all these options.

4.2.1 Temperature–salinity–pressure parameterisations

All equilibrium constants and total salt contents, calculated from salinity, temperature, and pressure, are virtually identical (absolute tolerance 10−12, relative tolerance 10−16, in pK* values or in µmol kg−1) to those in both CO2SYS-MATLAB v2.0.5 and v3.2.0. These tests are run across a range of practical salinity from 0 to 50, temperature from −1 to 50 C, and pressure from 0 to 105 dbar, including values of exactly zero in every case. Every pH scale and parameterisation option is included (Tables 1 and 2).

4.2.2 Solving the marine carbonate system

If PyCO2SYS is adjusted to match CO2SYS-MATLAB v2.0.5, i.e. if the following points are true, then the differences between PyCO2SYS and CO2SYS-MATLAB calculations are virtually zero (no greater than 10−10 %, excluding the Revelle factor as noted above):

  1. approximate slopes are used for the pH solvers, including only carbonate–borate–water alkalinity, instead of using automatic differentiation to determine these exactly (Sect. 3.1.1);

  2. pH solver tolerance is set to 10−4 instead of 10−8 (Sect. 3.1.2);

  3. the original approach to prevent overshoot from solver jumps in pH that are too great is used (Sect. 3.1.2);

  4. the iterative pH solver continues updating all elements until all pH changes fall beneath the tolerance threshold (Sect. 3.1.2);

  5. the pH scale conversion simplification is reinstated (Sect. 3.1.3);

  6. initial pH guesses are always set to 8 instead of using our extended Munhoven (2013) approach (Sect. 3.2);

The Revelle factor is an exception, but this is due to minor errors in its encoding in CO2SYS-MATLAB (Appendix G). If we replicate these errors in PyCO2SYS, then we do return virtually identical Revelle factor values.

If the adjustments above, other than fixing the pH scale conversion simplification, are not made, then the differences between PyCO2SYS and CO2SYS-MATLAB v2.0.5 are up to the order of 10−5 %: greater, but still negligible for all practical purposes.

Fixing the pH scale conversion simplification too (Sect. 3.1.3) makes no difference to calculations for which the user-defined input pH scale is total but causes discrepancies between PyCO2SYS and CO2SYS-MATLAB v2.0.5 of up to 50 % in the “free” hydrogen ion content and 10−2 % in other results when other input pH scale options are selected. The differences are amplified at low pH, as the assumptions of the pH scale conversion simplification do not hold (Sect. 3.1.3).

Repeating the exercise above for CO2SYS-MATLAB v3.2.0 has similar results, with differences negligible for all practical purposes. Only adjustments 1, 2, and 3 from the list above need to be made to PyCO2SYS in this case. With PyCO2SYS fully adjusted to match CO2SYS-MATLAB v3.2.0, differences in calculated values are still mostly less than 10−10 % and with one exception all less than 10−6 %. The exception, a difference still less than 10−3 %, is for the aqueous CO2 content under a limited set of input conditions and only with the new known parameter pair combinations added since CO2SYS-MATLAB v2.0.5. It arises because there are several different ways to calculate [CO2(aq)]: by difference from known TC, [HCO3-], and [CO32-]; from any one of these three variables, [H+], and K1* and K2* equilibrium constants using the equations in Appendix C (Sects. C2.6, C2.11, and C2.12); or from fCO2 or pCO2 and the CO2 solubility constant (K0*). While these approaches are identical in theory, in practice they return different results due to the limitations of solver tolerance and floating-point precision. PyCO2SYS and CO2SYS-MATLAB do not always use the same approach to calculate [CO2(aq)] in each situation (this also varies between CO2SYS-MATLAB versions), hence their greater – but still negligible – differences from each other. Whatever the known parameter pair, PyCO2SYS always follows the principles that (i) the values of parameters provided as arguments by the user should never be overwritten with recalculations, and (ii) the final unknown from TC, [CO32-], [HCO3-], and [CO2(aq)] should always be calculated from the other three by addition or by difference as appropriate.

4.2.3 Uncertainty propagation comparisons

PyCO2SYS reproduces all the derivatives reported by Orr et al. (2018) in their Tables 2 and 3 to within 10−3 % under the same input conditions and all the propagated uncertainties reported by Orr et al. (2018) in their Table 4 to within 10−4 %. We consider all these differences to be negligible.

Across all combinations of optional parameters, mean uncertainties in AT, TC, pCO2, fCO2, [HCO3-], [CO32-], [CO2(aq)], Ω(calcite), Ω(aragonite), and xCO2 propagated from the standard values suggested by Orr et al. (2018) are within 0.5 % of the corresponding uncertainty values calculated with CO2SYS-MATLAB v3.2.0 under the same input conditions. Greater differences in uncertainties calculated under output conditions arise because CO2SYS-MATLAB does not propagate the uncertainties from input-condition equilibrium constants through to output-condition results.

4.3 Simulated seawater titration

PyCO2SYS can be used to reproduce the closed-cell seawater titration datasets simulated by Dickson (1981). Each simulated dataset contains pH values for a seawater sample as it is titrated with incremental HCl additions across a pH range from approximately 8 to 3.

Dickson (1981) specified exact values for all stoichiometric equilibrium constants. PyCO2SYS allows these to be provided instead of being calculated internally from temperature and salinity (Sect. 2.2). The titration is then simulated by calculating how AT should change through the titration due to acid addition, accounting for dilution of AT, TC, and all other dissolved solutes by acid addition, and then solving the carbonate system for pH from the so-determined AT and TC. Being tested here is the ability to solve for pH from known AT and TC across a wide range of pH and AT values, including negative AT.

The first titration dataset, without phosphate, is reproduced perfectly by PyCO2SYS to the number of decimal places reported by Dickson (1981). The second titration, with 10 µmol kg−1 of total phosphate included, is reproduced perfectly by PyCO2SYS with the exception of three values at different titrant masses.

  • 0.45 g: pH either 6.588221 (Dickson) or 6.599221 (PyCO2SYS)

  • 0.60 g: pH either 6.366846 (Dickson) or 6.366486 (PyCO2SYS)

  • 1.25 g: pH either 5.549957 (Dickson) or 5.549951 (PyCO2SYS)

The other 48 data points in this titration agree perfectly. The noted discrepancies occur in non-consecutive data points and are therefore unlikely to all be associated with an error in a particular equilibrium. Coupled with the nature of the differences (underlined above), which is one or two specific digits switched or replaced rather than the entire number being different, we conclude that these differences most likely represent minor typographical errors and therefore that PyCO2SYS does accurately reproduce these simulations in full.

5 Discussion

5.1 Initial pH estimates

The aim of our revised scheme for initial pH estimates, following Munhoven (2013), was to find values that were closer to the final solution across a wide range of pH, thus providing a more suitable starting point for the iterative solvers and thereby reducing the number of iterations required to converge at the solution.

Figure 5Initial estimates (solid lines) and final solutions (dashed lines) of pH from known parameter pairs of total alkalinity (2.3 mmol kg−1) with a range of values for (a) dissolved inorganic carbon (TC), (b) aqueous CO2 fugacity, (c) bicarbonate ion content, and (d) carbonate ion content. The initial estimates track the final solution very closely across the range of typical seawater conditions. This is expected because these estimates were derived under the assumption that the carbonate and borate contributions are dominant in total alkalinity (Sect. 3.2), as is true for typical seawater. The default high and low pH values of 10 and 3 used where the initial estimate equations are not valid for the argument values (Eqs. 16 and F6) appear as flat sections in (a) and (c), respectively.


We find that the initial pH estimates determined according to the scheme described in Sect. 3.2 do follow a pattern similar to the final solutions across wide ranges of argument values, including at the extremes at which the initial-estimate equations become invalid and default pH values are used instead (Fig. 5). The number of iterations required to fall beneath the solver's tolerance threshold (10−8 in pH) is also reduced compared with the original approach of always using an initial pH of 8. Indeed, for typical ocean conditions we find that the iterative solver often does not alter the initial estimate at all. Suitable starting points for the iterative solvers are clearly being found.

5.2 Parameter pairs with multiple solutions

It is not strictly true that the marine carbonate system can always be solved from any pair of its parameters. Some combinations have multiple solutions. For example, both the AT[CO32-] and TC[HCO3-] pairs can correspond to two different pH values (Deffeyes1965; Zeebe and Wolf-Gladrow2001; Munhoven2021). In this section, we show how PyCO2SYS is designed to return the root corresponding to typical seawater. However, it is important to realise that these alternative pH values are real solutions that could be made up in the laboratory or be found in nature; they are not simply mathematical anomalies to be ignored. We therefore used PyCO2SYS to explore the compositions of these alternatives.

5.2.1 Total alkalinity and carbonate ion content

The iterative AT–pH solvers can be thought of as working by evaluating AT at a sequence of different possible pH values until the pH that returns the true AT is found. This pH is known as the “root” of the AT–pH equation. The difference between the true AT and these estimates from pH is the “residual” alkalinity, which is zero at the root. We find that the equations for initial pH estimates and final pH values have very similar roots and similar residuals in the region around these roots (Fig. 6). This similarity is why the initial pH estimates provide such suitable starting points for the final solvers.

Figure 6Residuals between known AT (2.3 mmol kg−1) and (i) carbonate–borate alkalinity (solid lines; ACB) from Eqs. (F1), (4), (11), and (17), as well as (ii) total alkalinity (dashed lines; AT) from Eq. (B1), calculated across a range of pH, with a second known parameter of (a) dissolved inorganic carbon (2.15 mmol kg−1), (b) CO2 fugacity (600 µatm), (c) bicarbonate ion content (2011 µmol kg−1), and (d) carbonate ion content (116 µmol kg−1), all at a salinity of 35 and temperature of 15 C. Each possible pH value returns a different residual alkalinity, and the true pH root is where the residual alkalinity is zero. Both the initial estimates and the final solutions find this zero-residual pH root using the ACB and AT equations, respectively (Sect. 3.1.1 and 3.2). The similarity between the ACB and AT residual curves, particularly around zero residual alkalinity, shows that the initial estimates provide excellent starting values for the subsequent iterative solvers. In (d), the final iterative solver has two possible roots, where residual alkalinity is zero. However, the initial estimate has only one root, corresponding to the lower-pH final root. This ensures that the final solver will always converge to the lower-pH root, which is usually appropriate for the seawater system.


For the AT[CO32-] parameter pair, there are generally two real pH roots and thus two possible equilibrium states of the marine carbonate system (Fig. 6d). We used PyCO2SYS to conceptualise the two pH roots for the AT[CO32-] parameter pair, as follows. The lower-pH root corresponds to typical seawater: a relatively high-TC system, wherein bicarbonate ions are the main component of TC, and carbonate alkalinity ([HCO3-]+2[CO32-]) is the main component of AT. The higher-pH root corresponds to a low-TC system, wherein virtually all of TC is in the form of carbonate ion, and AT is dominated by non-carbonate species (Fig. 7).

Figure 7Main components of (a) total alkalinity (AT) and (b) dissolved inorganic carbon (TC) at the two possible pH roots for a known parameter pair of AT (2300 µmol kg−1) and [CO32-] (120 µmol kg−1). The low-pH root (left) represents typical seawater, with relatively high TC (2143 µmol kg−1) and both AT and TC dominated by the bicarbonate ion (HCO3-). The high-pH root (right) has the same AT and [CO32-], but AT is dominated by hydroxide (OH), and TC is much lower (122 µmol kg−1) and comprised almost entirely of CO32-. These calculations were carried out at 15 C, with a practical salinity of 35 and zero nutrients. If nutrients were present, then like borate (B(OH)4-) they would have different contributions to AT at the different pH roots; pH is on the total scale (Appendix A).


Which root the solver finds depends on the initial pH estimate and the residual alkalinity–pH slope at that point (Eq. 1). This is an advantage of the improved initial pH estimates in PyCO2SYS when working with seawater and similar systems: the initial-estimate equation has only a single real root (Fig. 6d). Because the initial estimate is based on equations for a system that only includes carbonate and borate alkalinity (Sect. 3.2.3), the carbonate system contribution to total alkalinity will always dominate, so the single root of the initial estimate will coincide with the lower-pH true root, which is appropriate for seawater. The solver will thus more robustly find the correct root each time.

In typical open-ocean work this is largely academic: the true pH is typically around 8 and the higher root greater than 10, so a constant initial pH estimate of 8 would also return the correct root. But in more unusual environments, the new algorithm introduced here could help ensure that the solver identifies the correct root. It is possible for the user to specify a different initial pH estimate to control which root PyCO2SYS obtains (as we did to create Fig. 7).

5.2.2 Dissolved inorganic carbon and bicarbonate ion content

As noted previously (e.g. Zeebe and Wolf-Gladrow2001), there are also two possible pH solutions for the TC[HCO3-] parameter pair. We conceptualise these roots as follows: the remaining portion of TC not accounted for by HCO3- is either dominantly composed of CO32- if the solution's pH is closer to pK2* (i.e. higher) or of CO2(aq) if the pH is closer to pK1* (i.e. lower).

Solving from TC and [HCO3-] is more straightforward than from AT and [CO32-] because the unknown pH can be determined from a second-order polynomial, which can be calculated directly using the quadratic formula rather than needing to use an iterative solver. Here, the root found does not depend upon the value of some initial pH estimate. Instead, the quadratic formula generates two possible roots, which must be chosen between. The usual approach, as advised by e.g. Zeebe and Wolf-Gladrow (2001), is to take the higher-pH root, and this is the default behaviour of PyCO2SYS. However, PyCO2SYS can be set to return the other root instead, which we used to illustrate the differing chemistry of the two possibilities (Fig. 8). Munhoven (2021) discusses root-selection strategy for this parameter pair combination in more detail.

Figure 8Main components of (a) total alkalinity (AT) and (b) dissolved inorganic carbon (TC) at the two possible pH roots for a known parameter pair of TC (2100 µmol kg−1) and bicarbonate ion content ([HCO3-]; 1900 µmol kg−1). The high-pH root (left) represents typical seawater, in which most of the TC not accounted for by [HCO3-] is composed of [CO32-]; AT is relatively high (2364 µmol kg−1) and fCO2 low (331 µatm). In the low-pH root (right), the non-HCO3- portion of TC is instead dominated by CO2(aq); alkalinity is lower (1932 µmol kg−1) and fCO2 high (5008 µatm). These calculations were carried out at 15 C, with a practical salinity of 35 and zero nutrients; pH is on the total scale (Appendix A).


5.3 Pressure corrections for pCO2

In PyCO2SYS, pCO2 (and by extension, fCO2 and xCO2) is always evaluated at a total pressure near 1 atm; it is not corrected for the pressure of the overlying water column (Sect. 2.2). This approach is consistent with all existing implementations of CO2SYS. In practice, it means that these values represent the approximate pCO2 that seawater would have if it were brought to the surface ocean without changing the solution composition – “approximate” because this calculation should use potential temperature, rather than in situ temperature, to retrieve the true value expected after adiabatic decompression (Orr and Epitalon2015). PyCO2SYS does not calculate potential temperature, but this could be provided by the user in place of in situ temperature.

Although a pressure correction for pCO2 (i.e. a pressure correction for K0* and the fugacity factor; Appendix C1.2) is theoretically possible (Weiss1974; Orr and Epitalon2015), it could be argued that this is unnecessary. First, the vast majority of pCO2 measurements are carried out only at the surface ocean (e.g. Bakker et al.2016), in part due to practical constraints of the “gold-standard” equilibrator-based methodology. Second, the concept of pCO2 has utility only in the context of air–sea CO2 exchange, which takes place only at the surface ocean.

However, recent developments in sensor technology are beginning to enable direct measurements of in situ pCO2 at depth in the ocean (Clarke et al.2017). There is also growing interest in calculating in situ pCO2 values at depth for intercomparison exercises in which the marine carbonate system has been overdetermined by measuring more than two of its core parameters (e.g. Raimondi et al.2019), and the relevant pressure correction is implemented in software tools such as seacarb and mocsy (Orr and Epitalon2015; Orr et al.2015). Therefore, we anticipate an increasing need for pressure-corrected pCO2 values, and while we have kept the approach in PyCO2SYS consistent with other CO2SYS software for now, we consider a robust implementation of these calculations to be an important target for future code development.

5.4 Computational speed

One does not choose to write code in Python for its computational speed. Therefore, while optimising performance was not ignored in developing PyCO2SYS, it was not a main focus. We compared the computational speed of PyCO2SYS against that of CO2SYS-MATLAB v3.2.0 across a few different tasks for reference purposes. We ran CO2SYS-MATLAB in both MATLAB itself (expensive, proprietary software) and GNU Octave, a free and open-source MATLAB clone.

Table 5Comparison of computational speed for various tasks with PyCO2SYS and CO2SYS-MATLAB running in both MATLAB and GNU Octave. Values shown are the mean ± standard deviation of seven runs. The tasks are described in Sect. 5.4.

Download Print Version | Download XLSX

The different tasks are described in the subsequent sections, and the results are summarised in Table 5. Details of the computer and software used for testing are provided in Appendix I.

Overall, the PyCO2SYS computation time has the same order of magnitude as CO2SYS-MATLAB, but it is generally somewhat slower. However, the difference is negligible in practice for relatively small datasets (up to about 105 data points) but may become more noticeable in larger calculations. Potential future improvements to PyCO2SYS's computational speed are discussed in Sect. 5.5.

5.4.1 All combinations

The “all combinations” task was the validation test described in Sect. 4.2.2: that is, a single call to the (Py)CO2SYS function that includes one calculation using every possible combination of parameter pair and optional setting (e.g. choices of parameterisations for the equilibrium constants) for a total of 40 800 data points. Both input and output conditions were computed.

CO2SYS-MATLAB completed this task in a very similar time in both MATLAB and GNU Octave, with the latter slightly faster, and PyCO2SYS took about 1.5 times longer (Table 5). However, this difference would generally be negligible, as all three implementations of the test had an average run time of less than 1 s.

5.4.2 GLODAP

In this task, (Py)CO2SYS was run across the entire GLODAPv2.2021 Merged Master File (Lauvset et al.2021) with AT and TC as the known parameter pair. This file contains a little over 1.3 million data points for each variable. The results in Table 5 show the mean and standard deviation of seven runs in each case.

This calculation is an example in which results would only be required under one set of temperature and pressure conditions rather than needing to evaluate both input and output conditions. This allows PyCO2SYS to be used more efficiently, as it only calculates output-condition results if they are explicitly requested (Sect. 2.3), whereas CO2SYS-MATLAB always calculates its results at both input and output conditions.

The results in Table 5 show that CO2SYS-MATLAB running in MATLAB was the fastest, with GNU Octave taking longer by a factor of about 1.3. When calculating only under input conditions, PyCO2SYS took longer by a factor of about 1.8 than CO2SYS-MATLAB running in MATLAB and by 3.8 if both the input- and output-condition calculations were carried out.

5.5 Outlook

The Autograd package that PyCO2SYS uses for automatic differentiation is still being maintained, and its most recent release (v1.3, July 2019) is stable, but it is no longer in active development. Its successor, JAX (Bradbury et al.2018), has further benefits including “just-in-time” code compilation and parallelisation. These features could speed up computation speed in PyCO2SYS, especially the components involving automatic differentiation, potentially by several orders of magnitude. However, JAX cannot currently run natively on the Microsoft Windows operating system, which would greatly restrict the usability of PyCO2SYS for the oceanographic research community. This limitation is due to JAX's dependence on the separate XLA (accelerated linear algebra) compiler rather than being an intrinsic issue with JAX itself. Should this compatibility issue be resolved in the future, we envision updating PyCO2SYS to use JAX instead of Autograd. This should be relatively straightforward thanks to the close similarities between the API (application programming interface) of these packages.

As future developments are made to PyCO2SYS, we will aim to maintain consistency with other CO2SYS-family tools but cannot guarantee that all new features or updates will be added simultaneously across all implementations. In practice, the workload required to achieve this is not currently feasible, and we would not wish to hold back development because of the time required to replicate changes across multiple implementations. That said, the results should remain consistent enough that users can select which implementation to use based on their preferred software environment rather than the other way around.

This ambition could also extend beyond the CO2SYS family of software. Independently developed tools for solving the marine carbonate system exist in other languages, such as seacarb in R (Gattuso et al.2021) and mocsy in Fortran (Orr and Epitalon2015). These give sufficiently consistent results with each other that the selection of which tool to use does not affect scientific interpretation (Orr et al.2015), and we have shown that PyCO2SYS is, and will remain, no exception. Even so, development and validation of PyCO2SYS so far have focused on comparisons with only CO2SYS-family software for practical reasons. Now that the basis of PyCO2SYS is established, we would welcome more direct interaction with the groups developing these other tools, working towards a set of marine–carbonate system-solving tools that return identical results regardless of the software platform. There can be a great advantage in having independent implementations led by different groups of researchers and developers. For example, this approach can help catch bugs and typographical errors, especially if each group extracts equations and parameterisations from the original literature instead of copying existing code. Working together, the groups would have a greater pool of knowledge and experience to identify errors in the literature (see e.g. Lewis and Wallace1998, their Appendix A), which are often unpublished and known only through personal communications. But calculations must be regularly compared with each other if this advantage is to be realised.

Thanks largely to the efforts of Orr et al. (2018), many tools now have an uncertainty propagation capability, as does PyCO2SYS. However, we still lack meaningful and statistically equivalent estimates for the actual uncertainties in the equilibrium constants. The software therefore stands ahead of our knowledge: as more work is done to robustly quantify these uncertainties, the tools are already in place to propagate them through to all marine carbonate system calculations.

As development of PyCO2SYS continues, we do not anticipate changing its fundamental approach to solving the marine carbonate system, but we will try to incorporate the latest research, including keeping up to date with new parameterisations, for example of stoichiometric equilibrium constants (e.g. Sulpis et al.2020; Schockman and Byrne2021). Integration with a speciation model that can determine the equilibrium constants based on chemical activities, rather than parameterising these based on salinity, is an area of interest (Turner et al.2016) but would likely require such substantial changes to constitute a separate software tool. We do envision further additions to the main chemical speciation function in PyCO2SYS, for example to better represent the impact of organic contributions to alkalinity (e.g. Cantrell et al.1990; Muller and Bleie2008; Kuliński et al.2014; Abril et al.2015; Ulfsbo et al.2015) – noting that a simplified representation of such extra components can already be modelled in PyCO2SYS (Sect. 3.3.1).

Through all these efforts, we aim to ensure that PyCO2SYS remains a reliable and comprehensive tool for analysing seawater chemistry from samples and experiments in the laboratory through to the changing marine carbonate system across the global ocean.

Appendix A: pH scales and conversions

The pH scales in PyCO2SYS are free (pHF), total (pHT), seawater (pHS), and NBS (pHN), defined following e.g. Zeebe and Wolf-Gladrow (2001) and Velo et al. (2010).


Here, γH+ is the chemical activity coefficient for H+ (Table 3). Note that in PyCO2SYS, [H+] in all these definitions is a substance content (Sect. 2.1). The pH values and stoichiometric equilibrium constants (K*) are thus converted between these different pH scales using the following factors.


Here, γH+ is the hydrogen ion activity calculated from temperature and salinity following either Peng et al. (1987) or Takahashi et al. (1982) (see Table 3). The different scales are denoted by the subscript and superscript letters, with F for free, T for total, S for seawater, and N for NBS. Convert from any pH scale A to any other pH scale B using these factors:

(A8) pH B = pH A + p Y A B = pH A - log 10 Y A B ,

or alternatively and equivalently

(A9) [ H + ] B = Y A B [ H + ] A .

The equations above are used in the same way to convert K* values between pH scales.

Appendix B: Total alkalinity and its components

Each equation here is written assuming that [H+] and all equilibrium constants (K*) are supplied on the same pH scale as each other.

B1 Total alkalinity

Total alkalinity (AT) is calculated as the sum of all its components (Dickson1981; Wolf-Gladrow et al.2007; Sharp and Byrne2020).

(B1) A T = A w + A C + A B + A P + A Si + A NH 3 + A H 2 S + A SO 4 + A F + A α + A β

Equations for all the individual alkalinity components (AC, AB, etc.) are given in the subsequent sections in terms of pH-independent total substance contents (TC, TB, etc.) and [H+].

B2 Water

(BR1) H 2 O OH - + H + ; K w * = [ OH - ] [ H + ]

(B2) A w = [ OH - ] - [ H + ] = K w * [ H + ] - [ H + ]

B3 Carbonic acid

(B3) T C = [ CO 2 ( aq ) ] + [ HCO 3 - ] + [ CO 3 2 - ]


AC can be expressed in terms of [H+] and any of TC, fCO2, [HCO3-], or [CO32-].


Undissociated H2CO3 is considered negligible and thus not explicitly modelled, but rather implicitly included as part of the [CO2(aq)] term (Zeebe and Wolf-Gladrow2001).

B4 Boric acid

(B9) T B = [ B ( OH ) 3 ] + [ B ( OH ) 4 - ]

(BR4) B ( OH ) 3 + H 2 O B ( OH ) 4 - + H + ; K B * = [ B ( OH ) 4 - ] [ H + ] [ B ( OH ) 3 ]

(B10) A B = [ B ( OH ) 4 - ] = T B K B * K B * + [ H + ]

B5 Phosphoric acid

(B11) T P = [ H 3 PO 4 ] + [ H 2 PO 4 - ] + [ HPO 4 2 - ] + [ PO 4 3 - ]


(B12) A P = [ HPO 4 2 - ] + 2 [ PO 4 3 - ] - [ H 3 PO 4 ] = T P ( K P 1 * K P 2 * [ H + ] + 2 K P 1 * K P 2 * K P 3 * - [ H + ] 3 ) K P 1 * K P 2 * K P 3 * + K P 1 * K P 2 * [ H + ] + K P 1 * [ H + ] 2 + [ H + ] 3

B6 Orthosilicic acid

(B13) T Si = [ H 4 SiO 4 ] + [ H 3 SiO 3 - ]

(BR8) H 4 SiO 4 H 3 SiO 4 - + H + ; K Si * = [ H 3 SiO 4 - ] [ H + ] [ H 4 SiO 4 ]

(B14) A Si = [ H 3 SiO 4 - ] = T Si K Si * K Si * + [ H + ]

Further deprotonation of H3SiO4- is considered negligible and thus not modelled.

B7 Ammonium

(B15) T NH 3 = [ NH 3 ] + [ NH 4 + ]

(BR9) NH 4 + NH 3 + H + ; K NH 3 * = [ NH 3 ] [ H + ] [ NH 4 + ]

(B16) A NH 3 = [ NH 3 ] = T NH 3 K NH 3 * K NH 3 * + [ H + ]

B8 Sulfide

(B17) T H 2 S = [ H 2 S ] + [ HS - ]

(BR10) H 2 S HS - + H + ; K H 2 S * = [ HS - ] [ H + ] [ H 2 S ]

(B18) A H 2 S = [ HS - ] = T H 2 S K H 2 S * K H 2 S * + [ H + ]

Further deprotonation of HS is considered negligible and thus not modelled (Schoonen and Barnes1988).

B9 Sulfate

(B19) T SO 4 = [ HSO 4 - ] + [ SO 4 2 - ]

(BR11) HSO 4 - SO 4 2 - + H + ; K SO 4 * = [ SO 4 2 - ] [ H + ] [ HSO 4 - ]

(B20) A SO 4 = - [ HSO 4 - ] = - T SO 4 1 + K SO 4 * / [ H + ]

Undissociated H2SO4 is considered negligible and thus not modelled.

B10 Fluoride

(B21) T F = [ HF ] + [ F - ]

(BR12) HF F - + H + ; K F * = [ F - ] [ H + ] [ HF ]

(B22) A F = - [ HF ] = - T F 1 + K F * / [ H + ]

B11 Arbitrary additional components

(B23) T α = [ H α ] + [ α - ]

(BR13) H α α - + H + ; K α * = [ α - ] [ H + ] [ H α ]

(B24) A α = - [ H α ] for  - log 10 ( K α * ) 4.5 + [ α - ] for  - log 10 ( K α * ) > 4.5

The reactions and equations for the second additional component β and its alkalinity contribution Aβ are identical to those given for α above. PyCO2SYS automatically determines how to modify the alkalinity equation following Eq. (B24) based on the user-provided Kα* and Kβ* values, with a zero level of protons corresponding to a pK* of 4.5 (Wolf-Gladrow et al.2007).

Though the definition of alkalinity (Dickson1981) states that species are separated into proton acceptors and donors based on their dissociation constant at zero ionic strength and 25 C, we use the user-defined dissociation constants at the given conditions because one cannot convert arbitrary dissociation constants to their alkalinity-relevant values. Interpretations of results when arbitrary components are supplied to PyCO2SYS with pK* values close to 4.5 should consider this nuance.

Appendix C: Solving the core marine carbonate system

Here, we lay out all the equations that are used to convert between different carbonate system parameters in PyCO2SYS. These follow long-established approaches from the literature (Zeebe and Wolf-Gladrow2001; Dickson et al.2007). The equations are organised based on which parameter pair is initially known.

C1 General considerations

C1.1 pH to [H+] conversions

As the stoichiometric equilibrium constants are converted to the user-specified pH scale, i.e. consistent with the pH values, pH and [H+] are interconverted in the equations throughout this section using

(C1) pH = - log 10 [ H + ]

regardless of which pH scale is being used.

C1.2 Known pCO2, xCO2 or [CO2(aq)]

If one of pCO2, xCO2, or [CO2(aq)] is in the known parameter pair, then its values are first converted to fCO2 as follows.

For known pCO2,

(C2) f CO 2 = G p CO 2 ,

where G is the fugacity factor (Table 2), typically near 0.997.

For known xCO2,

(C3) f CO 2 = G P v x CO 2 ,

where Pv is the humidity correction (Table 3):

(C4) P v = P a - p w ,

in which Pa is total atmospheric pressure (assumed to be 1 atm unless a different value is provided by the user) and pw is the water vapour pressure (Weiss and Price1980).

For known [CO2(aq)],

(C5) f CO 2 = [ CO 2 ( aq ) ] K 0 * ,

where K0* is the solubility factor for CO2 (Table 3).

The calculation steps given below for fCO2 are then followed to solve the core marine carbonate system. Afterwards, pCO2, xCO2, and [CO2(aq)] are calculated where they were not in the original known parameter pair: pCO2 and xCO2 are calculated using Eqs. (C2) and (C3), while [CO2(aq)] is calculated by difference using the definition of TC in Eq. (B3).

C2 Solving routines

C2.1 From AT and TC

An initial pH estimate is determined as described in Appendix F. The estimate is then revised using the iterative approach of Sect. 3.1, in which the AT(pHn,v) term in Eq. (2) is calculated from Eq. (B1) for AT substituting in Eq. (B5) for the AC term. Equation (2) is the automatically differentiated with respect to pH to obtain the ΔAT term in Eq. (1).

The components of TC are then calculated from TC and the final pH value.


C2.2 From AT and pH

First, we determine AC from known AT and pH by using Eq. (B1). TC is then calculated from AC.

(C9) T C = A C ( [ H + ] 2 + K 1 * [ H + ] + K 1 * K 2 * ) K 1 * ( [ H + ] + 2 K 2 * )

The components of TC are then calculated from TC and pH using Eqs. (C6), (C7), and (C8).

There is an upper limit on pH for each given AT value, above which negative AC would be required to balance Eq. (B1). PyCO2SYS prints a warning if such an impossible pairing is used and returns NaN (not a number) for TC (and all other results calculated from it) instead of a negative value.

C2.3 From AT and fCO2

An initial pH estimate is determined as described in Appendix F. The estimate is then revised using the iterative approach of Sect. 3.1, in which the AT(pHn,v) term in Eq. (2) is calculated from Eq. (B1) for AT substituting in Eq. (B6) for the AC term. Equation (2) is the automatically differentiated with respect to pH to obtain the ΔAT term in Eq. (1).

TC is then calculated from AT and pH following Sect. C2.2 and its remaining unknown components with Eqs. (C7) and (C8).

C2.4 From AT and [CO32-]

An initial pH estimate is determined as described in Appendix F. The estimate is then revised using the iterative approach of Sect. 3.1, in which the AT(pHn,v) term in Eq. (2) is calculated from Eq. (B1) for AT substituting in Eq. (B8) for the AC term. Equation (2) is the automatically differentiated with respect to pH to obtain the ΔAT term in Eq. (1). The lower of the two pH roots is returned by default, as discussed in Sect. 5.2.2.

TC is then calculated from AT and pH following Sect. C2.2 and its remaining unknown components with Eqs. (C6) and (C7).

C2.5 From AT and [HCO3-]

An initial pH estimate is determined as described in Appendix F. The estimate is then revised using the iterative approach of Sect. 3.1, in which the AT(pHn,v) term in Eq. (2) is calculated from Eq. (B1) for AT substituting in Eq. (B7) for the AC term. Equation (2) is the automatically differentiated with respect to pH to obtain the ΔAT term in Eq. (1).

TC is then calculated from AT and pH following Sect. C2.2 and its remaining unknown components with Eqs. (C6) and (C8).

C2.6 From TC and pH

First, AT is calculated from TC and pH using Eq. (B1). The components of TC are then calculated from TC and pH using Eqs. (C6), (C7), and (C8).

C2.7 From TC and fCO2

First, pH is calculated from TC and fCO2 using

(C10) [ H + ] = K 1 * r + ( K 1 * r ) 2 + 4 ( 1 - r ) K 1 * K 2 * r 2 ( 1 - r ) ,


(C11) r = K 0 * f CO 2 / T C .

AT and the remaining unknown components of TC are then calculated from TC and pH using Eqs. (B1), (C7), and (C8).

C2.8 From TC and [CO32-]

First, pH is calculated from TC and [CO32-] using

(C12) [ H + ] = - K 1 * + K 1 * 2 - 4 K 1 * K 2 * ( 1 - T C / [ CO 3 2 - ] ) 2 .

AT and the remaining unknown components of TC are then calculated from TC and pH using Eqs. (B1), (C6), and (C7).

C2.9 From TC and [HCO3-]

First, pH is calculated from TC and [HCO3-] using

(C13) [ H + ] = T C - [ HCO 3 - ] - ( [ HCO 3 - ] - T C ) 2 - 4 [ HCO 3 - ] 2 K 2 * / K 1 * 2 [ HCO 3 - ] / K 1 * .

AT and the remaining unknown components of TC are then calculated from TC and pH using Eqs. (B1), (C6), and (C8).

C2.10 From pH and fCO2

First, TC is calculated from pH and fCO2 using

(C14) T C = K 0 * f CO 2 [ H + ] 2 + K 1 * [ H + ] + K 1 * K 2 * [ H + ] 2 .

AT and the remaining unknown components of TC are then calculated from TC and pH using Eqs. (B1), (C7), and (C8).

C2.11 From pH and [CO32-]

First, fCO2 is calculated from pH and [CO32-] using

(C15) f CO 2 = [ CO 3 2 - ] [ H + ] 2 K 0 * K 1 * K 2 * .

TC is then calculated from pH and fCO2 using Eq. (C14). Finally, AT and [HCO3-] are calculated from TC and pH using Eqs. (B1) and (C7), respectively.

C2.12 From pH and [HCO3-]

First, TC is calculated from pH and [HCO3-] using

(C16) T C = [ HCO 3 - ] 1 + [ H + ] K 1 * + K 2 * [ H + ] .

AT and the remaining unknown components of TC are then calculated from TC and pH using Eqs. (B1), (C6), and (C8).

C2.13 From fCO2 and [CO32-]

First, pH is calculated from fCO2 and [CO32-] using

(C17) [ H + ] = K 0 * K 1 * K 2 * f CO 2 [ CO 3 2 - ] .

TC is then calculated from pH and fCO2 using Eq. (C14). Finally, AT and [HCO3-] are calculated from TC and pH using Eqs. (B1) and (C7), respectively.

C2.14 From fCO2 and [HCO3-]

First, [CO32-] is calculated from fCO2 and [HCO3-] using

(C18) [ CO 3 2 - ] = [ HCO 3 - ] 2 K 2 * K 0 * K 1 * f CO 2 .

The pH is then calculated from fCO2 and [CO32-] using Eq. (C17). Next, TC is calculated from pH and fCO2 using Eq. (C14). Finally, AT is calculated from TC and pH using Eq. (B1).

C2.15 From [CO32-] and [HCO3-]

First, fCO2 is calculated from [CO32-] and [HCO3-] using

(C19) f CO 2 = [ HCO 3 - ] 2 K 2 * K 0 * K 1 * [ CO 3 2 - ] .

The pH is then calculated from fCO2 and [CO32-] using Eq. (C17). Next, TC is calculated from pH and fCO2 using Eq. (C14). Finally, AT is calculated from TC and pH using Eq. (B1).

Appendix D: Other marine carbonate system variables

Calcite and aragonite saturation states (Ω) are calculated from the definition

(D1) Ω = [ Ca 2 + ] [ CO 3 2 - ] K sp * ,

where Ksp* is the solubility product, a function of salinity, temperature, and pressure that is different for each mineral (Table 2).

The “substrate : inhibitor ratio” of Bach (2015) is calculated from the bicarbonate and free hydrogen ion contents.

(D2) SIR = [ HCO 3 - ] [ H + ]

Note that in Eq. (D2), the [H+] term is always calculated on the free pH scale of Eq. (A1).

Appendix E: Buffer factors with automatic differentiation

E1 Buffer factors of Egleston et al. (2010)

To evaluate the buffer factors of Egleston et al. (2010) with automatic differentiation (AD), we first evaluated the following partial differentials (with the subscripted variable held constant):

  • (TC/pH)AT by AD of Eq. (C9) with respect to pH;

  • (AT/pH)TC by AD of Eq. (B1), substituting AC by Eq. (B5), with respect to pH;

  • (ln[CO2(aq)]/pH)TC by taking the natural log of the product of K0* and Eq. (C6), then AD with respect to pH;

  • (ln[CO2(aq)]/pH)AT by taking the natural log of the product of K0* and Eq. (C6), substituting TC by Eq. (9), then AD with respect to pH.

The buffer factors γTC, γAT, βTC, and βAT are thus defined (Egleston et al.2010) and calculated in PyCO2SYS.


Here, e is Euler's number (2.71828).

For the saturation-state buffers ωTC and ωAT we also evaluate

  • (lnΩ/[CO32-]) by AD of the natural log of Ω(aragonite), calculated with Eq. (D1), with respect to [CO32-] (note that this is the same value as for Ω(calcite), due to the logarithm and the fact that these terms differ by the constant ratio of their solubility products);

  • ([CO32-]/pH)TC by AD of Eq. (C8) with respect to pH;

  • ([CO32-]/pH)AT by AD of Eq. (C8), substituting TC by Eq. (C9), with respect to pH.

The buffer factors are then given by the following.


The approach taken here avoids AD evaluations over the iterative solvers because, while possible, that is computationally slower than over non-iterative functions.

E2 Revelle factor

The Revelle factor (RFBroecker et al.1979) is computed from TC and γTC, with the latter evaluated as described in Sect. E1, following Egleston et al. (2010).

(E7) R F = f CO 2 T C T C f CO 2 = T C γ T C

E3 Isocapnic quotient and ψ

To evaluate the isocapnic quotient (Q) of Humphreys et al. (2018), we first evaluate the derivatives:

  • (TC/pH)fCO2 by AD of Eq. (C14) with respect to pH;

  • (AT/pH)fCO2 by AD of Eq. (B1) using Eq. (B6) for the AC term with respect to pH.

The isocapnic quotient is defined and calculated in PyCO2SYS as follows.

(E8) Q = A T T C f CO 2 = A T pH f CO 2 T C pH f CO 2 - 1

Finally, the “released CO2 : precipitated carbonate ratio” (ψ) of Frankignoulle et al. (1994) is calculated following Humphreys et al. (2018).

(E9) ψ = 2 Q - 1
Appendix F: Initial pH estimate when solving from AT and TC

For clarity in the equations in this section, we abbreviate [H+] as h.

Following Munhoven (2013), carbonate–borate alkalinity (ACB) from Eq. (3) as a function of TC and h is

(F1) A CB ( h , T C ) = T C K 1 * ( h + 2 K 2 * ) h 2 + K 1 * h + K 1 * K 2 * + T B K B * h + K B * .

This can be rearranged into a third-order polynomial in h:

(F2) P T C ( h ) = h 3 + h 2 g 2 ( T C ) + h g 1 ( T C ) + g 0 ( T C ) = 0 ,

with the following.


The initial h value is determined by

(F6) h 0 ( T C ) = 10 - 3 for  A T 0 h min + - P T C ( h min ) g 2 2 - 3 g 1 for  A T > 0 10 - 10 for  A T 2 T C + T B ,

where hmin is defined in Eq. (10). Negative ACB is impossible because its equation contains only positive terms, so the equations above cannot be applied if AT is indeed negative. The default h0 of 10−3 mol kg−1, corresponding to a pH of 3, is therefore used for that case (e.g. after the alkalinity end point in an acidimetric titration). The maximum possible ACB is 2TC+TB, where TC is entirely CO32- and TB is entirely B(OH)4-. Where AT is actually higher than this limit of this simplified expression, we expect a high pH (given the dominance of CO32- within TC) and therefore use an initial estimate pH of 10. Otherwise, hmin in Eq. (F6) is found using Eq. (10). A default h0 of 10−7 mol kg−1 (pH 7) is used when g22-3g10 in Eq. (F6) (Munhoven2013).

Appendix G: Revelle factor calculation errors in older versions of CO2SYS-MATLAB

Older versions of CO2SYS-MATLAB, including v2.0.5 (Orr et al.2018) from which PyCO2SYS was originally converted, have minor errors in how the Revelle factor is evaluated. These have been corrected in PyCO2SYS (also in CO2SYS-MATLAB v3.2.0 and CO2SYS-Excel v3), leading to small differences in the calculated values. These differences are on the order of 0.1; for context, the Revelle factor typically has a value on the order of 10. The differences are thus notable from a computational perspective (i.e. many orders of magnitude greater than solver tolerance and floating-point errors) but still mostly negligible in practical applications.

Rather than being corrected explicitly in PyCO2SYS, these errors are corrected automatically thanks to the approach of using automatic differentiation instead of finite-difference derivatives. The key errors in the original CO2SYS-MATLAB implementation of the finite-difference approach are the following.

  1. An incorrect reference TC value is used in the final evaluation. Rather than using the “central” TC value, the change in pCO2 is divided by the adjusted (TC−ΔTC).

  2. Under output conditions, the Peng correction is not included in the evaluation of the Revelle factor (Sect. 2.2).

The lower accuracy of the finite-difference method relative to automatic differentiation, particularly given the relatively large ΔTC used in the original finite-difference implementation (i.e. 1 µmol kg−1), explains the differences between the two approaches that remain after the errors above have been corrected.

Appendix H: Fixed Δa values for uncertainty analysis

Table H1Fixed Δa values for uncertainty analysis.

Download Print Version | Download XLSX

Appendix I: Set-up for computational speed testing

The computational speed tests described in Sect. 5.4 were run on an HP Spectre x360 laptop with an Intel Core i7-8565U CPU (1.80 GHz) and 16 GB of RAM. The operating system was Windows 10.

The Python tests were run using Python v3.9.7, Autograd v1.3, NumPy v1.21.2, and PyCO2SYS v1.8.0.

The MATLAB tests were run using MATLAB R2019b (Update 9) and CO2SYS-MATLAB v3.2.0.

The GNU Octave tests were run using GNU Octave v6.3.0 via its command-line interface and CO2SYS-MATLAB v3.2.0.

Code availability

The current version of PyCO2SYS is freely available from its GitHub repository at (last access: 23 December 2021) under the GNU General Public License v3. Installation is recommended from the Python Package Index (PyPI) via pip, and documentation is available online (, last access: 23 December 2021). The exact version of PyCO2SYS used to produce the results discussed in this paper (v1.8.0), including input data and scripts to run the model and perform all validation tests described here, is archived on Zenodo (, Humphreys et al.2021).

Data availability

No data sets were used in this article.

Author contributions

MPH was responsible for conceptualisation, methodology, software, validation, writing the original draft, and visualisation. ERL was responsible for software and writing (review and editing). JDS was responsible for software, validation, and writing (review and editing). DP was responsible for software and writing (review and editing).​​​​​​​

Competing interests

The contact author has declared that neither they nor their co-authors have any competing interests.


Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.


We thank Doug Wallace for providing useful comments on this paper, and we acknowledge his important role in the creation of the original CO2SYS software. We further acknowledge the developers of all subsequent versions of CO2SYS upon whose work PyCO2SYS was built. We thank Luke Gregor, Daniel Sandborn, and Abigail Schiller for code contributions including extending the range of data types with which PyCO2SYS can be used. We are grateful to Guy Munhoven and James Orr for their detailed and constructive reviews.

Review statement

This paper was edited by Paul Halloran and reviewed by James Orr and Guy Munhoven.


Abril, G., Bouillon, S., Darchambeau, F., Teodoru, C. R., Marwick, T. R., Tamooh, F., Ochieng Omengo, F., Geeraert, N., Deirmendjian, L., Polsenaere, P., and Borges, A. V.: Technical Note: Large overestimation of pCO2 calculated from pH and alkalinity in acidic, organic-rich freshwaters, Biogeosciences, 12, 67–78,, 2015. a

Álvarez, M., Sanleón-Bartolomé, H., Tanhua, T., Mintrop, L., Luchetta, A., Cantoni, C., Schroeder, K., and Civitarese, G.: The CO2 system in the Mediterranean Sea: a basin wide perspective, Ocean Sci., 10, 69–92,, 2014. a

Bach, L. T.: Reconsidering the role of carbonate ion concentration in calcification by marine organisms, Biogeosciences, 12, 4939–4951,, 2015. a, b

Bakker, D. C. E., Pfeil, B., Landa, C. S., Metzl, N., O'Brien, K. M., Olsen, A., Smith, K., Cosca, C., Harasawa, S., Jones, S. D., Nakaoka, S., Nojiri, Y., Schuster, U., Steinhoff, T., Sweeney, C., Takahashi, T., Tilbrook, B., Wada, C., Wanninkhof, R., Alin, S. R., Balestrini, C. F., Barbero, L., Bates, N. R., Bianchi, A. A., Bonou, F., Boutin, J., Bozec, Y., Burger, E. F., Cai, W.-J., Castle, R. D., Chen, L., Chierici, M., Currie, K., Evans, W., Featherstone, C., Feely, R. A., Fransson, A., Goyet, C., Greenwood, N., Gregor, L., Hankin, S., Hardman-Mountford, N. J., Harlay, J., Hauck, J., Hoppema, M., Humphreys, M. P., Hunt, C. W., Huss, B., Ibánhez, J. S. P., Johannessen, T., Keeling, R., Kitidis, V., Körtzinger, A., Kozyr, A., Krasakopoulou, E., Kuwata, A., Landschützer, P., Lauvset, S. K., Lefèvre, N., Lo Monaco, C., Manke, A., Mathis, J. T., Merlivat, L., Millero, F. J., Monteiro, P. M. S., Munro, D. R., Murata, A., Newberger, T., Omar, A. M., Ono, T., Paterson, K., Pearce, D., Pierrot, D., Robbins, L. L., Saito, S., Salisbury, J., Schlitzer, R., Schneider, B., Schweitzer, R., Sieger, R., Skjelvan, I., Sullivan, K. F., Sutherland, S. C., Sutton, A. J., Tadokoro, K., Telszewski, M., Tuma, M., van Heuven, S. M. A. C., Vandemark, D., Ward, B., Watson, A. J., and Xu, S.: A multi-decade record of high-quality fCO2 data in version 3 of the Surface Ocean CO2 Atlas (SOCAT), Earth Syst. Sci. Data, 8, 383–413,, 2016. a

Bradbury, J., Frostig, R., Hawkins, P., Johnson, M. J., Leary, C., Maclaurin, D., and Wanderman-Milne, S.: JAX: composable transformations of Python+NumPy programs, GitHub, available at: (last access: 23 December 2021), 2018. a

Branson, O.: oscarbranson/cbsyst: beta, Zenodo [code],, 2018. a

Broecker, W. S., Takahashi, T., Simpson, H. J., and Peng, T.-H.: Fate of Fossil Fuel Carbon Dioxide and the Global Carbon Budget, Science, 206, 409–418,, 1979. a

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

Cai, W.-J., Huang, W.-J., Luther, G. W., Pierrot, D., Li, M., Testa, J., Xue, M., Joesoef, A., Mann, R., Brodeur, J., Xu, Y.-Y., Chen, B., Hussain, N., Waldbusser, G. G., Cornwell, J., and Kemp, W. M.: Redox reactions and weak buffering capacity lead to acidification in the Chesapeake Bay, Nat. Commun., 8, 369,, 2017. a

Cantrell, K. J., Serkiz, S. M., and Perdue, E. M.: Evaluation of acid neutralizing capacity data for solutions containing natural organic acids, Geochim. Cosmochim. Acta, 54, 1247–1254,, 1990. a

Clarke, J. S., Achterberg, E. P., Connelly, D. P., Schuster, U., and Mowlem, M.: Developments in marine pCO2 measurement technology; towards sustained in situ observations, Trends Anal. Chem., 88, 53–61,, 2017. a

Clegg, S. L. and Whitfield, M.: A chemical model of seawater including dissolved ammonia and the stoichiometric dissociation constant of ammonia in estuarine water and seawater from 2 to 40 C, Geochim. Cosmochim. Acta, 59, 2403–2421,, 1995. a

Culkin, F.: The major constituents of sea water, in: Chemical Oceanography, edited by: Riley, J. P. and Skirrow, G., Academic Press, London, UK, vol. 1, 121–161, 1965. a, b

Deffeyes, K. S.: Carbonate Equilibria: A Graphic and Algebraic Approach, Limnol. Oceanogr., 10, 412–426,, 1965. a

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

Dickson, A. G.: An exact definition of total alkalinity and a procedure for the estimation of alkalinity and total inorganic carbon from titration data, Deep-Sea Res. Pt. A, 28, 609–623,, 1981. a, b, c, d, e, f

Dickson, A. G.: Standard potential of the reaction: AgCl(s)+0.5H2(g)=Ag(s)+HCl(aq), and the standard acidity constant of the ion HSO4- in synthetic sea water from 273.15 to 318.15 K, J. Chem. Thermodyn., 22, 113–127,, 1990a. a

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

Dickson, A. G. and Riley, J. P.: The estimation of acid dissociation constants in sea-water media from potentiometric titrations with strong base. II. The dissociation of phosphoric acid, Mar. Chem., 7, 101–109,, 1979. a

Dickson, A. G., Sabine, C. L., and Christian, J. R. (Eds.): Guide to Best Practices for Ocean CO2 Measurements, PICES Special Publication 3, North Pacific Marine Science Organization, Sidney, BC, Canada, 2007. a

Dickson, A. G., Camões, M. F., Spitzer, P., Fisicaro, P., Stoica, D., Pawlowicz, R., and Feistel, R.: Metrological challenges for measurements of key climatological observables. Part 3: seawater pH, Metrologia, 53, R26–R39,, 2015. a

Doney, S. C., Fabry, V. J., Feely, R. A., and Kleypas, J. A.: Ocean Acidification: The Other CO2 Problem, Annu. Rev. Marine Sci., 1, 169–192,, 2009. a

Edmond, J. M. and Gieskes, J. M. T. M.: On the calculation of the degree of saturation of sea water with respect to calcium carbonate under in situ conditions, Geochim. Cosmochim. Acta, 34, 1261–1291,, 1970. a

Egleston, E. S., Sabine, C. L., and Morel, F. M. M.: Revelle revisited: Buffer factors that quantify the response of ocean chemistry to changes in DIC and alkalinity, Global Biogeochem. Cy., 24, GB1002,, 2010. a, b, c, d, e, f, g, h, i, j

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

Frankignoulle, M., Canon, C., and Gattuso, J.-P.: Marine calcification as a source of carbon dioxide: Positive feedback of increasing atmospheric CO2, Limnol. Oceanogr., 39, 458–462,, 1994. a, b, c

Friedlingstein, P., O'Sullivan, M., Jones, M. W., Andrew, R. M., Hauck, J., Olsen, A., Peters, G. P., Peters, W., Pongratz, J., Sitch, S., Le Quéré, C., Canadell, J. G., Ciais, P., Jackson, R. B., Alin, S., Aragão, L. E. O. C., Arneth, A., Arora, V., Bates, N. R., Becker, M., Benoit-Cattin, A., Bittig, H. C., Bopp, L., Bultan, S., Chandra, N., Chevallier, F., Chini, L. P., Evans, W., Florentie, L., Forster, P. M., Gasser, T., Gehlen, M., Gilfillan, D., Gkritzalis, T., Gregor, L., Gruber, N., Harris, I., Hartung, K., Haverd, V., Houghton, R. A., Ilyina, T., Jain, A. K., Joetzjer, E., Kadono, K., Kato, E., Kitidis, V., Korsbakken, J. I., Landschützer, P., Lefèvre, N., Lenton, A., Lienert, S., Liu, Z., Lombardozzi, D., Marland, G., Metzl, N., Munro, D. R., Nabel, J. E. M. S., Nakaoka, S.-I., Niwa, Y., O'Brien, K., Ono, T., Palmer, P. I., Pierrot, D., Poulter, B., Resplandy, L., Robertson, E., Rödenbeck, C., Schwinger, J., Séférian, R., Skjelvan, I., Smith, A. J. P., Sutton, A. J., Tanhua, T., Tans, P. P., Tian, H., Tilbrook, B., van der Werf, G., Vuichard, N., Walker, A. P., Wanninkhof, R., Watson, A. J., Willis, D., Wiltshire, A. J., Yuan, W., Yue, X., and Zaehle, S.: Global Carbon Budget 2020, Earth Syst. Sci. Data, 12, 3269–3340,, 2020. a

Gattuso, J.-P., Epitalon, J.-M., Lavigne, H., Orr, J., Gentili, B., Hagens, M., Hofmann, A., Mueller, J.-D., Proye, A., Rae, J., and Soetaert, K.: seacarb: Seawater Carbonate Chemistry, available at:, last access: 23 December 2021. a

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

Hagens, M. and Middelburg, J. J.: Generalised expressions for the response of pH to changes in ocean chemistry, Geochim. Cosmochim. Acta, 187, 334–349,, 2016. a

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 Chem. Scand., 27, 931–944,, 1973a. a, b

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

Harris, C. R., Millman, K. J., van der Walt, S. J., Gommers, R., Virtanen, P., Cournapeau, D., Wieser, E., Taylor, J., Berg, S., Smith, N. J., Kern, R., Picus, M., Hoyer, S., van Kerkwijk, M. H., Brett, M., Haldane, A., del Río, J. F., Wiebe, M., Peterson, P., Gérard-Marchant, P., Sheppard, K., Reddy, T., Weckesser, W., Abbasi, H., Gohlke, C., and Oliphant, T. E.: Array programming with NumPy, Nature, 585, 357–362,, 2020. a

Humphreys, M. P., Daniels, C. J., Wolf-Gladrow, D. A., Tyrrell, T., and Achterberg, E. P.: On the influence of marine biogeochemical processes over CO2 exchange between the atmosphere and ocean, Mar. Chem., 199, 1–11,, 2018. a, b, c, d, e

Humphreys, M. P., Gregor, L., Pierrot, D., van Heuven, S. M. A. C., Lewis, E. R., and Wallace, D. W. R.: PyCO2SYS: marine carbonate system calculations in Python (v1.3.0), Zenodo [code],, 2020. a

Humphreys, M. P., Schiller, A. J., Sandborn, D., Gregor, L., Pierrot, D., van Heuven, S. M. A. C., Lewis, E. R., and Wallace, D. W. R.: PyCO2SYS: marine carbonate system calculations in Python (v1.8.0), Zenodo [code],, 2021. a, b, c, d

Ingle, S. E.: Solubility of calcite in the ocean, Mar. Chem., 3, 301–319,, 1975. a, b, c

Ingle, S. E., Culberson, C. H., Hawley, J. E., and Pytkowicz, R. M.: The solubility of calcite in seawater at atmospheric pressure and 35 ‰ salinity, Mar. Chem., 1, 295–307,, 1973. a

IUPAC: Compendium of Chemical Terminology, 2nd edn. (the “Gold Book”), Blackwell Scientific Publications, Oxford, UK,, 1997. a

JCGM: JCGM 100:2008 Evaluation of measurement data – Guide to the expression of uncertainty in measurement, Bureau International des Poids et Mesures, Sèvres, France, available at: (last access: 23 December 2021), 2008. a

Kester, D. R. and Pytkowicz, R. M.: Determination of the Apparent Dissociation Constants of Phosphoric Acid in Seawater, Limnol. Oceanogr., 12, 243–252,, 1967. a

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

Kuliński, K., Schneider, B., Hammer, K., Machulik, U., and Schulz-Bull, D.: The influence of dissolved organic matter on the acid–base system of the Baltic Sea, J. Mar. Syst., 132, 106–115,, 2014. a

Lauvset, S. K., Lange, N., Tanhua, T., Bittig, H. C., Olsen, A., Kozyr, A., Álvarez, M., Becker, S., Brown, P. J., Carter, B. R., Cotrim da Cunha, L., Feely, R. A., van Heuven, S., Hoppema, M., Ishii, M., Jeansson, E., Jutterström, S., Jones, S. D., Karlsen, M. K., Lo Monaco, C., Michaelis, P., Murata, A., Pérez, F. F., Pfeil, B., Schirnick, C., Steinfeldt, R., Suzuki, T., Tilbrook, B., Velo, A., Wanninkhof, R., Woosley, R. J., and Key, R. M.: An updated version of the global interior ocean biogeochemical data product, GLODAPv2.2021, Earth Syst. Sci. Data, 13, 5565–5589,, 2021. a

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. Acta, 74, 1801–1811,, 2010. a

Lewis, E. and Wallace, D. W. R.: Program Developed for CO2 System Calculations, ORNL/CDIAC-105, Carbon Dioxide Information Analysis Center, Oak Ridge National Laboratory, U.S. Department of Energy, Oak Ridge, TN, USA,, 1998. a, b, c, d, e, f

Li, Y.-H., Takahashi, T., and Broecker, W. S.: Degree of saturation of CaCO3 in the oceans, J. Geophys. Res., 74, 5507–5525,, 1969. a

Lueker, T. J., Dickson, A. G., and Keeling, C. D.: Ocean pCO2 calculated from dissolved inorganic carbon, alkalinity, and equations for K1 and K2: validation based on laboratory measurements of CO2 in gas and seawater at equilibrium, Mar. Chem., 70, 105–119,, 2000. a

Maclaurin, D.: Autograd: Automatic Differentiation for Python, in: Modeling, Inference and Optimization with Composable Differentiable Procedures, PhD thesis, Harvard University, Cambridge, MA, USA, 41–57, available at: (last access: 23 December 2021), 2016. a

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. a, b, c, d

Millero, F. J.: The thermodynamics of the carbonate system in seawater, Geochim. Cosmochim. Acta, 43, 1651–1661,, 1979. a, b, c, d

Millero, F. J.: Influence of pressure on chemical processes in the sea, in: Chemical Oceanography, edited by: Riley, J. P. and Chester, R., Academic Press, 1983. a, b, c, d, e, f

Millero, F. J.: Thermodynamics of the carbon dioxide system in the oceans, Geochim. Cosmochim. Acta, 59, 661–677,, 1995. a, b, c, d, e, f, g, h, i, j

Millero, F. J.: The Carbonate System in Marine Environments, in: Chemical Processes in Marine Environments, edited by: Gianguzza, A., Pelizetti, E., and Sammartano, S., Environmental Science, Springer, Berlin, Heidelberg, 9–41,, 2000. a

Millero, F. J.: Carbonate constants for estuarine waters, Mar. Freshw. Res., 61, 139–142,, 2010. a

Millero, F. J., Pierrot, D., Lee, K., Wanninkhof, R., Feely, R., Sabine, C. L., Key, R. M., and Takahashi, T.: Dissociation constants for carbonic acid determined from field measurements, Deep-Sea Res. Pt. I, 49, 1705–1723,, 2002. a

Millero, F. J., Graham, T. B., Huang, F., Bustos-Serrano, H., and Pierrot, D.: Dissociation constants of carbonic acid in seawater as a function of salinity and temperature, Mar. Chem., 100, 80–94,, 2006. a

Millero, F. J., Feistel, R., Wright, D. G., and McDougall, T. J.: The composition of Standard Seawater and the definition of the Reference-Composition Salinity Scale, Deep-Sea Res. Pt. I, 55, 50–72,, 2008. a

Mojica Prieto, F. J. and Millero, F. J.: The values of pK1+ pK2 for the dissociation of carbonic acid in seawater, Geochim. Cosmochim. Acta, 66, 2529–2540,, 2002. a

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

Muller, F. L. L. and Bleie, B.: Estimating the organic acid contribution to coastal seawater alkalinity by potentiometric titrations in a closed cell, Anal. Chim. Acta, 619, 183–191,, 2008. a

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,, 2013. a, b, c, d, e, f, g, h, i

Munhoven, G.: SolveSAPHE-r2 (v2.0.1): revisiting and extending the Solver Suite for Alkalinity-PH Equations for usage with CO2, HCO3 or CO32 input data, Geosci. Model Dev., 14, 4225–4240,, 2021. a, b, c

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

Orr, J. C. and Epitalon, J.-M.: Improved routines to model the ocean carbonate system: mocsy 2.0, Geosci. Model Dev., 8, 485–499,, 2015. a, b, c, d, e

Orr, J. C., Epitalon, J.-M., and Gattuso, J.-P.: Comparison of ten packages that compute ocean carbonate chemistry, Biogeosciences, 12, 1483–1510,, 2015. a, b, c, d, e

Orr, J. C., Najjar, R. G., Aumont, O., Bopp, L., Bullister, J. L., Danabasoglu, G., Doney, S. C., Dunne, J. P., Dutay, J.-C., Graven, H., Griffies, S. M., John, J. G., Joos, F., Levin, I., Lindsay, K., Matear, R. J., McKinley, G. A., Mouchet, A., Oschlies, A., Romanou, A., Schlitzer, R., Tagliabue, A., Tanhua, T., and Yool, A.: Biogeochemical protocols and diagnostics for the CMIP6 Ocean Model Intercomparison Project (OMIP), Geosci. Model Dev., 10, 2169–2199,, 2017. a

Orr, J. C., Epitalon, J.-M., Dickson, A. G., and Gattuso, J.-P.: Routine uncertainty propagation for the marine carbon dioxide system, Mar. Chem., 207, 84–107,, 2018. a, b, c, d, e, f, g, h, i, j, k, l, m, n

Park, P. K.: Oceanic CO2 system: an evaluation of ten methods of investigation, Limnol. Oceanogr., 14, 179–186,, 1969. a

Peng, T.-H., Takahashi, T., Broecker, W. S., and Olafsson, J.: Seasonal variability of carbon dioxide, nutrients and oxygen in the northern North Atlantic surface water: observations and a model, Tellus B, 39, 439–458,, 1987. a, b

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

Pierrot, D., Lewis, E., and Wallace, D. W. R.: MS Excel Program Developed for CO2 System Calculations, available at: (last access: 23 December 2021), 2006. a

Pierrot, D., Epitalon, J.-M., Orr, J. C., Lewis, E. R., and Wallace, D.: MS Excel program developed for CO2 system calculations – version 3.0, available at:, last access: 23 December 2021. a, b

Raimondi, L., Matthews, J. B. R., Atamanchuk, D., Azetsu-Scott, K., and Wallace, D. W. R.: The internal consistency of the marine carbon dioxide system for high latitude shipboard and in situ monitoring, Mar. Chem., 213, 49–70,, 2019. a

Revelle, R. and Suess, H. E.: Carbon Dioxide Exchange Between Atmosphere and Ocean and the Question of an Increase of Atmospheric CO2 during the Past Decades, Tellus, 9, 18–27,, 1957. a

Richier, S., Achterberg, E. P., Humphreys, M. P., Poulton, A. J., Suggett, D. J., Tyrrell, T., and Moore, C. M.: Geographical CO2 sensitivity of phytoplankton correlates with ocean buffer capacity, Glob. Change Biol., 24, 4438–4452,, 2018. a

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

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

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

Schockman, K. M. and Byrne, R. H.: Spectrophotometric Determination of the Bicarbonate Dissociation Constant in Seawater, Geochim. Cosmochim. Acta, 300, 231–245,, 2021. a, b, c

Schoonen, M. A. A. and Barnes, H. L.: An approximation of the second dissociation constant for H2S, Geochim. Cosmochim. Acta, 52, 649–654,, 1988. a

Sharp, J. D. and Byrne, R. H.: Carbonate ion concentrations in seawater: Spectrophotometric determination at ambient temperatures and evaluation of propagated calculation uncertainties, Mar. Chem., 209, 70–80,, 2019. a

Sharp, J. D. and Byrne, R. H.: Interpreting measurements of total alkalinity in marine and estuarine waters in the presence of proton-binding organic matter, Deep-Sea Res. Pt. I, 165, 103 338,, 2020. a, b

Sharp, J. D., Pierrot, D., Humphreys, M. P., Epitalon, J.-M., Orr, J. C., Lewis, E. R., and Wallace, D. W.: CO2SYSv3 for MATLAB, Zenodo [code],, 2020. a, b, c, d, e

Sillén, L. G., Martell, A. E., and Bjerrum, J.: Stability constants of metal-ion complexes, special publication, 17th edn., Chemical Society, London, UK, 1964. a

Sulpis, O., Lauvset, S. K., and Hagens, M.: Current estimates of K1* and K2* appear inconsistent with measured CO2 system parameters in cold oceanic regions, Ocean Sci., 16, 847–862,, 2020. a, b, c

Takahashi, T., Williams, R. T., and Bos, D. L.: Carbonate Chemistry, in: GEOSECS Pacific Expedition: Hydrographic Data, vol. 3, 77–105, National Science Foundation, Washington, D.C., 1982. a, b, c, d, e

Takahashi, T., Sutherland, S. C., Wanninkhof, R., Sweeney, C., Feely, R. A., Chipman, D. W., Hales, B., Friederich, G., Chavez, F., Sabine, C., Watson, A., Bakker, D. C., Schuster, U., Metzl, N., Yoshikawa-Inoue, H., Ishii, M., Midorikawa, T., Nojiri, Y., Körtzinger, A., Steinhoff, T., Hoppema, M., Olafsson, J., Arnarson, T. S., Tilbrook, B., Johannessen, T., Olsen, A., Bellerby, R., Wong, C., Delille, B., Bates, N., and de Baar, H. J.: Climatological mean and decadal change in surface ocean pCO2, and net sea–air CO2 flux over the global oceans, Deep-Sea Res. Pt. II, 56, 554–577,, 2009.  a

Turner, D. R., Achterberg, E. P., Chen, C.-T. A., Clegg, S. L., Hatje, V., Maldonado, M. T., Sander, S. G., Berg, V. D., G, C. M., and Wells, M.: Toward a Quality-Controlled and Accessible Pitzer Model for Seawater and Related Systems, Front. Mar. Sci., 3, 139,, 2016. a

Ulfsbo, A., Kuliński, K., Anderson, L. G., and Turner, D. R.: Modelling organic alkalinity in the Baltic Sea using a Humic-Pitzer approach, Mar. Chem., 168, 18–26,, 2015. a

Uppström, L. R.: The boron/chlorinity ratio of deep-sea water from the Pacific Ocean, Deep-Sea Res., 21, 161–162,, 1974. a

van Heuven, S., Pierrot, D., Rae, J. W. B., Lewis, E., and Wallace, D. W. R.: CO2SYS v 1.1, MATLAB program developed for CO2 system calculations, ORNL/CDIAC-105b, Carbon Dioxide Information Analysis Center, Oak Ridge National Laboratory, U.S. Department of Energy, Oak Ridge, TN, USA, available at: (last access: 23 December 2021), 2011. a, b, c

Velo, A., Pérez, F. F., Lin, X., Key, R. M., Tanhua, T., de la Paz, M., Olsen, A., van Heuven, S., Jutterström, S., and Ríos, A. F.: CARINA data synthesis project: pH data scale unification and cruise adjustments, Earth Syst. Sci. Data, 2, 133–155,, 2010. a, b

Waters, J., Millero, F. J., and Woosley, R. J.: Corrigendum to “The free proton concentration scale for seawater pH”, [MARCHE: 149 (2013) 8–22], Mar. Chem., 165, 66–67,, 2014. a, b

Waters, J. F. and Millero, F. J.: The free proton concentration scale for seawater pH, Mar. Chem., 149, 8–22,, 2013. a, b, c

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

Weiss, R. F. and Price, B. A.: Nitrous oxide solubility in water and seawater, Mar. Chem., 8, 347–359,, 1980. a, b

Wolf-Gladrow, D. A., Zeebe, R. E., Klaas, C., Körtzinger, A., and Dickson, A. G.: Total alkalinity: The explicit conservative expression and its application to biogeochemical processes, Mar. Chem., 106, 287–300,, 2007. a, b

Xu, Y.-Y., Pierrot, D., and Cai, W.-J.: Ocean carbonate system computation for anoxic waters using an updated CO2SYS program, Mar. Chem., 195, 90–93,, 2017. a, b

Yao, W. and Millero, F. J.: The chemistry of the anoxic waters in the Framvaren Fjord, Norway, Aquat. Geochem., 1, 53–88,, 1995. a, b, c

Zeebe, R. E. and Wolf-Gladrow, D.: CO2 in Seawater: Equilibrium, Kinetics, Isotopes, Elsevier Oceanography Series 65, Elsevier B.V., Amsterdam, the Netherlands, 2001. a, b, c, d, e, f, g, h

Short summary
The ocean helps to mitigate our impact on Earth's climate by absorbing about a quarter of the carbon dioxide (CO2) released by human activities each year. However, once absorbed, chemical reactions between CO2 and water reduce seawater pH (ocean acidification), which may have adverse effects on marine ecosystems. Our Python package, PyCO2SYS, models the chemical reactions of CO2 in seawater, allowing us to quantify the corresponding changes in pH and related chemical properties.