PyCO2SYS v1.8: marine carbonate system calculations in Python

Oceanic dissolved inorganic carbon (T_C) is the largest pool of carbon that substantially interacts with the atmosphere on human timescales. Oceanic T_C is increasing through uptake of anthropogenic carbon dioxide (CO₂), and seawater pH is decreasing as a consequence. Both the exchange of CO₂ 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 T_C, total alkalinity (A_T), pH, and seawater CO₂ fugacity ($f_{{\mathrm{CO}}_{\mathrm{2}}}$; or its partial pressure, $p_{{\mathrm{CO}}_{\mathrm{2}}}$, or its dry-air mole fraction, $x_{{\mathrm{CO}}_{\mathrm{2}}}$) – 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 (https://github.com/mvdh7/PyCO2SYS, last access: 23 December 2021) under the GNU General Public License v3, archived on Zenodo (Humphreys et al., 2021), and documented online (https://pyco2sys.readthedocs.io/en/latest/, last access: 23 December 2021).

1 Introduction

The ocean absorbs about a quarter of the anthropogenic carbon dioxide (CO₂) currently emitted each year (Friedlingstein et al., 2020). This absorption is a double-edged sword. Removing CO₂ from the atmosphere reduces the impact of these emissions on Earth's climate. However, CO₂ 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 CO₂ 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 (Millero, 2000). The core parameters are the substance contents of aqueous CO₂, the bicarbonate and carbonate ions formed by its hydration and dissociation (${\mathrm{HCO}}_{\mathrm{3}}^{-}$ and ${\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}$), and the sum of these three components (dissolved inorganic carbon, T_C); total alkalinity (A_T; Dickson, 1981); the fugacity, partial pressure, or dry-air mole fraction of CO₂ in seawater ($f_{{\mathrm{CO}}_{\mathrm{2}}}$, $p_{{\mathrm{CO}}_{\mathrm{2}}}$, or $x_{{\mathrm{CO}}_{\mathrm{2}}}$; Weiss, 1974); 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 (Park, 1969; Zeebe and Wolf-Gladrow, 2001).

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. Branson, 2018). 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 Wallace, 1998) 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 Byrne, 2019; 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⁻¹), 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 H₂O. Although sometimes referred to as molinity, the correct term is substance content (IUPAC, 1997), 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⁻¹) (Appendix A; Zeebe and Wolf-Gladrow, 2001; 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 1 Parameterisations 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).

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 2 Parameterisations 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: $K_{\mathrm{P}\mathrm{1}}^{*}$, $K_{\mathrm{P}\mathrm{2}}^{*}$, and $K_{\mathrm{P}\mathrm{3}}^{*}$ (Appendix B). ^d Copies the pressure correction for boric acid. ^e A constant value of 1 is used in this case, i.e. $p_{{\mathrm{CO}}_{\mathrm{2}}}=f_{{\mathrm{CO}}_{\mathrm{2}}}$.

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, $K_{\mathrm{1}}^{*}$ and $K_{\mathrm{2}}^{*}$, 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 3 Parameterisations 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.

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 Wallace, 1998). 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 $K_{{\mathrm{SO}}_{\mathrm{4}}}^{*}$, and (iii) the hydrogen fluoride dissociation constant $K_{\mathrm{HF}}^{*}$ (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 $K_{{\mathrm{SO}}_{\mathrm{4}}}^{*}$, $K_{\mathrm{HF}}^{*}$, $K_{\mathrm{sp}}^{*}$(calcite), $K_{\mathrm{sp}}^{*}$(aragonite), or $K_{\mathrm{0}}^{*}$. For $K_{{\mathrm{SO}}_{\mathrm{4}}}^{*}$ and $K_{\mathrm{HF}}^{*}$, 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 CO₂ solubility factor $K_{\mathrm{0}}^{*}$ (Weiss, 1974). This value and calculations of $f_{{\mathrm{CO}}_{\mathrm{2}}}$, $p_{{\mathrm{CO}}_{\mathrm{2}}}$, and $x_{{\mathrm{CO}}_{\mathrm{2}}}$ 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 A_T and T_C 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 A_T, T_C, pH (on any scale), $p_{{\mathrm{CO}}_{\mathrm{2}}}$, $f_{{\mathrm{CO}}_{\mathrm{2}}}$, $x_{{\mathrm{CO}}_{\mathrm{2}}}$, [CO₂(aq)], $\left[{\mathrm{HCO}}_{\mathrm{3}}^{-}\right]$, and $\left[{\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}\right]$. Any pair of these can be provided, except for two of $p_{{\mathrm{CO}}_{\mathrm{2}}}$, $f_{{\mathrm{CO}}_{\mathrm{2}}}$, $x_{{\mathrm{CO}}_{\mathrm{2}}}$, and [CO₂(aq)], which would not be valid because these are all directly proportional to each other at a given temperature, salinity, and atmospheric pressure.

Figure 1 Overview 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 [CO₂(aq)], $p_{{\mathrm{CO}}_{\mathrm{2}}}$, $f_{{\mathrm{CO}}_{\mathrm{2}}}$, or $x_{{\mathrm{CO}}_{\mathrm{2}}}$ 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 A_T and T_C, 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. A_T plus T_C or one of its components) are solved using a scheme that has been updated from previous versions of CO2SYS (Sect. 3.1). The A_T and T_C 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 A_T and T_C 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

$$\begin{array}{}\text{(1)}& {\mathrm{pH}}_{n+\mathrm{1}}={\mathrm{pH}}_n-{\displaystyle \frac{\mathrm{\Delta }{A}_{\mathrm{T}}({\mathrm{pH}}_n,v)}{\mathrm{\Delta }{A}_{\mathrm{T}}^{\prime }({\mathrm{pH}}_n,v)}},\end{array}$$

where $A_{\mathrm{T}}^{\prime }=\mathrm{d}{A}_{\mathrm{T}}/\mathrm{dpH}$ and

$$\begin{array}{}\text{(2)}& \mathrm{\Delta }{A}_{\mathrm{T}}({\mathrm{pH}}_n,v)=A_{\mathrm{T}}({\mathrm{pH}}_n,v)-A_{\mathrm{T}}\left(\mathrm{known}\right),\end{array}$$

in which v is any of T_C, $f_{{\mathrm{CO}}_{\mathrm{2}}}$, $\left[{\mathrm{HCO}}_{\mathrm{3}}^{-}\right]$, or $\left[{\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}\right]$. A_T(pH_n,v) is determined as described in Sects. C2.1 (when v is T_C) and C2.3–C2.5 (when v is one of $f_{{\mathrm{CO}}_{\mathrm{2}}}$, $\left[{\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}\right]$, or $\left[{\mathrm{HCO}}_{\mathrm{3}}^{-}\right]$).

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 (Maclaurin, 2016). 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⁻⁴ in CO2SYS-MATLAB, 10⁻⁸ 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 Wallace, 1998): “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 $K_{{\mathrm{SO}}_{\mathrm{4}}}^{*}$ and $K_{\mathrm{F}}^{*}$, which are always on the free scale (Sect. 2.2). Calculations of all alkalinity components except $\left[{\mathrm{HSO}}_{\mathrm{4}}^{-}\right]$ and [HF] have therefore always been correct. However, because $K_{{\mathrm{SO}}_{\mathrm{4}}}^{*}$ and $K_{\mathrm{F}}^{*}$ are always on the free scale, pH must be converted to this scale in order to determine the contributions of $\left[{\mathrm{HSO}}_{\mathrm{4}}^{-}\right]$ 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 [${\mathrm{HSO}}_{\mathrm{4}}^{-}$] and [HF] are each on the order of 10⁻¹⁰ µmol kg⁻¹ relative to A_T on the order of 2000 µmol kg⁻¹) 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 A_T. 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 A_T and T_C by considering only the contributions of carbonate and borate species to A_T, simplifying the A_T equation.

$$\begin{array}{}\text{(3)}& A_{\mathrm{CB}}=\left[{\mathrm{HCO}}_{\mathrm{3}}^{-}\right]+\mathrm{2}\left[{\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}\right]+\left[\mathrm{B}(\mathrm{OH}{)}_{\mathrm{4}}^{-}\right]\end{array}$$

