the Creative Commons Attribution 4.0 License.
the Creative Commons Attribution 4.0 License.
PyCO2SYS v1.8: marine carbonate system calculations in Python
Matthew P. Humphreys
Ernie R. Lewis
Jonathan D. Sharp
Denis Pierrot
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_{2}), and seawater pH is decreasing as a consequence. Both the exchange of CO_{2} 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_{2} fugacity (${f}_{{\mathrm{CO}}_{\mathrm{2}}}$; or its partial pressure, ${p}_{{\mathrm{CO}}_{\mathrm{2}}}$, or its dryair 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).
The ocean absorbs about a quarter of the anthropogenic carbon dioxide (CO_{2}) currently emitted each year (Friedlingstein et al., 2020). This absorption is a doubleedged sword. Removing CO_{2} from the atmosphere reduces the impact of these emissions on Earth's climate. However, CO_{2} 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_{2} 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_{2}, 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 dryair mole fraction of CO_{2} 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 WolfGladrow, 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, opensource, 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 MSDOS (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 ascloseaspossible translation of CO2SYSMATLAB 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 wellestablished in the research field, we provide a relatively brief summary of the components of PyCO2SYS that are identical to CO2SYSMATLAB in Sect. 2, before describing the areas where PyCO2SYS differs in more detail in Sect. 3. Equations that were inherited from CO2SYSMATLAB 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 CO2SYSMATLAB, before concluding with our perspectives on the outlook for PyCO2SYS and other related software.
The components of PyCO2SYS that have been inherited directly from CO2SYSMATLAB v2.0.5 (Orr et al., 2018), with only the minimal changes needed to translate to Python plus aesthetic code restructuring, are described in this section.
2.1 Units and pH scales
The abundances of all solutes and total alkalinity provided as arguments to PyCO2SYS or returned from it as results are in units of micromoles per kilogram (µmol kg^{−1}), with kilograms (kg) being of the total solution. This means that they are neither concentrations nor molarity values, which are both per unit volume rather than mass, nor are they molality values, which are per kilogram of H_{2}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^{−1}) (Appendix A; Zeebe and WolfGladrow, 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)^{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 zerosalinity freshwater. ^{e} Including the corrections of Waters et al. (2014).
^{a} Depending on user input. ^{b} In GEOSECSTakahashi, phosphate is not included in the definition of total alkalinity; in GEOSECSPeng, 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).

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