Following Munhoven (2013) and as also implemented elsewhere (e.g. Orr and Epitalon, 2015), 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 T_C 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 (https://doi.org/10.5281/zenodo.3780139, Humphreys et al., 2020​​​​​​​) before the publication of that study.

3.2.1 Solving from A_T and $f_{{\mathrm{CO}}_{\mathrm{2}}}$

For clarity in the equations in this section, we abbreviate [CO₂(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 A_T plus any of $p_{{\mathrm{CO}}_{\mathrm{2}}}$, $x_{{\mathrm{CO}}_{\mathrm{2}}}$, or [CO₂(aq)].

First, $f_{{\mathrm{CO}}_{\mathrm{2}}}$ is converted to s using Eq. (C5). Carbonate–borate alkalinity (A_CB) as a function of s and h is

$$\begin{array}{}\text{(4)}& A_{\mathrm{CB}}(h,s)={\displaystyle \frac{K_{\mathrm{1}}^{*}s(h+\mathrm{2}{K}_{\mathrm{2}}^{*})}{h^{\mathrm{2}}}}+{\displaystyle \frac{K_{\mathrm{B}}^{*}{T}_{\mathrm{B}}}{h+K_{\mathrm{B}}^{*}}}.\end{array}$$

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

$$\begin{array}{}\text{(5)}& P_s(h,s)=h^{\mathrm{3}}+h^{\mathrm{2}}{g}_{\mathrm{2}}\left(s\right)+h{g}_{\mathrm{1}}\left(s\right)+g_{\mathrm{0}}\left(s\right)=\mathrm{0},\end{array}$$

with the following.

$$\begin{array}{}\text{(6)}& {\displaystyle }{g}_{\mathrm{2}}\left(s\right)=K_{\mathrm{B}}^{*}\left(\mathrm{1}-{\displaystyle \frac{T_{\mathrm{B}}}{A_{\mathrm{CB}}}}\right)-{\displaystyle \frac{K_{\mathrm{1}}^{*}s}{A_{\mathrm{CB}}}}\text{(7)}& {\displaystyle }{g}_{\mathrm{1}}\left(s\right)={\displaystyle \frac{(\mathrm{2}{K}_{\mathrm{2}}^{*}+K_{\mathrm{B}}^{*})K_{\mathrm{1}}^{*}s}{-A_{\mathrm{CB}}}}\text{(8)}& {\displaystyle }{g}_{\mathrm{0}}\left(s\right)={\displaystyle \frac{\mathrm{2}{K}_{\mathrm{1}}^{*}{K}_{\mathrm{2}}^{*}{K}_{\mathrm{B}}^{*}s}{-A_{\mathrm{CB}}}}\end{array}$$

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

$$\begin{array}{}\text{(9)}& h_{\mathrm{0}}\left(s\right)=\left\{\begin{array}{ll}{\mathrm{10}}^{-\mathrm{3}}& \text{for }{A}_{\mathrm{T}}\le \mathrm{0}\\ h_{\min }+\sqrt{-\frac{P_s\left(h_{\min }\right)}{\sqrt{g_{\mathrm{2}}^{\mathrm{2}}-\mathrm{3}{g}_{\mathrm{1}}}}}& \text{for }{A}_{\mathrm{T}}>\mathrm{0}.\end{array}\right.\end{array}$$

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

$$\begin{array}{}\text{(10)}& h_{\min }=\left\{\begin{array}{ll}\left(-g_{\mathrm{2}}+\sqrt{g_{\mathrm{2}}^{\mathrm{2}}-\mathrm{3}{g}_{\mathrm{1}}}\right)/\mathrm{3}& \text{for }{g}_{\mathrm{2}}<\mathrm{0}\\ -g_{\mathrm{1}}/\left(g_{\mathrm{2}}+\sqrt{g_{\mathrm{2}}^{\mathrm{2}}-\mathrm{3}{g}_{\mathrm{1}}}\right)& \text{for }{g}_{\mathrm{2}}\ge \mathrm{0}\end{array}\right.\end{array}$$

When A_T is positive, the square-rooted term $g_{\mathrm{2}}^{\mathrm{2}}-\mathrm{3}{g}_{\mathrm{1}}$ is always greater than zero, and thus h_min has a real value. However, there is an additional constraint: A_CB cannot be greater than 2T_C+T_B (Munhoven, 2013). If A_T is actually greater than this limit, then we use a default h₀ of 10⁻⁷ mol kg⁻¹ instead (pH 7).

3.2.2 Solving from A_T and [${\mathrm{HCO}}_{\mathrm{3}}^{-}$]

For clarity in the equations in this section, we abbreviate [${\mathrm{HCO}}_{\mathrm{3}}^{-}$] as b and [H⁺] as h.

Carbonate–borate alkalinity as a function of b is

$$\begin{array}{}\text{(11)}& A_{\mathrm{CB}}(h,b)=b+{\displaystyle \frac{\mathrm{2}{K}_{\mathrm{2}}^{*}b}{h}}+{\displaystyle \frac{K_{\mathrm{B}}^{*}{T}_{\mathrm{B}}}{h+K_{\mathrm{B}}^{*}}}.\end{array}$$

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

$$\begin{array}{}\text{(12)}& P_b(h,b)=h^{\mathrm{2}}{g}_{\mathrm{2}}\left(b\right)+h{g}_{\mathrm{1}}\left(b\right)+g_{\mathrm{1}}\left(b\right)=\mathrm{0},\end{array}$$

with the following.

$$\begin{array}{}\text{(13)}& {\displaystyle }{g}_{\mathrm{2}}\left(b\right)=b-A_{\mathrm{CB}}\text{(14)}& {\displaystyle }{g}_{\mathrm{1}}\left(b\right)=K_{\mathrm{B}}^{*}(b+T_{\mathrm{B}}-A_{\mathrm{CB}})+\mathrm{2}{K}_{\mathrm{2}}^{*}b\text{(15)}& {\displaystyle }{g}_{\mathrm{0}}\left(b\right)=\mathrm{2}{K}_{\mathrm{2}}^{*}{K}_{\mathrm{B}}^{*}b\end{array}$$

The initial h value is estimated following

$$\begin{array}{}\text{(16)}& h_{\mathrm{0}}\left(b\right)=\left\{\begin{array}{ll}\frac{-g_{\mathrm{1}}-\sqrt{g_{\mathrm{1}}^{\mathrm{2}}-\mathrm{4}{g}_{\mathrm{0}}{g}_{\mathrm{2}}}}{\mathrm{2}{g}_{\mathrm{2}}}& \text{for }b<A_{\mathrm{T}}\\ {\mathrm{10}}^{-\mathrm{3}}& \text{for }b\ge A_{\mathrm{T}}.\end{array}\right.\end{array}$$

When b<A_T, the square-rooted term $g_{\mathrm{1}}^{\mathrm{2}}-\mathrm{4}{g}_{\mathrm{0}}{g}_{\mathrm{2}}$ is always positive, and thus h₀(b) has a real value. Otherwise, b can only be greater than A_T if the negative components of A_T 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 A_T and [${\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}$]

For clarity in the equations in this section, we abbreviate [${\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}$] here as c and [H⁺] as h.

Carbonate–borate alkalinity as a function of c is

$$\begin{array}{}\text{(17)}& A_{\mathrm{CB}}(h,c)={\displaystyle \frac{ch}{K_{\mathrm{2}}^{*}}}+\mathrm{2}c+{\displaystyle \frac{K_{\mathrm{B}}^{*}{T}_{\mathrm{B}}}{h+K_{\mathrm{B}}^{*}}}.\end{array}$$

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

$$\begin{array}{}\text{(18)}& P_c(h,c)=h^{\mathrm{2}}{g}_{\mathrm{2}}\left(c\right)+h{g}_{\mathrm{1}}\left(c\right)+g_{\mathrm{0}}\left(c\right)=\mathrm{0},\end{array}$$

with the following.

$$\begin{array}{}\text{(19)}& {\displaystyle }{g}_{\mathrm{2}}\left(c\right)=c\text{(20)}& {\displaystyle }{g}_{\mathrm{1}}\left(c\right)=K_{\mathrm{B}}^{*}c+K_{\mathrm{2}}^{*}(\mathrm{2}c-A_{\mathrm{CB}})\text{(21)}& {\displaystyle }{g}_{\mathrm{0}}\left(c\right)=K_{\mathrm{2}}^{*}{K}_{\mathrm{B}}^{*}(\mathrm{2}c+T_{\mathrm{B}}-A_{\mathrm{CB}})\end{array}$$

The initial h value is estimated following

$$\begin{array}{}\text{(22)}& h_{\mathrm{0}}\left(c\right)=\left\{\begin{array}{ll}\frac{-g_{\mathrm{1}}+\sqrt{g_{\mathrm{1}}^{\mathrm{2}}-\mathrm{4}{g}_{\mathrm{0}}{g}_{\mathrm{2}}}}{\mathrm{2}{g}_{\mathrm{2}}}& \text{for }{A}_{\mathrm{T}}>\mathrm{2}c+T_{\mathrm{B}}\\ {\mathrm{10}}^{-\mathrm{3}}& \text{for }{A}_{\mathrm{T}}\le \mathrm{2}c+T_{\mathrm{B}}.\end{array}\right.\end{array}$$

When $\mathrm{2}c+T_{\mathrm{B}}<A_{\mathrm{T}}$, the square-rooted term $g_{\mathrm{1}}^{\mathrm{2}}-\mathrm{4}{g}_{\mathrm{0}}{g}_{\mathrm{2}}$ is always positive, and thus h₀(c) has a real value. Otherwise, 2c+T_B can only be greater than A_T if the negative components of A_T 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⁻¹ K⁻¹. PyCO2SYS uses the 2018 CODATA-recommended value by default instead (i.e. 8.314462618 J mol⁻¹ K⁻¹), consistent with CO2SYS-MATLAB v3.2.0 and CO2SYS-Excel v3. This has a minor effect on conversions between $p_{{\mathrm{CO}}_{\mathrm{2}}}$, $f_{{\mathrm{CO}}_{\mathrm{2}}}$, and $x_{{\mathrm{CO}}_{\mathrm{2}}}$ (less than 10⁻⁴ %), as well as on the pressure corrections for the equilibrium constants (less than 10⁻³ % 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. ${\mathrm{HCO}}_{\mathrm{3}}^{-}$) 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 $p_{{\mathrm{CO}}_{\mathrm{2}}}$ corresponding to a fractional change in T_C at constant A_T (Revelle and Suess, 1957). Frankignoulle (1994) derived a broader set of buffer factors for the marine carbonate system, quantifying the responses of several different parameters to changes in T_C and A_T; 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 T_C required to return to the original seawater $p_{{\mathrm{CO}}_{\mathrm{2}}}$ after the action of calcification (which reduces A_T and T_C in a 2 : 1 ratio) or CaCO₃ dissolution (the reverse). Humphreys et al. (2018) introduced the “isocapnic quotient” (Q), which is the ratio of A_T to T_C change that does not affect seawater $p_{{\mathrm{CO}}_{\mathrm{2}}}$, thus generalising the concept of ψ for application to all biogeochemical processes that affect A_T and T_C (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. Orr, 2011; Á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 ${\mathit{\gamma }}_{T_{\mathrm{C}}}$ 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 $p_{{\mathrm{CO}}_{\mathrm{2}}}$, $f_{{\mathrm{CO}}_{\mathrm{2}}}$, and $x_{{\mathrm{CO}}_{\mathrm{2}}}$, 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 $p_{{\mathrm{CO}}_{\mathrm{2}}}$ and $f_{{\mathrm{CO}}_{\mathrm{2}}}$ 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 $p_{{\mathrm{CO}}_{\mathrm{2}}}$, $f_{{\mathrm{CO}}_{\mathrm{2}}}$, and $x_{{\mathrm{CO}}_{\mathrm{2}}}$ 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, $p_{{\mathrm{CO}}_{\mathrm{2}}}$, $f_{{\mathrm{CO}}_{\mathrm{2}}}$, $x_{{\mathrm{CO}}_{\mathrm{2}}}$, and [CO₂(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 $p_{{\mathrm{CO}}_{\mathrm{2}}}$, $f_{{\mathrm{CO}}_{\mathrm{2}}}$, $x_{{\mathrm{CO}}_{\mathrm{2}}}$, and [CO₂(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 $p_{{\mathrm{CO}}_{\mathrm{2}}}$ 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 $p_{{\mathrm{CO}}_{\mathrm{2}}}$ 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 2 Schematic 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 T_C 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 Wallace, 1998) 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.

$$\begin{array}{}\text{(23)}& {\displaystyle \frac{\partial r\left(a\right)}{\partial a}}\approx {\displaystyle \frac{r(a+\mathrm{\Delta }a)-r\left(a\right)}{\mathrm{\Delta }a}}\end{array}$$

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 3 An example figure used to select a suitable Δa value for uncertainty propagation, in this case for $\mathrm{\Delta }{K}_{\mathrm{1}}^{*}$. The pH was calculated from A_T and T_C, and then the value of $K_{\mathrm{1}}^{*}$ 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⁻¹²), 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

$$\begin{array}{}\text{(24)}& {\mathit{\sigma }}^{\mathrm{2}}\left(r\right)=\sum _i{\left({\displaystyle \frac{\partial r}{\partial a_i}}\right)}^{\mathrm{2}}{\mathit{\sigma }}^{\mathrm{2}}\left(a_i\right),\end{array}$$

where σ is the uncertainty as a standard deviation (thus, σ² 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 (JCGM, 2008).

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 (https://github.com/mvdh7/PyCO2SYS), 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 (https://doi.org/10.5281/zenodo.5602840, 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 4 The status badge for the validation tests, which are publicly visible at PyCO2SYS's GitHub repository (https://github.com/mvdh7/PyCO2SYS#pyco2sys, 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⁻⁸ in pH). The maximum absolute difference in each parameter across all possible input pair combinations is acceptably small (Table 4).

Table 4 Results 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⁻¹, total phosphate 1 µmol kg⁻¹, total ammonia 2 µmol kg⁻¹, total sulfide 3 µmol kg⁻¹. The pH solver tolerance in PyCO2SYS is 10⁻⁸ in terms of pH.

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 T_B) 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⁻¹² % for the Egleston et al. (2010) buffers, 10⁻⁹ % for ψ and Q, and 10⁻⁷ % 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 $T_{{\mathrm{SO}}_{\mathrm{4}}}$ and T_F 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=10⁴) 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⁻¹⁰). 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⁻¹², relative tolerance 10⁻¹⁶, in pK^* values or in µmol kg⁻¹) 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 10⁵ 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⁻¹⁰ %, 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⁻⁴ instead of 10⁻⁸ (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⁻⁵ %: 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⁻² % 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⁻¹⁰ % and with one exception all less than 10⁻⁶ %. The exception, a difference still less than 10⁻³ %, is for the aqueous CO₂ 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 [CO₂(aq)]: by difference from known T_C, $\left[{\mathrm{HCO}}_{\mathrm{3}}^{-}\right]$, and $\left[{\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}\right]$; from any one of these three variables, [H⁺], and $K_{\mathrm{1}}^{*}$ and $K_{\mathrm{2}}^{*}$ equilibrium constants using the equations in Appendix C (Sects. C2.6, C2.11, and C2.12); or from $f_{{\mathrm{CO}}_{\mathrm{2}}}$ or $p_{{\mathrm{CO}}_{\mathrm{2}}}$ and the CO₂ solubility constant ($K_{\mathrm{0}}^{*}$). 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 [CO₂(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 T_C, $\left[{\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}\right]$, $\left[{\mathrm{HCO}}_{\mathrm{3}}^{-}\right]$, and [CO₂(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⁻³ % under the same input conditions and all the propagated uncertainties reported by Orr et al. (2018) in their Table 4 to within 10⁻⁴ %. We consider all these differences to be negligible.

Across all combinations of optional parameters, mean uncertainties in A_T, T_C, $p_{{\mathrm{CO}}_{\mathrm{2}}}$, $f_{{\mathrm{CO}}_{\mathrm{2}}}$, $\left[{\mathrm{HCO}}_{\mathrm{3}}^{-}\right]$, $\left[{\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}\right]$, [CO₂(aq)], Ω(calcite), Ω(aragonite), and $x_{{\mathrm{CO}}_{\mathrm{2}}}$ 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 A_T should change through the titration due to acid addition, accounting for dilution of A_T, T_C, and all other dissolved solutes by acid addition, and then solving the carbonate system for pH from the so-determined A_T and T_C. Being tested here is the ability to solve for pH from known A_T and T_C across a wide range of pH and A_T values, including negative A_T.

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⁻¹ 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 5 Initial estimates (solid lines) and final solutions (dashed lines) of pH from known parameter pairs of total alkalinity (2.3 mmol kg⁻¹) with a range of values for (a) dissolved inorganic carbon (T_C), (b) aqueous CO₂ 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⁻⁸ 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 A_T–$\left[{\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}\right]$ and T_C–$\left[{\mathrm{HCO}}_{\mathrm{3}}^{-}\right]$ pairs can correspond to two different pH values (Deffeyes, 1965; Zeebe and Wolf-Gladrow, 2001; Munhoven, 2021). 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 A_T–pH solvers can be thought of as working by evaluating A_T at a sequence of different possible pH values until the pH that returns the true A_T is found. This pH is known as the “root” of the A_T–pH equation. The difference between the true A_T 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 6 Residuals between known A_T (2.3 mmol kg⁻¹) and (i) carbonate–borate alkalinity (solid lines; A_CB) from Eqs. (F1), (4), (11), and (17), as well as (ii) total alkalinity (dashed lines; A_T) from Eq. (B1), calculated across a range of pH, with a second known parameter of (a) dissolved inorganic carbon (2.15 mmol kg⁻¹), (b) CO₂ fugacity (600 µatm), (c) bicarbonate ion content (2011 µmol kg⁻¹), and (d) carbonate ion content (116 µmol kg⁻¹), 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 A_CB and A_T equations, respectively (Sect. 3.1.1 and 3.2). The similarity between the A_CB and A_T 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 A_T–$\left[{\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}\right]$ 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 A_T–$\left[{\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}\right]$ parameter pair, as follows. The lower-pH root corresponds to typical seawater: a relatively high-T_C system, wherein bicarbonate ions are the main component of T_C, and carbonate alkalinity ($\left[{\mathrm{HCO}}_{\mathrm{3}}^{-}\right]+\mathrm{2}\left[{\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}\right]$) is the main component of A_T. The higher-pH root corresponds to a low-T_C system, wherein virtually all of T_C is in the form of carbonate ion, and A_T is dominated by non-carbonate species (Fig. 7).

Figure 7 Main components of (a) total alkalinity (A_T) and (b) dissolved inorganic carbon (T_C) at the two possible pH roots for a known parameter pair of A_T (2300 µmol kg⁻¹) and $\left[{\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}\right]$ (120 µmol kg⁻¹). The low-pH root (left) represents typical seawater, with relatively high T_C (2143 µmol kg⁻¹) and both A_T and T_C dominated by the bicarbonate ion (${\mathrm{HCO}}_{\mathrm{3}}^{-}$). The high-pH root (right) has the same A_T and $\left[{\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}\right]$, but A_T is dominated by hydroxide (OH⁻), and T_C is much lower (122 µmol kg⁻¹) and comprised almost entirely of ${\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}$. These calculations were carried out at 15 ^∘C, with a practical salinity of 35 and zero nutrients. If nutrients were present, then like borate ($\mathrm{B}(\mathrm{OH}{)}_{\mathrm{4}}^{-}$) they would have different contributions to A_T 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-Gladrow, 2001), there are also two possible pH solutions for the T_C–$\left[{\mathrm{HCO}}_{\mathrm{3}}^{-}\right]$ parameter pair. We conceptualise these roots as follows: the remaining portion of T_C not accounted for by ${\mathrm{HCO}}_{\mathrm{3}}^{-}$ is either dominantly composed of ${\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}$ if the solution's pH is closer to p$K_{\mathrm{2}}^{*}$ (i.e. higher) or of CO₂(aq) if the pH is closer to p$K_{\mathrm{1}}^{*}$ (i.e. lower).

Solving from T_C and $\left[{\mathrm{HCO}}_{\mathrm{3}}^{-}\right]$ is more straightforward than from A_T and $\left[{\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}\right]$ 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 8 Main components of (a) total alkalinity (A_T) and (b) dissolved inorganic carbon (T_C) at the two possible pH roots for a known parameter pair of T_C (2100 µmol kg⁻¹) and bicarbonate ion content ($\left[{\mathrm{HCO}}_{\mathrm{3}}^{-}\right]$; 1900 µmol kg⁻¹). The high-pH root (left) represents typical seawater, in which most of the T_C not accounted for by $\left[{\mathrm{HCO}}_{\mathrm{3}}^{-}\right]$ is composed of $\left[{\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}\right]$; A_T is relatively high (2364 µmol kg⁻¹) and $f_{{\mathrm{CO}}_{\mathrm{2}}}$ low (331 µatm). In the low-pH root (right), the non-${\mathrm{HCO}}_{\mathrm{3}}^{-}$ portion of T_C is instead dominated by CO₂(aq); alkalinity is lower (1932 µmol kg⁻¹) and $f_{{\mathrm{CO}}_{\mathrm{2}}}$ 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 $p_{{\mathrm{CO}}_{\mathrm{2}}}$

In PyCO2SYS, $p_{{\mathrm{CO}}_{\mathrm{2}}}$ (and by extension, $f_{{\mathrm{CO}}_{\mathrm{2}}}$ and $x_{{\mathrm{CO}}_{\mathrm{2}}}$) 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 $p_{{\mathrm{CO}}_{\mathrm{2}}}$ 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 Epitalon, 2015). 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 $p_{{\mathrm{CO}}_{\mathrm{2}}}$ (i.e. a pressure correction for $K_{\mathrm{0}}^{*}$ and the fugacity factor; Appendix C1.2) is theoretically possible (Weiss, 1974; Orr and Epitalon, 2015), it could be argued that this is unnecessary. First, the vast majority of $p_{{\mathrm{CO}}_{\mathrm{2}}}$ 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 $p_{{\mathrm{CO}}_{\mathrm{2}}}$ has utility only in the context of air–sea CO₂ exchange, which takes place only at the surface ocean.

However, recent developments in sensor technology are beginning to enable direct measurements of in situ $p_{{\mathrm{CO}}_{\mathrm{2}}}$ at depth in the ocean (Clarke et al., 2017). There is also growing interest in calculating in situ $p_{{\mathrm{CO}}_{\mathrm{2}}}$ 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 Epitalon, 2015; Orr et al., 2015). Therefore, we anticipate an increasing need for pressure-corrected $p_{{\mathrm{CO}}_{\mathrm{2}}}$ 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 5 Comparison 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.

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 10⁵ 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 A_T and T_C 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 Epitalon, 2015). 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 Wallace, 1998, 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 Byrne, 2021). 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 Bleie, 2008; 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 (pH_F), total (pH_T), seawater (pH_S), and NBS (pH_N), defined following e.g. Zeebe and Wolf-Gladrow (2001) and Velo et al. (2010).

$$\begin{array}{}\text{(A1)}& {\displaystyle }{\mathrm{pH}}_{\mathrm{F}}=-{\log }_{\mathrm{10}}\mathit{\left\{}\right[{\mathrm{H}}^{+}\left]\mathit{\right\}}\text{(A2)}& {\displaystyle }{\mathrm{pH}}_{\mathrm{T}}=-{\log }_{\mathrm{10}}\mathit{\left\{}\right[{\mathrm{H}}^{+}\left]\right(\mathrm{1}+T_{{\mathrm{SO}}_{\mathrm{4}}}/K_{{\mathrm{SO}}_{\mathrm{4}}}^{*}\left)\mathit{\right\}}\text{(A3)}& {\displaystyle }{\mathrm{pH}}_{\mathrm{S}}=-{\log }_{\mathrm{10}}\mathit{\left\{}\right[{\mathrm{H}}^{+}\left]\right(\mathrm{1}+T_{{\mathrm{SO}}_{\mathrm{4}}}/K_{{\mathrm{SO}}_{\mathrm{4}}}^{*}+T_{\mathrm{F}}/K_{\mathrm{F}}^{*}\left)\mathit{\right\}}\text{(A4)}& {\displaystyle }{\mathrm{pH}}_{\mathrm{N}}=-{\log }_{\mathrm{10}}\mathit{\left\{}\right[{\mathrm{H}}^{+}\left]\right(\mathrm{1}+T_{{\mathrm{SO}}_{\mathrm{4}}}/K_{{\mathrm{SO}}_{\mathrm{4}}}^{*}+T_{\mathrm{F}}/K_{\mathrm{F}}^{*}\left){\mathit{\gamma }}_{{\mathrm{H}}^{+}}\mathit{\right\}}\end{array}$$

Here, ${\mathit{\gamma }}_{{\mathrm{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.

$$\begin{array}{}\text{(A5)}& {\displaystyle }{Y}_{\mathrm{F}}^{\mathrm{T}}=\mathrm{1}+T_{{\mathrm{SO}}_{\mathrm{4}}}/K_{{\mathrm{SO}}_{\mathrm{4}}}^{*};Y_{\mathrm{T}}^{\mathrm{F}}=\mathrm{1}/Y_{\mathrm{F}}^{\mathrm{T}}\text{(A6)}& {\displaystyle }{Y}_{\mathrm{F}}^{\mathrm{S}}=\mathrm{1}+T_{{\mathrm{SO}}_{\mathrm{4}}}/K_{{\mathrm{SO}}_{\mathrm{4}}}^{*}+T_{\mathrm{F}}/K_{\mathrm{F}}^{*};Y_{\mathrm{S}}^{\mathrm{F}}=\mathrm{1}/Y_{\mathrm{F}}^{\mathrm{S}}\text{(A7)}& {\displaystyle }{Y}_{\mathrm{S}}^{\mathrm{N}}={\mathit{\gamma }}_{{\mathrm{H}}^{+}};Y_{\mathrm{N}}^{\mathrm{S}}=\mathrm{1}/Y_{\mathrm{S}}^{\mathrm{N}}\end{array}$$

Here, ${\mathit{\gamma }}_{{\mathrm{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:

$$\begin{array}{}\text{(A8)}& {\mathrm{pH}}_B={\mathrm{pH}}_A+\mathrm{p}{Y}_A^B={\mathrm{pH}}_A-{\log }_{\mathrm{10}}\left(Y_A^B\right),\end{array}$$

or alternatively and equivalently

$$\begin{array}{}\text{(A9)}& [{\mathrm{H}}^{+}{]}_B=Y_A^B[{\mathrm{H}}^{+}{]}_A.\end{array}$$

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 (A_T) is calculated as the sum of all its components (Dickson, 1981; Wolf-Gladrow et al., 2007; Sharp and Byrne, 2020).

$$\begin{array}{}\text{(B1)}& \begin{array}{rl}{A}_{\mathrm{T}}& =A_{\mathrm{w}}+A_{\mathrm{C}}+A_{\mathrm{B}}+A_{\mathrm{P}}+A_{\mathrm{Si}}+A_{{\mathrm{NH}}_{\mathrm{3}}}+A_{{\mathrm{H}}_{\mathrm{2}}\mathrm{S}}\\ & +A_{{\mathrm{SO}}_{\mathrm{4}}}+A_{\mathrm{F}}+A_{\mathit{\alpha }}+A_{\mathit{\beta }}\end{array}\end{array}$$

Equations for all the individual alkalinity components (A_C, A_B, etc.) are given in the subsequent sections in terms of pH-independent total substance contents (T_C, T_B, etc.) and [H⁺].

B2 Water

$$\begin{array}{}\text{(BR1)}& {\mathrm{H}}_{\mathrm{2}}\mathrm{O}\rightleftharpoons {\mathrm{OH}}^{-}+{\mathrm{H}}^{+};K_{\mathrm{w}}^{*}=\left[{\mathrm{OH}}^{-}\right]\left[{\mathrm{H}}^{+}\right]\end{array}$$

$$\begin{array}{}\text{(B2)}& A_{\mathrm{w}}=\left[{\mathrm{OH}}^{-}\right]-\left[{\mathrm{H}}^{+}\right]={\displaystyle \frac{K_{\mathrm{w}}^{*}}{\left[{\mathrm{H}}^{+}\right]}}-\left[{\mathrm{H}}^{+}\right]\end{array}$$

B3 Carbonic acid

$$\begin{array}{}\text{(B3)}& T_{\mathrm{C}}=\left[{\mathrm{CO}}_{\mathrm{2}}\left(\mathrm{aq}\right)\right]+\left[{\mathrm{HCO}}_{\mathrm{3}}^{-}\right]+\left[{\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}\right]\end{array}$$

$$\begin{array}{}\text{(BR2)}& {\displaystyle }\begin{array}{rl}& {\mathrm{CO}}_{\mathrm{2}}\left(\mathrm{aq}\right)+{\mathrm{H}}_{\mathrm{2}}\mathrm{O}\rightleftharpoons {\mathrm{HCO}}_{\mathrm{3}}^{-}+{\mathrm{H}}^{+};\\ & K_{\mathrm{1}}^{*}={\displaystyle \frac{\left[{\mathrm{HCO}}_{\mathrm{3}}^{-}\right]\left[{\mathrm{H}}^{+}\right]}{\left[{\mathrm{CO}}_{\mathrm{2}}\left(\mathrm{aq}\right)\right]}}\end{array}\text{(BR3)}& {\displaystyle }{\mathrm{HCO}}_{\mathrm{3}}^{-}\rightleftharpoons {\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}+{\mathrm{H}}^{+};K_{\mathrm{2}}^{*}={\displaystyle \frac{\left[{\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}\right]\left[{\mathrm{H}}^{+}\right]}{\left[{\mathrm{HCO}}_{\mathrm{3}}^{-}\right]}}\end{array}$$

A_C can be expressed in terms of [H⁺] and any of T_C, $f_{{\mathrm{CO}}_{\mathrm{2}}}$, $\left[{\mathrm{HCO}}_{\mathrm{3}}^{-}\right]$, or $\left[{\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}\right]$.

$$\begin{array}{}\text{(B4)}& {\displaystyle }{A}_{\mathrm{C}}=\left[{\mathrm{HCO}}_{\mathrm{3}}^{-}\right]+\mathrm{2}\left[{\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}\right]\text{(B5)}& {\displaystyle }{A}_{\mathrm{C}}\left(\right[{\mathrm{H}}^{+}],T_{\mathrm{C}})={\displaystyle \frac{T_{\mathrm{C}}{K}_{\mathrm{1}}^{*}\left(\right[{\mathrm{H}}^{+}]+\mathrm{2}{K}_{\mathrm{2}}^{*})}{K_{\mathrm{1}}^{*}{K}_{\mathrm{2}}^{*}+K_{\mathrm{1}}^{*}\left[{\mathrm{H}}^{+}\right]+[{\mathrm{H}}^{+}{]}^{\mathrm{2}}}}\text{(B6)}& {\displaystyle }{A}_{\mathrm{C}}\left(\right[{\mathrm{H}}^{+}],f_{{\mathrm{CO}}_{\mathrm{2}}})={\displaystyle \frac{f_{{\mathrm{CO}}_{\mathrm{2}}}{K}_{\mathrm{0}}^{*}{K}_{\mathrm{1}}^{*}\left(\right[{\mathrm{H}}^{+}]+\mathrm{2}{K}_{\mathrm{2}}^{*})}{[{\mathrm{H}}^{+}{]}^{\mathrm{2}}}}\text{(B7)}& {\displaystyle }{A}_{\mathrm{C}}\left(\right[{\mathrm{H}}^{+}],[{\mathrm{HCO}}_{\mathrm{3}}^{-}\left]\right)=\left[{\mathrm{HCO}}_{\mathrm{3}}^{-}\right]+{\displaystyle \frac{\mathrm{2}{K}_{\mathrm{2}}^{*}\left[{\mathrm{HCO}}_{\mathrm{3}}^{-}\right]}{\left[{\mathrm{H}}^{+}\right]}}\text{(B8)}& {\displaystyle }{A}_{\mathrm{C}}\left(\right[{\mathrm{H}}^{+}],[{\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}\left]\right)={\displaystyle \frac{\left[{\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}\right]\left[{\mathrm{H}}^{+}\right]}{K_{\mathrm{2}}^{*}}}+\mathrm{2}\left[{\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}\right]\end{array}$$

Undissociated H₂CO₃ is considered negligible and thus not explicitly modelled, but rather implicitly included as part of the [CO₂(aq)] term (Zeebe and Wolf-Gladrow, 2001).

B4 Boric acid

$$\begin{array}{}\text{(B9)}& T_{\mathrm{B}}=\left[\mathrm{B}(\mathrm{OH}{)}_{\mathrm{3}}\right]+\left[\mathrm{B}(\mathrm{OH}{)}_{\mathrm{4}}^{-}\right]\end{array}$$

$$\begin{array}{}\text{(BR4)}& \begin{array}{rl}& \mathrm{B}(\mathrm{OH}{)}_{\mathrm{3}}+{\mathrm{H}}_{\mathrm{2}}\mathrm{O}\rightleftharpoons \mathrm{B}(\mathrm{OH}{)}_{\mathrm{4}}^{-}+{\mathrm{H}}^{+};\\ & K_{\mathrm{B}}^{*}={\displaystyle \frac{\left[\mathrm{B}(\mathrm{OH}{)}_{\mathrm{4}}^{-}\right]\left[{\mathrm{H}}^{+}\right]}{\left[\mathrm{B}(\mathrm{OH}{)}_{\mathrm{3}}\right]}}\end{array}\end{array}$$

$$\begin{array}{}\text{(B10)}& A_{\mathrm{B}}=\left[\mathrm{B}(\mathrm{OH}{)}_{\mathrm{4}}^{-}\right]={\displaystyle \frac{T_{\mathrm{B}}{K}_{\mathrm{B}}^{*}}{K_{\mathrm{B}}^{*}+\left[{\mathrm{H}}^{+}\right]}}\end{array}$$

B5 Phosphoric acid

$$\begin{array}{}\text{(B11)}& T_{\mathrm{P}}=\left[{\mathrm{H}}_{\mathrm{3}}{\mathrm{PO}}_{\mathrm{4}}\right]+\left[{\mathrm{H}}_{\mathrm{2}}{\mathrm{PO}}_{\mathrm{4}}^{-}\right]+\left[{\mathrm{HPO}}_{\mathrm{4}}^{\mathrm{2}-}\right]+\left[{\mathrm{PO}}_{\mathrm{4}}^{\mathrm{3}-}\right]\end{array}$$

$$\begin{array}{}\text{(BR5)}& {\displaystyle }{\mathrm{H}}_{\mathrm{3}}{\mathrm{PO}}_{\mathrm{4}}\rightleftharpoons {\mathrm{H}}_{\mathrm{2}}{\mathrm{PO}}_{\mathrm{4}}^{-}+{\mathrm{H}}^{+};K_{\mathrm{P}\mathrm{1}}^{*}={\displaystyle \frac{\left[{\mathrm{H}}_{\mathrm{2}}{\mathrm{PO}}_{\mathrm{4}}^{-}\right]\left[{\mathrm{H}}^{+}\right]}{\left[{\mathrm{H}}_{\mathrm{3}}{\mathrm{PO}}_{\mathrm{4}}\right]}}\text{(BR6)}& {\displaystyle }{\mathrm{H}}_{\mathrm{2}}{\mathrm{PO}}_{\mathrm{4}}^{-}\rightleftharpoons {\mathrm{HPO}}_{\mathrm{4}}^{\mathrm{2}-}+{\mathrm{H}}^{+};K_{\mathrm{P}\mathrm{2}}^{*}={\displaystyle \frac{\left[{\mathrm{HPO}}_{\mathrm{4}}^{\mathrm{2}-}\right]\left[{\mathrm{H}}^{+}\right]}{\left[{\mathrm{H}}_{\mathrm{2}}{\mathrm{PO}}_{\mathrm{4}}^{-}\right]}}\text{(BR7)}& {\displaystyle }{\mathrm{HPO}}_{\mathrm{4}}^{\mathrm{2}-}\rightleftharpoons {\mathrm{PO}}_{\mathrm{4}}^{\mathrm{3}-}+{\mathrm{H}}^{+};K_{\mathrm{P}\mathrm{3}}^{*}={\displaystyle \frac{\left[{\mathrm{PO}}_{\mathrm{4}}^{\mathrm{3}-}\right]\left[{\mathrm{H}}^{+}\right]}{\left[{\mathrm{HPO}}_{\mathrm{4}}^{\mathrm{2}-}\right]}}\end{array}$$

$$\begin{array}{}\text{(B12)}& \begin{array}{rl}{A}_{\mathrm{P}}& =\left[{\mathrm{HPO}}_{\mathrm{4}}^{\mathrm{2}-}\right]+\mathrm{2}\left[{\mathrm{PO}}_{\mathrm{4}}^{\mathrm{3}-}\right]-\left[{\mathrm{H}}_{\mathrm{3}}{\mathrm{PO}}_{\mathrm{4}}\right]\\ & ={\displaystyle \frac{T_{\mathrm{P}}\left(K_{\mathrm{P}\mathrm{1}}^{*}{K}_{\mathrm{P}\mathrm{2}}^{*}\right[{\mathrm{H}}^{+}]+\mathrm{2}{K}_{\mathrm{P}\mathrm{1}}^{*}{K}_{\mathrm{P}\mathrm{2}}^{*}{K}_{\mathrm{P}\mathrm{3}}^{*}-[{\mathrm{H}}^{+}{]}^{\mathrm{3}})}{K_{\mathrm{P}\mathrm{1}}^{*}{K}_{\mathrm{P}\mathrm{2}}^{*}{K}_{\mathrm{P}\mathrm{3}}^{*}+K_{\mathrm{P}\mathrm{1}}^{*}{K}_{\mathrm{P}\mathrm{2}}^{*}\left[{\mathrm{H}}^{+}\right]+K_{\mathrm{P}\mathrm{1}}^{*}[{\mathrm{H}}^{+}{]}^{\mathrm{2}}+[{\mathrm{H}}^{+}{]}^{\mathrm{3}}}}\end{array}\end{array}$$

B6 Orthosilicic acid

$$\begin{array}{}\text{(B13)}& T_{\mathrm{Si}}=\left[{\mathrm{H}}_{\mathrm{4}}{\mathrm{SiO}}_{\mathrm{4}}\right]+\left[{\mathrm{H}}_{\mathrm{3}}{\mathrm{SiO}}_{\mathrm{3}}^{-}\right]\end{array}$$

$$\begin{array}{}\text{(BR8)}& \begin{array}{rl}& {\mathrm{H}}_{\mathrm{4}}{\mathrm{SiO}}_{\mathrm{4}}\rightleftharpoons {\mathrm{H}}_{\mathrm{3}}{\mathrm{SiO}}_{\mathrm{4}}^{-}+{\mathrm{H}}^{+};\\ & K_{\mathrm{Si}}^{*}={\displaystyle \frac{\left[{\mathrm{H}}_{\mathrm{3}}{\mathrm{SiO}}_{\mathrm{4}}^{-}\right]\left[{\mathrm{H}}^{+}\right]}{\left[{\mathrm{H}}_{\mathrm{4}}{\mathrm{SiO}}_{\mathrm{4}}\right]}}\end{array}\end{array}$$

$$\begin{array}{}\text{(B14)}& A_{\mathrm{Si}}=\left[{\mathrm{H}}_{\mathrm{3}}{\mathrm{SiO}}_{\mathrm{4}}^{-}\right]={\displaystyle \frac{T_{\mathrm{Si}}{K}_{\mathrm{Si}}^{*}}{K_{\mathrm{Si}}^{*}+\left[{\mathrm{H}}^{+}\right]}}\end{array}$$

Further deprotonation of ${\mathrm{H}}_{\mathrm{3}}{\mathrm{SiO}}_{\mathrm{4}}^{-}$ is considered negligible and thus not modelled.

B7 Ammonium

$$\begin{array}{}\text{(B15)}& T_{{\mathrm{NH}}_{\mathrm{3}}}=\left[{\mathrm{NH}}_{\mathrm{3}}\right]+\left[{\mathrm{NH}}_{\mathrm{4}}^{+}\right]\end{array}$$

$$\begin{array}{}\text{(BR9)}& {\mathrm{NH}}_{\mathrm{4}}^{+}\rightleftharpoons {\mathrm{NH}}_{\mathrm{3}}+{\mathrm{H}}^{+};K_{{\mathrm{NH}}_{\mathrm{3}}}^{*}={\displaystyle \frac{\left[{\mathrm{NH}}_{\mathrm{3}}\right]\left[{\mathrm{H}}^{+}\right]}{\left[{\mathrm{NH}}_{\mathrm{4}}^{+}\right]}}\end{array}$$

$$\begin{array}{}\text{(B16)}& A_{{\mathrm{NH}}_{\mathrm{3}}}=\left[{\mathrm{NH}}_{\mathrm{3}}\right]={\displaystyle \frac{T_{{\mathrm{NH}}_{\mathrm{3}}}{K}_{{\mathrm{NH}}_{\mathrm{3}}}^{*}}{K_{{\mathrm{NH}}_{\mathrm{3}}}^{*}+\left[{\mathrm{H}}^{+}\right]}}\end{array}$$

B8 Sulfide

$$\begin{array}{}\text{(B17)}& T_{{\mathrm{H}}_{\mathrm{2}}\mathrm{S}}=\left[{\mathrm{H}}_{\mathrm{2}}\mathrm{S}\right]+\left[{\mathrm{HS}}^{-}\right]\end{array}$$

$$\begin{array}{}\text{(BR10)}& {\mathrm{H}}_{\mathrm{2}}\mathrm{S}\rightleftharpoons {\mathrm{HS}}^{-}+{\mathrm{H}}^{+};K_{{\mathrm{H}}_{\mathrm{2}}\mathrm{S}}^{*}={\displaystyle \frac{\left[{\mathrm{HS}}^{-}\right]\left[{\mathrm{H}}^{+}\right]}{\left[{\mathrm{H}}_{\mathrm{2}}\mathrm{S}\right]}}\end{array}$$

$$\begin{array}{}\text{(B18)}& A_{{\mathrm{H}}_{\mathrm{2}}\mathrm{S}}=\left[{\mathrm{HS}}^{-}\right]={\displaystyle \frac{T_{{\mathrm{H}}_{\mathrm{2}}\mathrm{S}}{K}_{{\mathrm{H}}_{\mathrm{2}}\mathrm{S}}^{*}}{K_{{\mathrm{H}}_{\mathrm{2}}\mathrm{S}}^{*}+\left[{\mathrm{H}}^{+}\right]}}\end{array}$$

Further deprotonation of HS⁻ is considered negligible and thus not modelled (Schoonen and Barnes, 1988).

B9 Sulfate

$$\begin{array}{}\text{(B19)}& T_{{\mathrm{SO}}_{\mathrm{4}}}=\left[{\mathrm{HSO}}_{\mathrm{4}}^{-}\right]+\left[{\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}\right]\end{array}$$

$$\begin{array}{}\text{(BR11)}& {\mathrm{HSO}}_{\mathrm{4}}^{-}\rightleftharpoons {\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}+{\mathrm{H}}^{+};K_{{\mathrm{SO}}_{\mathrm{4}}}^{*}={\displaystyle \frac{\left[{\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}\right]\left[{\mathrm{H}}^{+}\right]}{\left[{\mathrm{HSO}}_{\mathrm{4}}^{-}\right]}}\end{array}$$

$$\begin{array}{}\text{(B20)}& A_{{\mathrm{SO}}_{\mathrm{4}}}=-\left[{\mathrm{HSO}}_{\mathrm{4}}^{-}\right]={\displaystyle \frac{-T_{{\mathrm{SO}}_{\mathrm{4}}}}{\mathrm{1}+K_{{\mathrm{SO}}_{\mathrm{4}}}^{*}/\left[{\mathrm{H}}^{+}\right]}}\end{array}$$

Undissociated H₂SO₄ is considered negligible and thus not modelled.

B10 Fluoride

$$\begin{array}{}\text{(B21)}& T_{\mathrm{F}}=\left[\mathrm{HF}\right]+\left[{\mathrm{F}}^{-}\right]\end{array}$$

$$\begin{array}{}\text{(BR12)}& \mathrm{HF}\rightleftharpoons {\mathrm{F}}^{-}+{\mathrm{H}}^{+};K_{\mathrm{F}}^{*}={\displaystyle \frac{\left[{\mathrm{F}}^{-}\right]\left[{\mathrm{H}}^{+}\right]}{\left[\mathrm{HF}\right]}}\end{array}$$

$$\begin{array}{}\text{(B22)}& A_{\mathrm{F}}=-\left[\mathrm{HF}\right]={\displaystyle \frac{-T_{\mathrm{F}}}{\mathrm{1}+K_{\mathrm{F}}^{*}/\left[{\mathrm{H}}^{+}\right]}}\end{array}$$

B11 Arbitrary additional components

$$\begin{array}{}\text{(B23)}& T_{\mathit{\alpha }}=\left[\mathrm{H}\mathit{\alpha }\right]+\left[{\mathit{\alpha }}^{-}\right]\end{array}$$

$$\begin{array}{}\text{(BR13)}& \mathrm{H}\mathit{\alpha }\rightleftharpoons {\mathit{\alpha }}^{-}+{\mathrm{H}}^{+};K_{\mathit{\alpha }}^{*}={\displaystyle \frac{\left[{\mathit{\alpha }}^{-}\right]\left[{\mathrm{H}}^{+}\right]}{\left[\mathrm{H}\mathit{\alpha }\right]}}\end{array}$$

$$\begin{array}{}\text{(B24)}& A_{\mathit{\alpha }}=\left\{\begin{array}{ll}-\left[\mathrm{H}\mathit{\alpha }\right]& \text{for }-{\log }_{\mathrm{10}}\left(K_{\mathit{\alpha }}^{*}\right)\le \mathrm{4.5}\\ +\left[{\mathit{\alpha }}^{-}\right]& \text{for }-{\log }_{\mathrm{10}}\left(K_{\mathit{\alpha }}^{*}\right)>\mathrm{4.5}\end{array}\right.\end{array}$$

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_{\mathit{\alpha }}^{*}$ and $K_{\mathit{\beta }}^{*}$ values, with a zero level of protons corresponding to a pK^* of 4.5 (Wolf-Gladrow et al., 2007).

Though the definition of alkalinity (Dickson, 1981) 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-Gladrow, 2001; 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

$$\begin{array}{}\text{(C1)}& \mathrm{pH}=-{\log }_{\mathrm{10}}\left[{\mathrm{H}}^{+}\right]\end{array}$$

regardless of which pH scale is being used.

C1.2 Known $p_{{\mathrm{CO}}_{\mathrm{2}}}$, $x_{{\mathrm{CO}}_{\mathrm{2}}}$ or [CO₂(aq)]

If one of $p_{{\mathrm{CO}}_{\mathrm{2}}}$, $x_{{\mathrm{CO}}_{\mathrm{2}}}$, or [CO₂(aq)] is in the known parameter pair, then its values are first converted to $f_{{\mathrm{CO}}_{\mathrm{2}}}$ as follows.

For known $p_{{\mathrm{CO}}_{\mathrm{2}}}$,

$$\begin{array}{}\text{(C2)}& f_{{\mathrm{CO}}_{\mathrm{2}}}=G{p}_{{\mathrm{CO}}_{\mathrm{2}}},\end{array}$$

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

For known $x_{{\mathrm{CO}}_{\mathrm{2}}}$,

$$\begin{array}{}\text{(C3)}& f_{{\mathrm{CO}}_{\mathrm{2}}}=G{P}_{\mathrm{v}}{x}_{{\mathrm{CO}}_{\mathrm{2}}},\end{array}$$

where P_v is the humidity correction (Table 3):

$$\begin{array}{}\text{(C4)}& P_{\mathrm{v}}=P_{\mathrm{a}}-p_{\mathrm{w}},\end{array}$$

in which P_a is total atmospheric pressure (assumed to be 1 atm unless a different value is provided by the user) and p_w is the water vapour pressure (Weiss and Price, 1980).

For known [CO₂(aq)],

$$\begin{array}{}\text{(C5)}& f_{{\mathrm{CO}}_{\mathrm{2}}}={\displaystyle \frac{\left[{\mathrm{CO}}_{\mathrm{2}}\left(\mathrm{aq}\right)\right]}{K_{\mathrm{0}}^{*}}},\end{array}$$

where $K_{\mathrm{0}}^{*}$ is the solubility factor for CO₂ (Table 3).

The calculation steps given below for $f_{{\mathrm{CO}}_{\mathrm{2}}}$ are then followed to solve the core marine carbonate system. Afterwards, $p_{{\mathrm{CO}}_{\mathrm{2}}}$, $x_{{\mathrm{CO}}_{\mathrm{2}}}$, and [CO₂(aq)] are calculated where they were not in the original known parameter pair: $p_{{\mathrm{CO}}_{\mathrm{2}}}$ and $x_{{\mathrm{CO}}_{\mathrm{2}}}$ are calculated using Eqs. (C2) and (C3), while [CO₂(aq)] is calculated by difference using the definition of T_C in Eq. (B3).

C2 Solving routines

C2.1 From A_T and T_C

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 A_T(pH_n,v) term in Eq. (2) is calculated from Eq. (B1) for A_T substituting in Eq. (B5) for the A_C term. Equation (2) is the automatically differentiated with respect to pH to obtain the $\mathrm{\Delta }{A}_{\mathrm{T}}^{\prime }$ term in Eq. (1).

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

$$\begin{array}{}\text{(C6)}& {\displaystyle }{f}_{{\mathrm{CO}}_{\mathrm{2}}}={\displaystyle \frac{T_{\mathrm{C}}[{\mathrm{H}}^{+}{]}^{\mathrm{2}}}{K_{\mathrm{0}}^{*}\left(\right[{\mathrm{H}}^{+}{]}^{\mathrm{2}}+K_{\mathrm{1}}^{*}\left[{\mathrm{H}}^{+}\right]+K_{\mathrm{1}}^{*}{K}_{\mathrm{2}}^{*})}}\text{(C7)}& {\displaystyle }\left[{\mathrm{HCO}}_{\mathrm{3}}^{-}\right]={\displaystyle \frac{T_{\mathrm{C}}{K}_{\mathrm{1}}^{*}\left[{\mathrm{H}}^{+}\right]}{[{\mathrm{H}}^{+}{]}^{\mathrm{2}}+K_{\mathrm{1}}^{*}[{\mathrm{H}}^{+}]+K_{\mathrm{1}}^{*}{K}_{\mathrm{2}}^{*}}}\text{(C8)}& {\displaystyle }\left[{\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}\right]={\displaystyle \frac{T_{\mathrm{C}}{K}_{\mathrm{1}}^{*}{K}_{\mathrm{2}}^{*}}{[{\mathrm{H}}^{+}{]}^{\mathrm{2}}+K_{\mathrm{1}}^{*}[{\mathrm{H}}^{+}]+K_{\mathrm{1}}^{*}{K}_{\mathrm{2}}^{*}}}\end{array}$$

C2.2 From A_T and pH

First, we determine A_C from known A_T and pH by using Eq. (B1). T_C is then calculated from A_C.

$$\begin{array}{}\text{(C9)}& T_{\mathrm{C}}={\displaystyle \frac{A_{\mathrm{C}}\left(\right[{\mathrm{H}}^{+}{]}^{\mathrm{2}}+K_{\mathrm{1}}^{*}\left[{\mathrm{H}}^{+}\right]+K_{\mathrm{1}}^{*}{K}_{\mathrm{2}}^{*})}{K_{\mathrm{1}}^{*}\left(\right[{\mathrm{H}}^{+}]+\mathrm{2}{K}_{\mathrm{2}}^{*})}}\end{array}$$

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

There is an upper limit on pH for each given A_T value, above which negative A_C 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 T_C (and all other results calculated from it) instead of a negative value.

C2.3 From A_T and $f_{{\mathrm{CO}}_{\mathrm{2}}}$

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 A_T(pH_n,v) term in Eq. (2) is calculated from Eq. (B1) for A_T substituting in Eq. (B6) for the A_C term. Equation (2) is the automatically differentiated with respect to pH to obtain the $\mathrm{\Delta }{A}_{\mathrm{T}}^{\prime }$ term in Eq. (1).

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

C2.4 From A_T and $\left[{\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}\right]$

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 A_T(pH_n,v) term in Eq. (2) is calculated from Eq. (B1) for A_T substituting in Eq. (B8) for the A_C term. Equation (2) is the automatically differentiated with respect to pH to obtain the $\mathrm{\Delta }{A}_{\mathrm{T}}^{\prime }$ term in Eq. (1). The lower of the two pH roots is returned by default, as discussed in Sect. 5.2.2.

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

C2.5 From A_T and $\left[{\mathrm{HCO}}_{\mathrm{3}}^{-}\right]$

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 A_T(pH_n,v) term in Eq. (2) is calculated from Eq. (B1) for A_T substituting in Eq. (B7) for the A_C term. Equation (2) is the automatically differentiated with respect to pH to obtain the $\mathrm{\Delta }{A}_{\mathrm{T}}^{\prime }$ term in Eq. (1).

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

C2.6 From T_C and pH

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

C2.7 From T_C and $f_{{\mathrm{CO}}_{\mathrm{2}}}$

First, pH is calculated from T_C and $f_{{\mathrm{CO}}_{\mathrm{2}}}$ using

$$\begin{array}{}\text{(C10)}& \left[{\mathrm{H}}^{+}\right]={\displaystyle \frac{K_{\mathrm{1}}^{*}r+\sqrt{(K_{\mathrm{1}}^{*}r{)}^{\mathrm{2}}+\mathrm{4}(\mathrm{1}-r)K_{\mathrm{1}}^{*}{K}_{\mathrm{2}}^{*}r}}{\mathrm{2}(\mathrm{1}-r)}},\end{array}$$

where

$$\begin{array}{}\text{(C11)}& r=K_{\mathrm{0}}^{*}\cdot f_{{\mathrm{CO}}_{\mathrm{2}}}/T_{\mathrm{C}}.\end{array}$$

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

C2.8 From T_C and $\left[{\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}\right]$

First, pH is calculated from T_C and $\left[{\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}\right]$ using

$$\begin{array}{}\text{(C12)}& \left[{\mathrm{H}}^{+}\right]={\displaystyle \frac{-K_{\mathrm{1}}^{*}+\sqrt{K_{\mathrm{1}}^{*\mathrm{2}}-\mathrm{4}{K}_{\mathrm{1}}^{*}{K}_{\mathrm{2}}^{*}(\mathrm{1}-T_{\mathrm{C}}/[{\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}\left]\right)}}{\mathrm{2}}}.\end{array}$$

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

C2.9 From T_C and $\left[{\mathrm{HCO}}_{\mathrm{3}}^{-}\right]$

First, pH is calculated from T_C and $\left[{\mathrm{HCO}}_{\mathrm{3}}^{-}\right]$ using

$$\begin{array}{}\text{(C13)}& \left[{\mathrm{H}}^{+}\right]={\displaystyle \frac{\begin{array}{c}{T}_{\mathrm{C}}-\left[{\mathrm{HCO}}_{\mathrm{3}}^{-}\right]-\\ \sqrt{\left(\right[{\mathrm{HCO}}_{\mathrm{3}}^{-}]-T_{\mathrm{C}}{)}^{\mathrm{2}}-\mathrm{4}[{\mathrm{HCO}}_{\mathrm{3}}^{-}{]}^{\mathrm{2}}{K}_{\mathrm{2}}^{*}/K_{\mathrm{1}}^{*}}\end{array}}{\mathrm{2}\left[{\mathrm{HCO}}_{\mathrm{3}}^{-}\right]/K_{\mathrm{1}}^{*}}}.\end{array}$$

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

C2.10 From pH and $f_{{\mathrm{CO}}_{\mathrm{2}}}$

First, T_C is calculated from pH and $f_{{\mathrm{CO}}_{\mathrm{2}}}$ using

$$\begin{array}{}\text{(C14)}& T_{\mathrm{C}}={\displaystyle \frac{K_{\mathrm{0}}^{*}\cdot f_{{\mathrm{CO}}_{\mathrm{2}}}\left([{\mathrm{H}}^{+}{]}^{\mathrm{2}}+K_{\mathrm{1}}^{*}[{\mathrm{H}}^{+}]+K_{\mathrm{1}}^{*}{K}_{\mathrm{2}}^{*}\right)}{[{\mathrm{H}}^{+}{]}^{\mathrm{2}}}}.\end{array}$$

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

C2.11 From pH and $\left[{\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}\right]$

First, $f_{{\mathrm{CO}}_{\mathrm{2}}}$ is calculated from pH and $\left[{\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}\right]$ using

$$\begin{array}{}\text{(C15)}& f_{{\mathrm{CO}}_{\mathrm{2}}}={\displaystyle \frac{\left[{\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}\right][{\mathrm{H}}^{+}{]}^{\mathrm{2}}}{K_{\mathrm{0}}^{*}{K}_{\mathrm{1}}^{*}{K}_{\mathrm{2}}^{*}}}.\end{array}$$

T_C is then calculated from pH and $f_{{\mathrm{CO}}_{\mathrm{2}}}$ using Eq. (C14). Finally, A_T and $\left[{\mathrm{HCO}}_{\mathrm{3}}^{-}\right]$ are calculated from T_C and pH using Eqs. (B1) and (C7), respectively.

C2.12 From pH and $\left[{\mathrm{HCO}}_{\mathrm{3}}^{-}\right]$

First, T_C is calculated from pH and $\left[{\mathrm{HCO}}_{\mathrm{3}}^{-}\right]$ using

$$\begin{array}{}\text{(C16)}& T_{\mathrm{C}}=\left[{\mathrm{HCO}}_{\mathrm{3}}^{-}\right]\left(\mathrm{1}+{\displaystyle \frac{\left[{\mathrm{H}}^{+}\right]}{K_{\mathrm{1}}^{*}}}+{\displaystyle \frac{K_{\mathrm{2}}^{*}}{\left[{\mathrm{H}}^{+}\right]}}\right).\end{array}$$

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

C2.13 From $f_{{\mathrm{CO}}_{\mathrm{2}}}$ and $\left[{\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}\right]$

First, pH is calculated from $f_{{\mathrm{CO}}_{\mathrm{2}}}$ and $\left[{\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}\right]$ using

$$\begin{array}{}\text{(C17)}& \left[{\mathrm{H}}^{+}\right]=\sqrt{{\displaystyle \frac{K_{\mathrm{0}}^{*}{K}_{\mathrm{1}}^{*}{K}_{\mathrm{2}}^{*}\cdot f_{{\mathrm{CO}}_{\mathrm{2}}}}{\left[{\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}\right]}}}.\end{array}$$

T_C is then calculated from pH and $f_{{\mathrm{CO}}_{\mathrm{2}}}$ using Eq. (C14). Finally, A_T and $\left[{\mathrm{HCO}}_{\mathrm{3}}^{-}\right]$ are calculated from T_C and pH using Eqs. (B1) and (C7), respectively.

C2.14 From $f_{{\mathrm{CO}}_{\mathrm{2}}}$ and $\left[{\mathrm{HCO}}_{\mathrm{3}}^{-}\right]$

First, $\left[{\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}\right]$ is calculated from $f_{{\mathrm{CO}}_{\mathrm{2}}}$ and $\left[{\mathrm{HCO}}_{\mathrm{3}}^{-}\right]$ using

$$\begin{array}{}\text{(C18)}& \left[{\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}\right]={\displaystyle \frac{[{\mathrm{HCO}}_{\mathrm{3}}^{-}{]}^{\mathrm{2}}{K}_{\mathrm{2}}^{*}}{K_{\mathrm{0}}^{*}{K}_{\mathrm{1}}^{*}\cdot f_{{\mathrm{CO}}_{\mathrm{2}}}}}.\end{array}$$

The pH is then calculated from $f_{{\mathrm{CO}}_{\mathrm{2}}}$ and $\left[{\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}\right]$ using Eq. (C17). Next, T_C is calculated from pH and $f_{{\mathrm{CO}}_{\mathrm{2}}}$ using Eq. (C14). Finally, A_T is calculated from T_C and pH using Eq. (B1).

C2.15 From $\left[{\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}\right]$ and $\left[{\mathrm{HCO}}_{\mathrm{3}}^{-}\right]$

First, $f_{{\mathrm{CO}}_{\mathrm{2}}}$ is calculated from $\left[{\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}\right]$ and $\left[{\mathrm{HCO}}_{\mathrm{3}}^{-}\right]$ using

$$\begin{array}{}\text{(C19)}& f_{{\mathrm{CO}}_{\mathrm{2}}}={\displaystyle \frac{[{\mathrm{HCO}}_{\mathrm{3}}^{-}{]}^{\mathrm{2}}{K}_{\mathrm{2}}^{*}}{K_{\mathrm{0}}^{*}{K}_{\mathrm{1}}^{*}\left[{\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}\right]}}.\end{array}$$

The pH is then calculated from $f_{{\mathrm{CO}}_{\mathrm{2}}}$ and $\left[{\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}\right]$ using Eq. (C17). Next, T_C is calculated from pH and $f_{{\mathrm{CO}}_{\mathrm{2}}}$ using Eq. (C14). Finally, A_T is calculated from T_C and pH using Eq. (B1).

Appendix D: Other marine carbonate system variables

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

$$\begin{array}{}\text{(D1)}& \mathrm{\Omega }={\displaystyle \frac{\left[{\mathrm{Ca}}^{\mathrm{2}}+\right]\left[{\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}\right]}{K_{\mathrm{sp}}^{*}}},\end{array}$$

where $K_{\mathrm{sp}}^{*}$ 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.

$$\begin{array}{}\text{(D2)}& \mathrm{SIR}={\displaystyle \frac{\left[{\mathrm{HCO}}_{\mathrm{3}}^{-}\right]}{\left[{\mathrm{H}}^{+}\right]}}\end{array}$$

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):

• $(\partial T_{\mathrm{C}}/\partial \mathrm{pH}{)}_{A_{\mathrm{T}}}$ by AD of Eq. (C9) with respect to pH;

• $(\partial A_{\mathrm{T}}/\partial \mathrm{pH}{)}_{T_{\mathrm{C}}}$ by AD of Eq. (B1), substituting A_C by Eq. (B5), with respect to pH;

• $(\partial \ln [{\mathrm{CO}}_{\mathrm{2}}\left(\mathrm{aq}\right)]/\partial \mathrm{pH}{)}_{T_{\mathrm{C}}}$ by taking the natural log of the product of $K_{\mathrm{0}}^{*}$ and Eq. (C6), then AD with respect to pH;

• $(\partial \ln [{\mathrm{CO}}_{\mathrm{2}}\left(\mathrm{aq}\right)]/\partial \mathrm{pH}{)}_{A_{\mathrm{T}}}$ by taking the natural log of the product of $K_{\mathrm{0}}^{*}$ and Eq. (C6), substituting T_C by Eq. (9), then AD with respect to pH.

The buffer factors ${\mathit{\gamma }}_{T_{\mathrm{C}}}$, ${\mathit{\gamma }}_{A_{\mathrm{T}}}$, ${\mathit{\beta }}_{T_{\mathrm{C}}}$, and ${\mathit{\beta }}_{A_{\mathrm{T}}}$ are thus defined (Egleston et al., 2010) and calculated in PyCO2SYS.

$$\begin{array}{}\text{(E1)}& {\displaystyle }\begin{array}{rl}{\mathit{\gamma }}_{T_{\mathrm{C}}}& ={\left({\displaystyle \frac{\partial \ln \left[{\mathrm{CO}}_{\mathrm{2}}\left(\mathrm{aq}\right)\right]}{\partial T_{\mathrm{C}}}}\right)}_{A_{\mathrm{T}}}^{-\mathrm{1}}\\ & ={\left({\displaystyle \frac{\partial T_{\mathrm{C}}}{\partial \mathrm{pH}}}\right)}_{A_{\mathrm{T}}}{\left({\displaystyle \frac{\partial \ln \left[{\mathrm{CO}}_{\mathrm{2}}\left(\mathrm{aq}\right)\right]}{\partial \mathrm{pH}}}\right)}_{A_{\mathrm{T}}}^{-\mathrm{1}}\end{array}\text{(E2)}& {\displaystyle }\begin{array}{rl}{\mathit{\gamma }}_{A_{\mathrm{T}}}& ={\left({\displaystyle \frac{\partial \ln \left[{\mathrm{CO}}_{\mathrm{2}}\left(\mathrm{aq}\right)\right]}{\partial A_{\mathrm{T}}}}\right)}_{T_{\mathrm{C}}}^{-\mathrm{1}}\\ & ={\left({\displaystyle \frac{\partial A_{\mathrm{T}}}{\partial \mathrm{pH}}}\right)}_{T_{\mathrm{C}}}{\left({\displaystyle \frac{\partial \ln \left[{\mathrm{CO}}_{\mathrm{2}}\left(\mathrm{aq}\right)\right]}{\partial \mathrm{pH}}}\right)}_{T_{\mathrm{C}}}^{-\mathrm{1}}\end{array}\text{(E3)}& {\displaystyle }{\mathit{\beta }}_{T_{\mathrm{C}}}={\left({\displaystyle \frac{\partial \ln \left[{\mathrm{H}}^{+}\right]}{\partial T_{\mathrm{C}}}}\right)}_{A_{\mathrm{T}}}^{-\mathrm{1}}=-{\log }_{\mathrm{10}}\left(e\right){\left({\displaystyle \frac{\partial T_{\mathrm{C}}}{\partial \mathrm{pH}}}\right)}_{A_{\mathrm{T}}}\text{(E4)}& {\displaystyle }{\mathit{\beta }}_{A_{\mathrm{T}}}={\left({\displaystyle \frac{\partial \ln \left[{\mathrm{H}}^{+}\right]}{\partial A_{\mathrm{T}}}}\right)}_{T_{\mathrm{C}}}^{-\mathrm{1}}=-{\log }_{\mathrm{10}}\left(e\right){\left({\displaystyle \frac{\partial A_{\mathrm{T}}}{\partial \mathrm{pH}}}\right)}_{T_{\mathrm{C}}}\end{array}$$

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

For the saturation-state buffers ${\mathit{\omega }}_{T_{\mathrm{C}}}$ and ${\mathit{\omega }}_{A_{\mathrm{T}}}$ we also evaluate

• $(\partial \ln \mathrm{\Omega }/\partial [{\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}\left]\right)$ by AD of the natural log of Ω(aragonite), calculated with Eq. (D1), with respect to $\left[{\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}\right]$ (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);

• $(\partial [{\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}]/\partial \mathrm{pH}{)}_{T_{\mathrm{C}}}$ by AD of Eq. (C8) with respect to pH;

• $(\partial [{\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}]/\partial \mathrm{pH}{)}_{A_{\mathrm{T}}}$ by AD of Eq. (C8), substituting T_C by Eq. (C9), with respect to pH.

The buffer factors are then given by the following.

$$\begin{array}{}\text{(E5)}& {\displaystyle }\begin{array}{rl}{\mathit{\omega }}_{T_{\mathrm{C}}}& ={\left({\displaystyle \frac{\partial \ln \mathrm{\Omega }}{\partial T_{\mathrm{C}}}}\right)}_{A_{\mathrm{T}}}^{-\mathrm{1}}\\ & ={\left({\displaystyle \frac{\partial T_{\mathrm{C}}}{\partial \mathrm{pH}}}\right)}_{A_{\mathrm{T}}}{\left({\displaystyle \frac{\partial \ln \mathrm{\Omega }}{\partial \left[{\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}\right]}}\right)}_{A_{\mathrm{T}}}^{-\mathrm{1}}{\left({\displaystyle \frac{\partial \left[{\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}\right]}{\partial \mathrm{pH}}}\right)}_{A_{\mathrm{T}}}^{-\mathrm{1}}\end{array}\text{(E6)}& {\displaystyle }\begin{array}{rl}{\mathit{\omega }}_{A_{\mathrm{T}}}& ={\left({\displaystyle \frac{\partial \ln \mathrm{\Omega }}{\partial A_{\mathrm{T}}}}\right)}_{T_{\mathrm{C}}}^{-\mathrm{1}}\\ & ={\left({\displaystyle \frac{\partial A_{\mathrm{T}}}{\partial \mathrm{pH}}}\right)}_{T_{\mathrm{C}}}{\left({\displaystyle \frac{\partial \ln \mathrm{\Omega }}{\partial \left[{\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}\right]}}\right)}_{T_{\mathrm{C}}}^{-\mathrm{1}}{\left({\displaystyle \frac{\partial \left[{\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}\right]}{\partial \mathrm{pH}}}\right)}_{T_{\mathrm{C}}}^{-\mathrm{1}}\end{array}\end{array}$$

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 (R_F; Broecker et al., 1979) is computed from T_C and ${\mathit{\gamma }}_{T_{\mathrm{C}}}$, with the latter evaluated as described in Sect. E1, following Egleston et al. (2010).

$$\begin{array}{}\text{(E7)}& R_{\mathrm{F}}=\left({\displaystyle \frac{\partial f_{{\mathrm{CO}}_{\mathrm{2}}}}{\partial T_{\mathrm{C}}}}\right)\left({\displaystyle \frac{T_{\mathrm{C}}}{f_{{\mathrm{CO}}_{\mathrm{2}}}}}\right)={\displaystyle \frac{T_{\mathrm{C}}}{{\mathit{\gamma }}_{T_{\mathrm{C}}}}}\end{array}$$

E3 Isocapnic quotient and ψ

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

• $(\partial T_{\mathrm{C}}/\partial \mathrm{pH}{)}_{f_{{\mathrm{CO}}_{\mathrm{2}}}}$ by AD of Eq. (C14) with respect to pH;

• $(\partial A_{\mathrm{T}}/\partial \mathrm{pH}{)}_{f_{{\mathrm{CO}}_{\mathrm{2}}}}$ by AD of Eq. (B1) using Eq. (B6) for the A_C term with respect to pH.

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

$$\begin{array}{}\text{(E8)}& Q={\left({\displaystyle \frac{\partial A_{\mathrm{T}}}{\partial T_{\mathrm{C}}}}\right)}_{f_{{\mathrm{CO}}_{\mathrm{2}}}}={\left({\displaystyle \frac{\partial A_{\mathrm{T}}}{\partial \mathrm{pH}}}\right)}_{f_{{\mathrm{CO}}_{\mathrm{2}}}}{\left({\displaystyle \frac{\partial T_{\mathrm{C}}}{\partial \mathrm{pH}}}\right)}_{f_{{\mathrm{CO}}_{\mathrm{2}}}}^{-\mathrm{1}}\end{array}$$

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

$$\begin{array}{}\text{(E9)}& \mathit{\psi }={\displaystyle \frac{\mathrm{2}}{Q}}-\mathrm{1}\end{array}$$

Appendix F: Initial pH estimate when solving from A_T and T_C

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

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

$$\begin{array}{}\text{(F1)}& A_{\mathrm{CB}}(h,T_{\mathrm{C}})={\displaystyle \frac{T_{\mathrm{C}}{K}_{\mathrm{1}}^{*}(h+\mathrm{2}{K}_{\mathrm{2}}^{*})}{h^{\mathrm{2}}+K_{\mathrm{1}}^{*}h+K_{\mathrm{1}}^{*}{K}_{\mathrm{2}}^{*}}}+{\displaystyle \frac{T_{\mathrm{B}}{K}_{\mathrm{B}}^{*}}{h+K_{\mathrm{B}}^{*}}}.\end{array}$$

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

$$\begin{array}{}\text{(F2)}& P_{T_{\mathrm{C}}}\left(h\right)=h^{\mathrm{3}}+h^{\mathrm{2}}{g}_{\mathrm{2}}\left(T_{\mathrm{C}}\right)+h{g}_{\mathrm{1}}\left(T_{\mathrm{C}}\right)+g_{\mathrm{0}}\left(T_{\mathrm{C}}\right)=\mathrm{0},\end{array}$$

with the following.

$$\begin{array}{}\text{(F3)}& {\displaystyle }{g}_{\mathrm{2}}\left(T_{\mathrm{C}}\right)=K_{\mathrm{B}}^{*}\left(\mathrm{1}-{\displaystyle \frac{T_{\mathrm{B}}}{A_{\mathrm{CB}}}}\right)-K_{\mathrm{1}}^{*}\left(\mathrm{1}-{\displaystyle \frac{T_{\mathrm{C}}}{A_{\mathrm{CB}}}}\right)\text{(F4)}& {\displaystyle }{g}_{\mathrm{1}}\left(T_{\mathrm{C}}\right)=K_{\mathrm{1}}^{*}\left[K_{\mathrm{B}}^{*}\left(\mathrm{1}-{\displaystyle \frac{T_{\mathrm{B}}+T_{\mathrm{C}}}{A_{\mathrm{CB}}}}\right)+K_{\mathrm{2}}^{*}\left(\mathrm{1}-{\displaystyle \frac{\mathrm{2}{T}_{\mathrm{C}}}{A_{\mathrm{CB}}}}\right)\right]\text{(F5)}& {\displaystyle }{g}_{\mathrm{0}}\left(T_{\mathrm{C}}\right)=K_{\mathrm{1}}^{*}{K}_{\mathrm{2}}^{*}{K}_{\mathrm{B}}^{*}\left(\mathrm{1}-{\displaystyle \frac{\mathrm{2}{T}_{\mathrm{C}}+T_{\mathrm{B}}}{A_{\mathrm{CB}}}}\right)\end{array}$$

The initial h value is determined by

$$\begin{array}{}\text{(F6)}& h_{\mathrm{0}}\left(T_{\mathrm{C}}\right)=\left\{\begin{array}{ll}{\mathrm{10}}^{-\mathrm{3}}& \text{for }{A}_{\mathrm{T}}\le \mathrm{0}\\ h_{\min }+\sqrt{-\frac{P_{T_{\mathrm{C}}}\left(h_{\min }\right)}{\sqrt{g_{\mathrm{2}}^{\mathrm{2}}-\mathrm{3}{g}_{\mathrm{1}}}}}& \text{for }{A}_{\mathrm{T}}>\mathrm{0}\\ {\mathrm{10}}^{-\mathrm{10}}& \text{for }{A}_{\mathrm{T}}\ge \mathrm{2}{T}_{\mathrm{C}}+T_{\mathrm{B}},\end{array}\right.\end{array}$$

where h_min is defined in Eq. (10). Negative A_CB is impossible because its equation contains only positive terms, so the equations above cannot be applied if A_T is indeed negative. The default h₀ of 10⁻³ mol kg⁻¹, 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 A_CB is 2T_C+T_B, where T_C is entirely ${\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}$ and T_B is entirely $\mathrm{B}(\mathrm{OH}{)}_{\mathrm{4}}^{-}$. Where A_T is actually higher than this limit of this simplified expression, we expect a high pH (given the dominance of ${\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}$ within T_C) and therefore use an initial estimate pH of 10. Otherwise, h_min in Eq. (F6) is found using Eq. (10). A default h₀ of 10⁻⁷ mol kg⁻¹ (pH 7) is used when $g_{\mathrm{2}}^{\mathrm{2}}-\mathrm{3}{g}_{\mathrm{1}}\le \mathrm{0}$ in Eq. (F6) (Munhoven, 2013).

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 T_C value is used in the final evaluation. Rather than using the “central” T_C value, the change in $p_{{\mathrm{CO}}_{\mathrm{2}}}$ is divided by the adjusted (T_C−ΔT_C).

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 ΔT_C used in the original finite-difference implementation (i.e. 1 µmol kg⁻¹), 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 H1 Fixed Δa values for uncertainty analysis.

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.