GEOSECS cases: GEOSECSTakahashi and GEOSECSPeng. GEOSECSPeng 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.

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.
^{a} As selected by the user. ^{b} Including the corrections of Waters et al. (2014). ^{c} This option was written into the code for CO2SYSMATLAB v2.0.5 and other versions, but commented out and therefore not directly usable. It is available in CO2SYSMATLAB v3.2.0.
Other internal settings are consistent across all cases (Table 3). These three cases have been present since the original CO2SYS for MSDOS (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:

calculated on the pH scale reported in the literature as a function of temperature and salinity at zero inwater pressure;

converted to the seawater pH scale (Appendix A);

corrected to the in situ pressure;

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_{2} 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 asis 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 pressuresensitive, 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 subzero 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 outputcondition results from PyCO2SYS then represent the true state of the sample in situ in its environment. The inputcondition results are of less environmental interest but may be useful for qualitycontrol 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 outputcondition 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_{2}(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_{2}(aq)], which would not be valid because these are all directly proportional to each other at a given temperature, salinity, and atmospheric pressure.
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.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
where ${A}_{\mathrm{T}}^{\prime}=\mathrm{d}{A}_{\mathrm{T}}/\mathrm{dpH}$ and
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 openocean 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 CO2SYSMATLAB v2.0.5, the solvers continue to iterate and update all values until the change in every element of the array satisfies the ΔpH tolerance threshold (10^{−4} in CO2SYSMATLAB, 10^{−8} in PyCO2SYS). This means that a given set of arguments could return slightly different results depending on what data appear in the other, supposedly independent, elements of the argument arrays. Although negligible for all practical purposes, these differences are detectable in code validation exercises. In PyCO2SYS (and in CO2SYSMATLAB 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 CO2SYSMATLAB, 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 CO2SYSMATLAB functions, carried through from the original MSDOS 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 CO2SYSMATLAB v3.2.0 (Sharp et al., 2020) and CO2SYSExcel v3 (Pierrot et al., 2021) assume that the userselected pH scale is total and thus apply the totaltofree scale conversion (Appendix A) regardless of what the userselected 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^{−10} µmol kg^{−1} relative to A_{T} on the order of 2000 µmol kg^{−1}) and then only when the userselected 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, CO2SYSMATLAB v3.2.0, and CO2SYSExcel v3, the correct conversion factor is used based on the userselected 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 openocean 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.
Following Munhoven (2013) and as also implemented elsewhere (e.g. Orr and Epitalon, 2015), PyCO2SYS and CO2SYSMATLAB 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_{2}(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_{2}(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
This can be rearranged into a thirdorder polynomial in h:
with the following.
Following a scheme equivalent to Munhoven (2013), the initial h value is determined by
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_{0} of 10^{−3} mol kg^{−1}, corresponding to a pH of 3, is therefore used in this case. Otherwise, h_{min} in Eq. (9) is found following Munhoven (2013).
When A_{T} is positive, the squarerooted 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_{0} of 10^{−7} mol kg^{−1} 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
This can be rearranged into a secondorder polynomial in h:
with the following.
The initial h value is estimated following
When b<A_{T}, the squarerooted term ${g}_{\mathrm{1}}^{\mathrm{2}}\mathrm{4}{g}_{\mathrm{0}}{g}_{\mathrm{2}}$ is always positive, and thus h_{0}(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
This can be rearranged into a secondorder polynomial in h:
with the following.
The initial h value is estimated following
When $\mathrm{2}c+{T}_{\mathrm{B}}<{A}_{\mathrm{T}}$, the squarerooted term ${g}_{\mathrm{1}}^{\mathrm{2}}\mathrm{4}{g}_{\mathrm{0}}{g}_{\mathrm{2}}$ is always positive, and thus h_{0}(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 CO2SYSMATLAB 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 GEOSECScompatible results. If values are provided, then they will be included in the alkalinity equation just as for the nonGEOSECS 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 userdefined additional components within PyCO2SYS; the user must ensure that they are already suitable for the conditions being analysed and on the userindicated pH scale.
3.3.2 Gas constant
Previous versions of CO2SYS used an old value for the universal gas constant (R) of 8.31451 J mol^{−1} K^{−1}. PyCO2SYS uses the 2018 CODATArecommended value by default instead (i.e. 8.314462618 J mol^{−1} K^{−1}), consistent with CO2SYSMATLAB v3.2.0 and CO2SYSExcel 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^{−4} %), as well as on the pressure corrections for the equilibrium constants (less than 10^{−3} % at 5000 dbar). It is detectable in comparisons with other versions of CO2SYS, but it is of no practical consequence.
3.3.3 Substrate : inhibitor ratio
Like CO2SYSMATLAB 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_{3} 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 finitecentraldifference 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 CO2SYSMATLAB and CO2SYSExcel. 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 nonnegligible 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 Nosolve 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_{2}(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 outputcondition values are not calculated.
Seawater ${p}_{{\mathrm{CO}}_{\mathrm{2}}}$, ${f}_{{\mathrm{CO}}_{\mathrm{2}}}$, ${x}_{{\mathrm{CO}}_{\mathrm{2}}}$, and [CO_{2}(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 outputcondition 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 ndarray
s (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 CO2SYSMATLAB since its v1.1 (van Heuven et al., 2011). However, all multidimensional arrays in CO2SYSMATLAB are flattened into onedimensional vectors and returned in the results in that same format.
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 MSDOS 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 CO2SYSMATLAB, its implementation of uncertainty propagation differs.
PyCO2SYS evaluates the derivatives using a finiteforwarddifference 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 centraldifference 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.
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 nonlinear, 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 openocean conditions and selected an appropriate value between these extremes (e.g. Fig. 3). The full list of Δa values is provided in Table H1.
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
where σ is the uncertainty as a standard deviation (thus, σ^{2} is a variance). However, Eq. (24) is only valid if the uncertainties in all arguments are independent from each other. Propagation of covarying 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 covarying argument uncertainties through to any result (JCGM, 2008).
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.
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 Roundrobin test
In a “roundrobin” test, we first determine all of the core carbonate system parameters from one pair and then solve the system again using every possible pair of determined parameters. Under typical seawater conditions, we find the same results for every parameter pair to within better than the tolerance of the iterative pH solvers (i.e. 10^{−8} in pH). The maximum absolute difference in each parameter across all possible input pair combinations is acceptably small (Table 4).
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^{−12} % for the Egleston et al. (2010) buffers, 10^{−9} % for ψ and Q, and 10^{−7} % for the Revelle factor. The Revelle factor is less wellmatched because its explicit value is computed using a finitedifference scheme (for consistency with CO2SYSMATLAB), 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 forwarddifference 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^{4}) as a standard deviation is either within 3 % of the directly calculated value if the latter is nonzero or negligibly small if it is zero (absolute value less than 10^{−10}). The 3 % cutoff is relatively high because of the relatively small number of simulations; the cutoff can be reduced if a greater number of simulations is used, but then the computation time for the test suite becomes impractically long.
4.2 Comparison with other CO2SYS software
We used CO2SYSMATLAB v2.0.5 (Orr et al., 2018) as the main alternative software to compare our results with. PyCO2SYS was originally created as an ascloseaspossible Python translation of this particular version, so any differences in the results should be both understood and intentional. Its predecessor, CO2SYSMATLAB 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. CO2SYSMATLAB 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 CO2SYSMATLAB 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 CO2SYSMATLAB 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 CO2SYSMATLAB v3.2.0 (Sharp et al., 2020), which does include all these options.
4.2.1 Temperature–salinity–pressure parameterisations
All equilibrium constants and total salt contents, calculated from salinity, temperature, and pressure, are virtually identical (absolute tolerance 10^{−12}, relative tolerance 10^{−16}, in pK^{*} values or in µmol kg^{−1}) to those in both CO2SYSMATLAB 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^{5} 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 CO2SYSMATLAB v2.0.5, i.e. if the following points are true, then the differences between PyCO2SYS and CO2SYSMATLAB calculations are virtually zero (no greater than 10^{−10} %, excluding the Revelle factor as noted above):

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

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

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

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

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

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 CO2SYSMATLAB (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 CO2SYSMATLAB v2.0.5 are up to the order of 10^{−5} %: greater, but still negligible for all practical purposes.
Fixing the pH scale conversion simplification too (Sect. 3.1.3) makes no difference to calculations for which the userdefined input pH scale is total but causes discrepancies between PyCO2SYS and CO2SYSMATLAB v2.0.5 of up to 50 % in the “free” hydrogen ion content and 10^{−2} % in other results when other input pH scale options are selected. The differences are amplified at low pH, as the assumptions of the pH scale conversion simplification do not hold (Sect. 3.1.3).
Repeating the exercise above for CO2SYSMATLAB 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 CO2SYSMATLAB v3.2.0, differences in calculated values are still mostly less than 10^{−10} % and with one exception all less than 10^{−6} %. The exception, a difference still less than 10^{−3} %, is for the aqueous CO_{2} content under a limited set of input conditions and only with the new known parameter pair combinations added since CO2SYSMATLAB v2.0.5. It arises because there are several different ways to calculate [CO_{2}(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_{2} 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 floatingpoint precision. PyCO2SYS and CO2SYSMATLAB do not always use the same approach to calculate [CO_{2}(aq)] in each situation (this also varies between CO2SYSMATLAB 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_{2}(aq)] should always be calculated from the other three by addition or by difference as appropriate.
4.2.3 Uncertainty propagation comparisons
PyCO2SYS reproduces all the derivatives reported by Orr et al. (2018) in their Tables 2 and 3 to within 10^{−3} % under the same input conditions and all the propagated uncertainties reported by Orr et al. (2018) in their Table 4 to within 10^{−4} %. We consider all these differences to be negligible.
Across all combinations of optional parameters, mean uncertainties in 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_{2}(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 CO2SYSMATLAB v3.2.0 under the same input conditions. Greater differences in uncertainties calculated under output conditions arise because CO2SYSMATLAB does not propagate the uncertainties from inputcondition equilibrium constants through to outputcondition results.
4.3 Simulated seawater titration
PyCO2SYS can be used to reproduce the closedcell 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 sodetermined 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^{−1} of total phosphate included, is reproduced perfectly by PyCO2SYS with the exception of three values at different titrant masses.

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

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

1.25 g: pH either 5.549957 (Dickson) or 5.549951 (PyCO2SYS)
The other 48 data points in this titration agree perfectly. The noted discrepancies occur in nonconsecutive 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.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.
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 initialestimate equations become invalid and default pH values are used instead (Fig. 5). The number of iterations required to fall beneath the solver's tolerance threshold (10^{−8} in pH) is also reduced compared with the original approach of always using an initial pH of 8. Indeed, for typical ocean conditions we find that the iterative solver often does not alter the initial estimate at all. Suitable starting points for the iterative solvers are clearly being found.
5.2 Parameter pairs with multiple solutions
It is not strictly true that the marine carbonate system can always be solved from any pair of its parameters. Some combinations have multiple solutions. For example, both the 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 WolfGladrow, 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.
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 lowerpH root corresponds to typical seawater: a relatively highT_{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 higherpH root corresponds to a lowT_{C} system, wherein virtually all of T_{C} is in the form of carbonate ion, and A_{T} is dominated by noncarbonate species (Fig. 7).
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 initialestimate 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 lowerpH true root, which is appropriate for seawater. The solver will thus more robustly find the correct root each time.
In typical openocean 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 WolfGladrow, 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_{2}(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 secondorder 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 WolfGladrow (2001), is to take the higherpH 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 rootselection strategy for this parameter pair combination in more detail.
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 “goldstandard” equilibratorbased methodology. Second, the concept of ${p}_{{\mathrm{CO}}_{\mathrm{2}}}$ has utility only in the context of air–sea CO_{2} 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 pressurecorrected ${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 CO2SYSMATLAB v3.2.0 across a few different tasks for reference purposes. We ran CO2SYSMATLAB in both MATLAB itself (expensive, proprietary software) and GNU Octave, a free and opensource MATLAB clone.
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 CO2SYSMATLAB, but it is generally somewhat slower. However, the difference is negligible in practice for relatively small datasets (up to about 10^{5} 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.
CO2SYSMATLAB 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 outputcondition results if they are explicitly requested (Sect. 2.3), whereas CO2SYSMATLAB always calculates its results at both input and output conditions.
The results in Table 5 show that CO2SYSMATLAB 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 CO2SYSMATLAB running in MATLAB and by 3.8 if both the input and outputcondition 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 “justintime” 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 CO2SYSfamily 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 CO2SYSfamily 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 systemsolving 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.
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 WolfGladrow (2001) and Velo et al. (2010).
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.
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:
or alternatively and equivalently
The equations above are used in the same way to convert K^{*} values between pH scales.
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; WolfGladrow et al., 2007; Sharp and Byrne, 2020).
Equations for all the individual alkalinity components (A_{C}, A_{B}, etc.) are given in the subsequent sections in terms of pHindependent total substance contents (T_{C}, T_{B}, etc.) and [H^{+}].
B2 Water
B3 Carbonic acid
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]$.
Undissociated H_{2}CO_{3} is considered negligible and thus not explicitly modelled, but rather implicitly included as part of the [CO_{2}(aq)] term (Zeebe and WolfGladrow, 2001).
B4 Boric acid
B5 Phosphoric acid
B6 Orthosilicic acid
Further deprotonation of ${\mathrm{H}}_{\mathrm{3}}{\mathrm{SiO}}_{\mathrm{4}}^{}$ is considered negligible and thus not modelled.
B7 Ammonium
B8 Sulfide
Further deprotonation of HS^{−} is considered negligible and thus not modelled (Schoonen and Barnes, 1988).
B9 Sulfate
Undissociated H_{2}SO_{4} is considered negligible and thus not modelled.
B10 Fluoride
B11 Arbitrary additional components
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 userprovided ${K}_{\mathit{\alpha}}^{*}$ and ${K}_{\mathit{\beta}}^{*}$ values, with a zero level of protons corresponding to a pK^{*} of 4.5 (WolfGladrow 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 userdefined dissociation constants at the given conditions because one cannot convert arbitrary dissociation constants to their alkalinityrelevant values. Interpretations of results when arbitrary components are supplied to PyCO2SYS with pK^{*} values close to 4.5 should consider this nuance.
Here, we lay out all the equations that are used to convert between different carbonate system parameters in PyCO2SYS. These follow longestablished approaches from the literature (Zeebe and WolfGladrow, 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 userspecified pH scale, i.e. consistent with the pH values, pH and [H^{+}] are interconverted in the equations throughout this section using
regardless of which pH scale is being used.
C1.2 Known ${p}_{{\mathrm{CO}}_{\mathrm{2}}}$, ${x}_{{\mathrm{CO}}_{\mathrm{2}}}$ or [CO_{2}(aq)]
If one of ${p}_{{\mathrm{CO}}_{\mathrm{2}}}$, ${x}_{{\mathrm{CO}}_{\mathrm{2}}}$, or [CO_{2}(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}}}$,
where G is the fugacity factor (Table 2), typically near 0.997.
For known ${x}_{{\mathrm{CO}}_{\mathrm{2}}}$,
where P_{v} is the humidity correction (Table 3):
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_{2}(aq)],
where ${K}_{\mathrm{0}}^{*}$ is the solubility factor for CO_{2} (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_{2}(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_{2}(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.
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}.
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
where
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
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
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
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
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
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
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
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
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).
Calcite and aragonite saturation states (Ω) are calculated from the definition
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.
Note that in Eq. (D2), the [H^{+}] term is always calculated on the free pH scale of Eq. (A1).
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 \mathrm{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 \mathrm{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.
Here, e is Euler's number (2.71828…).
For the saturationstate buffers ${\mathit{\omega}}_{{T}_{\mathrm{C}}}$ and ${\mathit{\omega}}_{{A}_{\mathrm{T}}}$ we also evaluate

$(\partial \mathrm{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.
The approach taken here avoids AD evaluations over the iterative solvers because, while possible, that is computationally slower than over noniterative 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).
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.
Finally, the “released CO_{2} : precipitated carbonate ratio” (ψ) of Frankignoulle et al. (1994) is calculated following Humphreys et al. (2018).
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
This can be rearranged into a thirdorder polynomial in h:
with the following.
The initial h value is determined by
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_{0} of 10^{−3} mol kg^{−1}, corresponding to a pH of 3, is therefore used for that case (e.g. after the alkalinity end point in an acidimetric titration). The maximum possible 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_{0} of 10^{−7} mol kg^{−1} (pH 7) is used when ${g}_{\mathrm{2}}^{\mathrm{2}}\mathrm{3}{g}_{\mathrm{1}}\le \mathrm{0}$ in Eq. (F6) (Munhoven, 2013).
Older versions of CO2SYSMATLAB, 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 CO2SYSMATLAB v3.2.0 and CO2SYSExcel 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 floatingpoint 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 finitedifference derivatives. The key errors in the original CO2SYSMATLAB implementation of the finitedifference approach are the following.

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}).

Under output conditions, the Peng correction is not included in the evaluation of the Revelle factor (Sect. 2.2).
The lower accuracy of the finitedifference method relative to automatic differentiation, particularly given the relatively large ΔT_{C} used in the original finitedifference implementation (i.e. 1 µmol kg^{−1}), explains the differences between the two approaches that remain after the errors above have been corrected.
The computational speed tests described in Sect. 5.4 were run on an HP Spectre x360 laptop with an Intel Core i78565U 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 CO2SYSMATLAB v3.2.0.
The GNU Octave tests were run using GNU Octave v6.3.0 via its commandline interface and CO2SYSMATLAB v3.2.0.
The current version of PyCO2SYS is freely available from its GitHub repository at https://github.com/mvdh7/PyCO2SYS (last access: 23 December 2021) under the GNU General Public License v3. Installation is recommended from the Python Package Index (PyPI) via pip
, and documentation is available online (https://PyCO2SYS.readthedocs.io, last access: 23 December 2021). The exact version of PyCO2SYS used to produce the results discussed in this paper (v1.8.0), including input data and scripts to run the model and perform all validation tests described here, is archived on Zenodo (https://doi.org/10.5281/zenodo.5602840, Humphreys et al., 2021).
No data sets were used in this article.
MPH was responsible for conceptualisation, methodology, software, validation, writing the original draft, and visualisation. ERL was responsible for software and writing (review and editing). JDS was responsible for software, validation, and writing (review and editing). DP was responsible for software and writing (review and editing).
The contact author has declared that neither they nor their coauthors have any competing interests.
Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
We thank Doug Wallace for providing useful comments on this paper, and we acknowledge his important role in the creation of the original CO2SYS software. We further acknowledge the developers of all subsequent versions of CO2SYS upon whose work PyCO2SYS was built. We thank Luke Gregor, Daniel Sandborn, and Abigail Schiller for code contributions including extending the range of data types with which PyCO2SYS can be used. We are grateful to Guy Munhoven and James Orr for their detailed and constructive reviews.
This paper was edited by Paul Halloran and reviewed by James Orr and Guy Munhoven.
Abril, G., Bouillon, S., Darchambeau, F., Teodoru, C. R., Marwick, T. R., Tamooh, F., Ochieng Omengo, F., Geeraert, N., Deirmendjian, L., Polsenaere, P., and Borges, A. V.: Technical Note: Large overestimation of pCO_{2} calculated from pH and alkalinity in acidic, organicrich freshwaters, Biogeosciences, 12, 67–78, https://doi.org/10.5194/bg12672015, 2015. a
Álvarez, M., SanleónBartolomé, H., Tanhua, T., Mintrop, L., Luchetta, A., Cantoni, C., Schroeder, K., and Civitarese, G.: The CO_{2} system in the Mediterranean Sea: a basin wide perspective, Ocean Sci., 10, 69–92, https://doi.org/10.5194/os10692014, 2014. a
Bach, L. T.: Reconsidering the role of carbonate ion concentration in calcification by marine organisms, Biogeosciences, 12, 4939–4951, https://doi.org/10.5194/bg1249392015, 2015. a, b
Bakker, D. C. E., Pfeil, B., Landa, C. S., Metzl, N., O'Brien, K. M., Olsen, A., Smith, K., Cosca, C., Harasawa, S., Jones, S. D., Nakaoka, S., Nojiri, Y., Schuster, U., Steinhoff, T., Sweeney, C., Takahashi, T., Tilbrook, B., Wada, C., Wanninkhof, R., Alin, S. R., Balestrini, C. F., Barbero, L., Bates, N. R., Bianchi, A. A., Bonou, F., Boutin, J., Bozec, Y., Burger, E. F., Cai, W.J., Castle, R. D., Chen, L., Chierici, M., Currie, K., Evans, W., Featherstone, C., Feely, R. A., Fransson, A., Goyet, C., Greenwood, N., Gregor, L., Hankin, S., HardmanMountford, N. J., Harlay, J., Hauck, J., Hoppema, M., Humphreys, M. P., Hunt, C. W., Huss, B., Ibánhez, J. S. P., Johannessen, T., Keeling, R., Kitidis, V., Körtzinger, A., Kozyr, A., Krasakopoulou, E., Kuwata, A., Landschützer, P., Lauvset, S. K., Lefèvre, N., Lo Monaco, C., Manke, A., Mathis, J. T., Merlivat, L., Millero, F. J., Monteiro, P. M. S., Munro, D. R., Murata, A., Newberger, T., Omar, A. M., Ono, T., Paterson, K., Pearce, D., Pierrot, D., Robbins, L. L., Saito, S., Salisbury, J., Schlitzer, R., Schneider, B., Schweitzer, R., Sieger, R., Skjelvan, I., Sullivan, K. F., Sutherland, S. C., Sutton, A. J., Tadokoro, K., Telszewski, M., Tuma, M., van Heuven, S. M. A. C., Vandemark, D., Ward, B., Watson, A. J., and Xu, S.: A multidecade record of highquality fCO_{2} data in version 3 of the Surface Ocean CO2 Atlas (SOCAT), Earth Syst. Sci. Data, 8, 383–413, https://doi.org/10.5194/essd83832016, 2016. a
Bradbury, J., Frostig, R., Hawkins, P., Johnson, M. J., Leary, C., Maclaurin, D., and WandermanMilne, S.: JAX: composable transformations of Python+NumPy programs, GitHub, available at: http://github.com/google/jax (last access: 23 December 2021), 2018. a
Branson, O.: oscarbranson/cbsyst: beta, Zenodo [code], https://doi.org/10.5281/zenodo.1402261, 2018. a
Broecker, W. S., Takahashi, T., Simpson, H. J., and Peng, T.H.: Fate of Fossil Fuel Carbon Dioxide and the Global Carbon Budget, Science, 206, 409–418, https://doi.org/10.1126/science.206.4417.409, 1979. a
Cai, W.J. and Wang, Y.: The chemistry, fluxes, and sources of carbon dioxide in the estuarine waters of the Satilla and Altamaha Rivers, Georgia, Limnol. Oceanogr., 43, 657–668, https://doi.org/10.4319/lo.1998.43.4.0657, 1998. a
Cai, W.J., Huang, W.J., Luther, G. W., Pierrot, D., Li, M., Testa, J., Xue, M., Joesoef, A., Mann, R., Brodeur, J., Xu, Y.Y., Chen, B., Hussain, N., Waldbusser, G. G., Cornwell, J., and Kemp, W. M.: Redox reactions and weak buffering capacity lead to acidification in the Chesapeake Bay, Nat. Commun., 8, 369, https://doi.org/10.1038/s41467017004177, 2017. a
Cantrell, K. J., Serkiz, S. M., and Perdue, E. M.: Evaluation of acid neutralizing capacity data for solutions containing natural organic acids, Geochim. Cosmochim. Acta, 54, 1247–1254, https://doi.org/10.1016/00167037(90)90150J, 1990. a
Clarke, J. S., Achterberg, E. P., Connelly, D. P., Schuster, U., and Mowlem, M.: Developments in marine pCO_{2} measurement technology; towards sustained in situ observations, Trends Anal. Chem., 88, 53–61, https://doi.org/10.1016/j.trac.2016.12.008, 2017. a
Clegg, S. L. and Whitfield, M.: A chemical model of seawater including dissolved ammonia and the stoichiometric dissociation constant of ammonia in estuarine water and seawater from −2 to 40 ^{∘}C, Geochim. Cosmochim. Acta, 59, 2403–2421, https://doi.org/10.1016/00167037(95)001352, 1995. a
Culkin, F.: The major constituents of sea water, in: Chemical Oceanography, edited by: Riley, J. P. and Skirrow, G., Academic Press, London, UK, vol. 1, 121–161, 1965. a, b
Deffeyes, K. S.: Carbonate Equilibria: A Graphic and Algebraic Approach, Limnol. Oceanogr., 10, 412–426, https://doi.org/10.4319/lo.1965.10.3.0412, 1965. a
Dickson, A. and Millero, F.: A comparison of the equilibrium constants for the dissociation of carbonic acid in seawater media, DeepSea Res. Pt. A, 34, 1733–1743, https://doi.org/10.1016/01980149(87)900215, 1987. a, b, c
Dickson, A. G.: An exact definition of total alkalinity and a procedure for the estimation of alkalinity and total inorganic carbon from titration data, DeepSea Res. Pt. A, 28, 609–623, https://doi.org/10.1016/01980149(81)901217, 1981. a, b, c, d, e, f
Dickson, A. G.: Standard potential of the reaction: $\mathrm{AgCl}\left(\mathrm{s}\right)+\mathrm{0.5}{\mathrm{H}}_{\mathrm{2}}\left(\mathrm{g}\right)=\mathrm{Ag}\left(\mathrm{s}\right)+\mathrm{HCl}\left(\mathrm{aq}\right)$, and the standard acidity constant of the ion HSO${}_{\mathrm{4}}^{}$ in synthetic sea water from 273.15 to 318.15 K, J. Chem. Thermodyn., 22, 113–127, https://doi.org/10.1016/00219614(90)90074Z, 1990a. a
Dickson, A. G.: Thermodynamics of the dissociation of boric acid in synthetic seawater from 273.15 to 318.15 K, DeepSea Res. Pt. A, 37, 755–766, https://doi.org/10.1016/01980149(90)90004F, 1990b. a
Dickson, A. G. and Riley, J. P.: The estimation of acid dissociation constants in seawater media from potentiometric titrations with strong base. II. The dissociation of phosphoric acid, Mar. Chem., 7, 101–109, https://doi.org/10.1016/03044203(79)900021, 1979. a
Dickson, A. G., Sabine, C. L., and Christian, J. R. (Eds.): Guide to Best Practices for Ocean CO_{2} Measurements, PICES Special Publication 3, North Pacific Marine Science Organization, Sidney, BC, Canada, 2007. a
Dickson, A. G., Camões, M. F., Spitzer, P., Fisicaro, P., Stoica, D., Pawlowicz, R., and Feistel, R.: Metrological challenges for measurements of key climatological observables. Part 3: seawater pH, Metrologia, 53, R26–R39, https://doi.org/10.1088/00261394/53/1/R26, 2015. a
Doney, S. C., Fabry, V. J., Feely, R. A., and Kleypas, J. A.: Ocean Acidification: The Other CO_{2} Problem, Annu. Rev. Marine Sci., 1, 169–192, https://doi.org/10.1146/annurev.marine.010908.163834, 2009. a
Edmond, J. M. and Gieskes, J. M. T. M.: On the calculation of the degree of saturation of sea water with respect to calcium carbonate under in situ conditions, Geochim. Cosmochim. Acta, 34, 1261–1291, https://doi.org/10.1016/00167037(70)900414, 1970. a
Egleston, E. S., Sabine, C. L., and Morel, F. M. M.: Revelle revisited: Buffer factors that quantify the response of ocean chemistry to changes in DIC and alkalinity, Global Biogeochem. Cy., 24, GB1002, https://doi.org/10.1029/2008GB003407, 2010. a, b, c, d, e, f, g, h, i, j
Frankignoulle, M.: A complete set of buffer factors for acid/base CO_{2} system in seawater, J. Mar. Syst., 5, 111–118, https://doi.org/10.1016/09247963(94)900264, 1994. a
Frankignoulle, M., Canon, C., and Gattuso, J.P.: Marine calcification as a source of carbon dioxide: Positive feedback of increasing atmospheric CO_{2}, Limnol. Oceanogr., 39, 458–462, https://doi.org/10.4319/lo.1994.39.2.0458, 1994. a, b, c
Friedlingstein, P., O'Sullivan, M., Jones, M. W., Andrew, R. M., Hauck, J., Olsen, A., Peters, G. P., Peters, W., Pongratz, J., Sitch, S., Le Quéré, C., Canadell, J. G., Ciais, P., Jackson, R. B., Alin, S., Aragão, L. E. O. C., Arneth, A., Arora, V., Bates, N. R., Becker, M., BenoitCattin, A., Bittig, H. C., Bopp, L., Bultan, S., Chandra, N., Chevallier, F., Chini, L. P., Evans, W., Florentie, L., Forster, P. M., Gasser, T., Gehlen, M., Gilfillan, D., Gkritzalis, T., Gregor, L., Gruber, N., Harris, I., Hartung, K., Haverd, V., Houghton, R. A., Ilyina, T., Jain, A. K., Joetzjer, E., Kadono, K., Kato, E., Kitidis, V., Korsbakken, J. I., Landschützer, P., Lefèvre, N., Lenton, A., Lienert, S., Liu, Z., Lombardozzi, D., Marland, G., Metzl, N., Munro, D. R., Nabel, J. E. M. S., Nakaoka, S.I., Niwa, Y., O'Brien, K., Ono, T., Palmer, P. I., Pierrot, D., Poulter, B., Resplandy, L., Robertson, E., Rödenbeck, C., Schwinger, J., Séférian, R., Skjelvan, I., Smith, A. J. P., Sutton, A. J., Tanhua, T., Tans, P. P., Tian, H., Tilbrook, B., van der Werf, G., Vuichard, N., Walker, A. P., Wanninkhof, R., Watson, A. J., Willis, D., Wiltshire, A. J., Yuan, W., Yue, X., and Zaehle, S.: Global Carbon Budget 2020, Earth Syst. Sci. Data, 12, 3269–3340, https://doi.org/10.5194/essd1232692020, 2020. a
Gattuso, J.P., Epitalon, J.M., Lavigne, H., Orr, J., Gentili, B., Hagens, M., Hofmann, A., Mueller, J.D., Proye, A., Rae, J., and Soetaert, K.: seacarb: Seawater Carbonate Chemistry, available at: https://CRAN.Rproject.org/package=seacarb, last access: 23 December 2021. a
Goyet, C. and Poisson, A.: New determination of carbonic acid dissociation constants in seawater as a function of temperature and salinity, DeepSea Res. Pt. A, 36, 1635–1654, https://doi.org/10.1016/01980149(89)900642, 1989. a
Hagens, M. and Middelburg, J. J.: Generalised expressions for the response of pH to changes in ocean chemistry, Geochim. Cosmochim. Acta, 187, 334–349, https://doi.org/10.1016/j.gca.2016.04.012, 2016. a
Hansson, I.: The Determination of Dissociation Constants of Carbonic Acid in Synthetic Sea Water in the Salinity Range of 20–40 ‰ and Temperature Range of 5–30 ^{∘}C, Acta Chem. Scand., 27, 931–944, https://doi.org/10.3891/acta.chem.scand.270931, 1973a. a, b
Hansson, I.: A new set of acidity constants for carbonic acid and boric acid in sea water, DeepSea Res., 20, 461–478, https://doi.org/10.1016/00117471(73)901009, 1973b. a, b
Harris, C. R., Millman, K. J., van der Walt, S. J., Gommers, R., Virtanen, P., Cournapeau, D., Wieser, E., Taylor, J., Berg, S., Smith, N. J., Kern, R., Picus, M., Hoyer, S., van Kerkwijk, M. H., Brett, M., Haldane, A., del Río, J. F., Wiebe, M., Peterson, P., GérardMarchant, P., Sheppard, K., Reddy, T., Weckesser, W., Abbasi, H., Gohlke, C., and Oliphant, T. E.: Array programming with NumPy, Nature, 585, 357–362, https://doi.org/10.1038/s4158602026492, 2020. a
Humphreys, M. P., Daniels, C. J., WolfGladrow, D. A., Tyrrell, T., and Achterberg, E. P.: On the influence of marine biogeochemical processes over CO_{2} exchange between the atmosphere and ocean, Mar. Chem., 199, 1–11, https://doi.org/10.1016/j.marchem.2017.12.006, 2018. a, b, c, d, e
Humphreys, M. P., Gregor, L., Pierrot, D., van Heuven, S. M. A. C., Lewis, E. R., and Wallace, D. W. R.: PyCO2SYS: marine carbonate system calculations in Python (v1.3.0), Zenodo [code], https://doi.org/10.5281/zenodo.3780139, 2020. a
Humphreys, M. P., Schiller, A. J., Sandborn, D., Gregor, L., Pierrot, D., van Heuven, S. M. A. C., Lewis, E. R., and Wallace, D. W. R.: PyCO2SYS: marine carbonate system calculations in Python (v1.8.0), Zenodo [code], https://doi.org/10.5281/zenodo.5602840, 2021. a, b, c, d
Ingle, S. E.: Solubility of calcite in the ocean, Mar. Chem., 3, 301–319, https://doi.org/10.1016/03044203(75)900109, 1975. a, b, c
Ingle, S. E., Culberson, C. H., Hawley, J. E., and Pytkowicz, R. M.: The solubility of calcite in seawater at atmospheric pressure and 35 ‰ salinity, Mar. Chem., 1, 295–307, https://doi.org/10.1016/03044203(73)900194, 1973. a
IUPAC: Compendium of Chemical Terminology, 2nd edn. (the “Gold Book”), Blackwell Scientific Publications, Oxford, UK, https://doi.org/10.1351/goldbook, 1997. a
JCGM: JCGM 100:2008 Evaluation of measurement data – Guide to the expression of uncertainty in measurement, Bureau International des Poids et Mesures, Sèvres, France, available at: https://www.bipm.org/utils/common/documents/jcgm/JCGM_100_2008_E.pdf (last access: 23 December 2021), 2008. a
Kester, D. R. and Pytkowicz, R. M.: Determination of the Apparent Dissociation Constants of Phosphoric Acid in Seawater, Limnol. Oceanogr., 12, 243–252, https://doi.org/10.4319/lo.1967.12.2.0243, 1967. a
Khoo, K. H., Ramette, R. W., Culberson, C. H., and Bates, R. G.: Determination of hydrogen ion concentrations in seawater from 5 to 40C: standard potentials at salinities from 20 to 45 per mille, Anal. Chem., 49, 29–34, https://doi.org/10.1021/ac50009a016, 1977. a
Kuliński, K., Schneider, B., Hammer, K., Machulik, U., and SchulzBull, D.: The influence of dissolved organic matter on the acid–base system of the Baltic Sea, J. Mar. Syst., 132, 106–115, https://doi.org/10.1016/j.jmarsys.2014.01.011, 2014. a
Lauvset, S. K., Lange, N., Tanhua, T., Bittig, H. C., Olsen, A., Kozyr, A., Álvarez, M., Becker, S., Brown, P. J., Carter, B. R., Cotrim da Cunha, L., Feely, R. A., van Heuven, S., Hoppema, M., Ishii, M., Jeansson, E., Jutterström, S., Jones, S. D., Karlsen, M. K., Lo Monaco, C., Michaelis, P., Murata, A., Pérez, F. F., Pfeil, B., Schirnick, C., Steinfeldt, R., Suzuki, T., Tilbrook, B., Velo, A., Wanninkhof, R., Woosley, R. J., and Key, R. M.: An updated version of the global interior ocean biogeochemical data product, GLODAPv2.2021, Earth Syst. Sci. Data, 13, 5565–5589, https://doi.org/10.5194/essd1355652021, 2021. a
Lee, K., Kim, T.W., Byrne, R. H., Millero, F. J., Feely, R. A., and Liu, Y.M.: The universal ratio of boron to chlorinity for the North Pacific and North Atlantic oceans, Geochim. Cosmochim. Acta, 74, 1801–1811, https://doi.org/10.1016/j.gca.2009.12.027, 2010. a
Lewis, E. and Wallace, D. W. R.: Program Developed for CO_{2} System Calculations, ORNL/CDIAC105, Carbon Dioxide Information Analysis Center, Oak Ridge National Laboratory, U.S. Department of Energy, Oak Ridge, TN, USA, https://doi.org/10.15485/1464255, 1998. a, b, c, d, e, f
Li, Y.H., Takahashi, T., and Broecker, W. S.: Degree of saturation of CaCO_{3} in the oceans, J. Geophys. Res., 74, 5507–5525, https://doi.org/10.1029/JC074i023p05507, 1969. a
Lueker, T. J., Dickson, A. G., and Keeling, C. D.: Ocean pCO_{2} calculated from dissolved inorganic carbon, alkalinity, and equations for K_{1} and K_{2}: validation based on laboratory measurements of CO_{2} in gas and seawater at equilibrium, Mar. Chem., 70, 105–119, https://doi.org/10.1016/S03044203(00)000220, 2000. a
Maclaurin, D.: Autograd: Automatic Differentiation for Python, in: Modeling, Inference and Optimization with Composable Differentiable Procedures, PhD thesis, Harvard University, Cambridge, MA, USA, 41–57, available at: https://dougalmaclaurin.com/phdthesis.pdf (last access: 23 December 2021), 2016. a
Mehrbach, C., Culberson, C. H., Hawley, J. E., and Pytkowicz, R. M.: Measurement of the Apparent Dissociation Constants of Carbonic Acid in Seawater at Atmospheric Pressure, Limnol. Oceanogr., 18, 897–907, https://doi.org/10.4319/lo.1973.18.6.0897, 1973. a, b, c, d
Millero, F. J.: The thermodynamics of the carbonate system in seawater, Geochim. Cosmochim. Acta, 43, 1651–1661, https://doi.org/10.1016/00167037(79)901844, 1979. a, b, c, d
Millero, F. J.: Influence of pressure on chemical processes in the sea, in: Chemical Oceanography, edited by: Riley, J. P. and Chester, R., Academic Press, 1983. a, b, c, d, e, f
Millero, F. J.: Thermodynamics of the carbon dioxide system in the oceans, Geochim. Cosmochim. Acta, 59, 661–677, https://doi.org/10.1016/00167037(94)00354O, 1995. a, b, c, d, e, f, g, h, i, j
Millero, F. J.: The Carbonate System in Marine Environments, in: Chemical Processes in Marine Environments, edited by: Gianguzza, A., Pelizetti, E., and Sammartano, S., Environmental Science, Springer, Berlin, Heidelberg, 9–41, https://doi.org/10.1007/9783662042076_2, 2000. a
Millero, F. J.: Carbonate constants for estuarine waters, Mar. Freshw. Res., 61, 139–142, https://doi.org/10.1071/MF09254, 2010. a
Millero, F. J., Pierrot, D., Lee, K., Wanninkhof, R., Feely, R., Sabine, C. L., Key, R. M., and Takahashi, T.: Dissociation constants for carbonic acid determined from field measurements, DeepSea Res. Pt. I, 49, 1705–1723, https://doi.org/10.1016/S09670637(02)000936, 2002. a
Millero, F. J., Graham, T. B., Huang, F., BustosSerrano, H., and Pierrot, D.: Dissociation constants of carbonic acid in seawater as a function of salinity and temperature, Mar. Chem., 100, 80–94, https://doi.org/10.1016/j.marchem.2005.12.001, 2006. a
Millero, F. J., Feistel, R., Wright, D. G., and McDougall, T. J.: The composition of Standard Seawater and the definition of the ReferenceComposition Salinity Scale, DeepSea Res. Pt. I, 55, 50–72, https://doi.org/10.1016/j.dsr.2007.10.001, 2008. a
Mojica Prieto, F. J. and Millero, F. J.: The values of pK_{1} + pK_{2} for the dissociation of carbonic acid in seawater, Geochim. Cosmochim. Acta, 66, 2529–2540, https://doi.org/10.1016/S00167037(02)008554, 2002. a
Morris, A. W. and Riley, J. P.: The bromide/chlorinity and sulphate/chlorinity ratio in sea water, DeepSea Res., 13, 699–705, https://doi.org/10.1016/00117471(66)906012, 1966. a, b
Muller, F. L. L. and Bleie, B.: Estimating the organic acid contribution to coastal seawater alkalinity by potentiometric titrations in a closed cell, Anal. Chim. Acta, 619, 183–191, https://doi.org/10.1016/j.aca.2008.05.018, 2008. a
Munhoven, G.: Mathematics of the total alkalinity–pH equation – pathway to robust and universal solution algorithms: the SolveSAPHE package v1.0.1, Geosci. Model Dev., 6, 1367–1388, https://doi.org/10.5194/gmd613672013, 2013. a, b, c, d, e, f, g, h, i
Munhoven, G.: SolveSAPHEr2 (v2.0.1): revisiting and extending the Solver Suite for AlkalinityPH Equations for usage with CO_{2}, HCO${}_{\mathrm{3}}^{\mathrm{}}$ or CO${}_{\mathrm{3}}^{\mathrm{2}\mathrm{}}$ input data, Geosci. Model Dev., 14, 4225–4240, https://doi.org/10.5194/gmd1442252021, 2021. a, b, c
Orr, J. C.: Recent and future changes in ocean carbonate chemistry, in: Ocean Acidification, edited by: Gattuso, J. P. and Hansson, L., Oxford University Press, 41–66, https://doi.org/10.1093/oso/9780199591091.003.0008, 2011. a
Orr, J. C. and Epitalon, J.M.: Improved routines to model the ocean carbonate system: mocsy 2.0, Geosci. Model Dev., 8, 485–499, https://doi.org/10.5194/gmd84852015, 2015. a, b, c, d, e
Orr, J. C., Epitalon, J.M., and Gattuso, J.P.: Comparison of ten packages that compute ocean carbonate chemistry, Biogeosciences, 12, 1483–1510, https://doi.org/10.5194/bg1214832015, 2015. a, b, c, d, e
Orr, J. C., Najjar, R. G., Aumont, O., Bopp, L., Bullister, J. L., Danabasoglu, G., Doney, S. C., Dunne, J. P., Dutay, J.C., Graven, H., Griffies, S. M., John, J. G., Joos, F., Levin, I., Lindsay, K., Matear, R. J., McKinley, G. A., Mouchet, A., Oschlies, A., Romanou, A., Schlitzer, R., Tagliabue, A., Tanhua, T., and Yool, A.: Biogeochemical protocols and diagnostics for the CMIP6 Ocean Model Intercomparison Project (OMIP), Geosci. Model Dev., 10, 2169–2199, https://doi.org/10.5194/gmd1021692017, 2017. a
Orr, J. C., Epitalon, J.M., Dickson, A. G., and Gattuso, J.P.: Routine uncertainty propagation for the marine carbon dioxide system, Mar. Chem., 207, 84–107, https://doi.org/10.1016/j.marchem.2018.10.006, 2018. a, b, c, d, e, f, g, h, i, j, k, l, m, n
Park, P. K.: Oceanic CO_{2} system: an evaluation of ten methods of investigation, Limnol. Oceanogr., 14, 179–186, https://doi.org/10.4319/lo.1969.14.2.0179, 1969. a
Peng, T.H., Takahashi, T., Broecker, W. S., and Olafsson, J.: Seasonal variability of carbon dioxide, nutrients and oxygen in the northern North Atlantic surface water: observations and a model, Tellus B, 39, 439–458, https://doi.org/10.3402/tellusb.v39i5.15361, 1987. a, b
Perez, F. F. and Fraga, F.: Association constant of fluoride and hydrogen ions in seawater, Mar. Chem., 21, 161–168, https://doi.org/10.1016/03044203(87)900363, 1987. a
Pierrot, D., Lewis, E., and Wallace, D. W. R.: MS Excel Program Developed for CO_{2} System Calculations, available at: https://cdiac.essdive.lbl.gov/ftp/co2sys/CO2SYS_calc_XLS_v2.1/ (last access: 23 December 2021), 2006. a
Pierrot, D., Epitalon, J.M., Orr, J. C., Lewis, E. R., and Wallace, D.: MS Excel program developed for CO_{2} system calculations – version 3.0, available at: https://github.com/dpierrot/co2sys_xl, last access: 23 December 2021. a, b
Raimondi, L., Matthews, J. B. R., Atamanchuk, D., AzetsuScott, K., and Wallace, D. W. R.: The internal consistency of the marine carbon dioxide system for high latitude shipboard and in situ monitoring, Mar. Chem., 213, 49–70, https://doi.org/10.1016/j.marchem.2019.03.001, 2019. a
Revelle, R. and Suess, H. E.: Carbon Dioxide Exchange Between Atmosphere and Ocean and the Question of an Increase of Atmospheric CO_{2} during the Past Decades, Tellus, 9, 18–27, https://doi.org/10.3402/tellusa.v9i1.9075, 1957. a
Richier, S., Achterberg, E. P., Humphreys, M. P., Poulton, A. J., Suggett, D. J., Tyrrell, T., and Moore, C. M.: Geographical CO_{2} sensitivity of phytoplankton correlates with ocean buffer capacity, Glob. Change Biol., 24, 4438–4452, https://doi.org/10.1111/gcb.14324, 2018. a
Riley, J. P.: The occurrence of anomalously high fluoride concentrations in the North Atlantic, DeepSea Res., 12, 219–220, https://doi.org/10.1016/00117471(65)900276, 1965. a, b
Riley, J. P. and Tongudai, M.: The major cation/chlorinity ratios in sea water, Chem. Geol., 2, 263–269, https://doi.org/10.1016/00092541(67)900265, 1967. a
Roy, R. N., Roy, L. N., Vogel, K. M., PorterMoore, C., Pearson, T., Good, C. E., Millero, F. J., and Campbell, D. M.: The dissociation constants of carbonic acid in seawater at salinities 5 to 45 and temperatures 0 to 45 ^{∘}C, Mar. Chem., 44, 249–267, https://doi.org/10.1016/03044203(93)902075, 1993. a
Schockman, K. M. and Byrne, R. H.: Spectrophotometric Determination of the Bicarbonate Dissociation Constant in Seawater, Geochim. Cosmochim. Acta, 300, 231–245, https://doi.org/10.1016/j.gca.2021.02.008, 2021. a, b, c
Schoonen, M. A. A. and Barnes, H. L.: An approximation of the second dissociation constant for H_{2}S, Geochim. Cosmochim. Acta, 52, 649–654, https://doi.org/10.1016/00167037(88)903262, 1988. a
Sharp, J. D. and Byrne, R. H.: Carbonate ion concentrations in seawater: Spectrophotometric determination at ambient temperatures and evaluation of propagated calculation uncertainties, Mar. Chem., 209, 70–80, https://doi.org/10.1016/j.marchem.2018.12.001, 2019. a
Sharp, J. D. and Byrne, R. H.: Interpreting measurements of total alkalinity in marine and estuarine waters in the presence of protonbinding organic matter, DeepSea Res. Pt. I, 165, 103 338, https://doi.org/10.1016/j.dsr.2020.103338, 2020. a, b
Sharp, J. D., Pierrot, D., Humphreys, M. P., Epitalon, J.M., Orr, J. C., Lewis, E. R., and Wallace, D. W.: CO2SYSv3 for MATLAB, Zenodo [code], https://doi.org/10.5281/zenodo.3952803, 2020. a, b, c, d, e
Sillén, L. G., Martell, A. E., and Bjerrum, J.: Stability constants of metalion complexes, special publication, 17th edn., Chemical Society, London, UK, 1964. a
Sulpis, O., Lauvset, S. K., and Hagens, M.: Current estimates of K_{1}* and K_{2}* appear inconsistent with measured CO_{2} system parameters in cold oceanic regions, Ocean Sci., 16, 847–862, https://doi.org/10.5194/os168472020, 2020. a, b, c
Takahashi, T., Williams, R. T., and Bos, D. L.: Carbonate Chemistry, in: GEOSECS Pacific Expedition: Hydrographic Data, vol. 3, 77–105, National Science Foundation, Washington, D.C., 1982. a, b, c, d, e
Takahashi, T., Sutherland, S. C., Wanninkhof, R., Sweeney, C., Feely, R. A., Chipman, D. W., Hales, B., Friederich, G., Chavez, F., Sabine, C., Watson, A., Bakker, D. C., Schuster, U., Metzl, N., YoshikawaInoue, H., Ishii, M., Midorikawa, T., Nojiri, Y., Körtzinger, A., Steinhoff, T., Hoppema, M., Olafsson, J., Arnarson, T. S., Tilbrook, B., Johannessen, T., Olsen, A., Bellerby, R., Wong, C., Delille, B., Bates, N., and de Baar, H. J.: Climatological mean and decadal change in surface ocean pCO_{2}, and net sea–air CO_{2} flux over the global oceans, DeepSea Res. Pt. II, 56, 554–577, https://doi.org/10.1016/j.dsr2.2008.12.009, 2009. a
Turner, D. R., Achterberg, E. P., Chen, C.T. A., Clegg, S. L., Hatje, V., Maldonado, M. T., Sander, S. G., Berg, V. D., G, C. M., and Wells, M.: Toward a QualityControlled and Accessible Pitzer Model for Seawater and Related Systems, Front. Mar. Sci., 3, 139, https://doi.org/10.3389/fmars.2016.00139, 2016. a
Ulfsbo, A., Kuliński, K., Anderson, L. G., and Turner, D. R.: Modelling organic alkalinity in the Baltic Sea using a HumicPitzer approach, Mar. Chem., 168, 18–26, https://doi.org/10.1016/j.marchem.2014.10.013, 2015. a
Uppström, L. R.: The boron/chlorinity ratio of deepsea water from the Pacific Ocean, DeepSea Res., 21, 161–162, https://doi.org/10.1016/00117471(74)900746, 1974. a
van Heuven, S., Pierrot, D., Rae, J. W. B., Lewis, E., and Wallace, D. W. R.: CO_{2}SYS v 1.1, MATLAB program developed for CO_{2} system calculations, ORNL/CDIAC105b, Carbon Dioxide Information Analysis Center, Oak Ridge National Laboratory, U.S. Department of Energy, Oak Ridge, TN, USA, available at: https://cdiac.essdive.lbl.gov/ftp/co2sys/CO2SYS_calc_MATLAB_v1.1/ (last access: 23 December 2021), 2011. a, b, c
Velo, A., Pérez, F. F., Lin, X., Key, R. M., Tanhua, T., de la Paz, M., Olsen, A., van Heuven, S., Jutterström, S., and Ríos, A. F.: CARINA data synthesis project: pH data scale unification and cruise adjustments, Earth Syst. Sci. Data, 2, 133–155, https://doi.org/10.5194/essd21332010, 2010. a, b
Waters, J., Millero, F. J., and Woosley, R. J.: Corrigendum to “The free proton concentration scale for seawater pH”, [MARCHE: 149 (2013) 8–22], Mar. Chem., 165, 66–67, https://doi.org/10.1016/j.marchem.2014.07.004, 2014. a, b
Waters, J. F. and Millero, F. J.: The free proton concentration scale for seawater pH, Mar. Chem., 149, 8–22, https://doi.org/10.1016/j.marchem.2012.11.003, 2013. a, b, c
Weiss, R. F.: Carbon dioxide in water and seawater: the solubility of a nonideal gas, Mar. Chem., 2, 203–215, https://doi.org/10.1016/03044203(74)900152, 1974. a, b, c, d, e, f
Weiss, R. F. and Price, B. A.: Nitrous oxide solubility in water and seawater, Mar. Chem., 8, 347–359, https://doi.org/10.1016/03044203(80)900249, 1980. a, b
WolfGladrow, D. A., Zeebe, R. E., Klaas, C., Körtzinger, A., and Dickson, A. G.: Total alkalinity: The explicit conservative expression and its application to biogeochemical processes, Mar. Chem., 106, 287–300, https://doi.org/10.1016/j.marchem.2007.01.006, 2007. a, b
Xu, Y.Y., Pierrot, D., and Cai, W.J.: Ocean carbonate system computation for anoxic waters using an updated CO2SYS program, Mar. Chem., 195, 90–93, https://doi.org/10.1016/j.marchem.2017.07.002, 2017. a, b
Yao, W. and Millero, F. J.: The chemistry of the anoxic waters in the Framvaren Fjord, Norway, Aquat. Geochem., 1, 53–88, https://doi.org/10.1007/BF01025231, 1995. a, b, c
Zeebe, R. E. and WolfGladrow, D.: CO_{2} in Seawater: Equilibrium, Kinetics, Isotopes, Elsevier Oceanography Series 65, Elsevier B.V., Amsterdam, the Netherlands, 2001. a, b, c, d, e, f, g, h
 Abstract
 Introduction
 Methods inherited from CO2SYS
 New developments in PyCO2SYS
 Validation
 Discussion
 Appendix A: pH scales and conversions
 Appendix B: Total alkalinity and its components
 Appendix C: Solving the core marine carbonate system
 Appendix D: Other marine carbonate system variables
 Appendix E: Buffer factors with automatic differentiation
 Appendix F: Initial pH estimate when solving from A_{T} and T_{C}
 Appendix G: Revelle factor calculation errors in older versions of CO2SYSMATLAB
 Appendix H: Fixed Δa values for uncertainty analysis
 Appendix I: Setup for computational speed testing
 Code availability
 Data availability
 Author contributions
 Competing interests
 Disclaimer
 Acknowledgements
 Review statement
 References
ocean acidification), which may have adverse effects on marine ecosystems. Our Python package, PyCO2SYS, models the chemical reactions of CO_{2} in seawater, allowing us to quantify the corresponding changes in pH and related chemical properties.
 Abstract
 Introduction
 Methods inherited from CO2SYS
 New developments in PyCO2SYS
 Validation
 Discussion
 Appendix A: pH scales and conversions
 Appendix B: Total alkalinity and its components
 Appendix C: Solving the core marine carbonate system
 Appendix D: Other marine carbonate system variables
 Appendix E: Buffer factors with automatic differentiation
 Appendix F: Initial pH estimate when solving from A_{T} and T_{C}
 Appendix G: Revelle factor calculation errors in older versions of CO2SYSMATLAB
 Appendix H: Fixed Δa values for uncertainty analysis
 Appendix I: Setup for computational speed testing
 Code availability
 Data availability
 Author contributions
 Competing interests
 Disclaimer
 Acknowledgements
 Review statement
 References