# Routine uncertainty propagation for the marine carbon dioxide system

James C. Orr, Jean-Marie Epitalon, Andrew G Dickson, Jean-Pierre Gattuso

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

## To cite this version:

James C. Orr, Jean-Marie Epitalon, Andrew G Dickson, Jean-Pierre Gattuso. Routine uncertainty propagation for the marine carbon dioxide system. Marine Chemistry, 2018, 207, pp.84-107. 10.1016/j.marchem.2018.10.006. hal-01973728

HAL Id: hal-01973728

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

Submitted on 8 Jan 2019

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

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

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

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

# Routine uncertainty propagation for the marine carbon dioxide system

James C. Orr<sup>a, *,</sup>, Jean-Marie Epitalon<sup>b</sup>, Andrew G. Dickson<sup>c</sup>, Jean-Pierre Gattuso<sup>d, e</sup>

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

## ARTICLE INFO

### Keywords:

- Carbonate chemistry

- Uncertainty propagation

- Carbon dioxide

- pH

- Alkalinity

- Dissolved inorganic carbon

- Ocean acidification

## ABSTRACT

Pairs of marine carbonate system variables are often used to calculate others, but those results are seldom reported with estimates of uncertainties. Although the procedure to propagate these uncertainties is well known, it has not been offered in public packages that compute marine carbonate chemistry, fundamental tools that are relied on by the community. To remedy this shortcoming, four of these packages were expanded to calculate sensitivities of computed variables with respect to each input variable and to use those sensitivities along with user-specified estimates of input uncertainties (standard uncertainties) to propagate uncertainties of calculated variables (combined standard uncertainties). Sensitivities from these packages agree with one another and with analytical solutions to within 0.01%; similar agreement among packages was found for the combined standard uncertainties. One package was used to quantify how propagated uncertainties vary among computed variables, seawater conditions, and the chosen pair of carbonate system variables that is used as input. The relative contributions to propagated uncertainties from the standard uncertainties of the input pair of measurements and various other input data (equilibrium constants etc) were explored with a new type of diagram. These error-space diagrams illustrate that further improvement beyond today's state-of-the-art measurement uncertainties for the input pair would generally be ineffective at reducing the combined standard uncertainties because the contribution from the constants is larger. Likewise, using much more uncertain measurements of the input pair does not always substantially worsen combined standard uncertainty. The constants that contribute most to combined standard uncertainties are generally $K_1$ and $K_2$, as expected. Yet more of the propagated uncertainty in the computed saturation states of aragonite and calcite comes from their solubility products. Thus percent relative combined standard uncertainties for the saturation states are larger than for the carbonate ion concentration. Routine propagation of these uncertainties should become standard practice.

## 1. Introduction

Scientists commonly constrain the marine carbonate system by measuring two of its variables and using that pair to calculate the others. Less common is to report the propagated uncertainties of the calculated variables. Although the approach to propagate those uncertainties is well established (Dickson and Riley, 1978), it must be implemented by scientists, individually, because it is not an option in related public software packages (Orr et al., 2015). This technical gap may explain why these uncertainties are not routinely reported, even for efforts such as the ocean acidification community's long-term data compilation initiative (Nisumaa et al., 2010; Yang et al., 2016).

Among the limited studies that have reported propagated uncertainties, approaches and assumptions differ. Most have used Gaussian uncertainty

propagation (Dickson and Riley, 1978; Millero, 1995, 2007; McLaughlin et al., 2015; Sutton et al., 2016); others have used a Monte Carlo approach (Fassbender et al., 2016a; Williams et al., 2017) or a simpler but related method. That simpler approach was taken by Lauvset and Gruber (2014) who suggested that the Gaussian approach underestimates propagated uncertainties when its linear approximations (first order Taylor series) are applied to solve the series of nonlinear equations inherent in the marine $CO_2$ system. That concern, if generally true, would undermine previous work. Yet the approach is not the only source of differences between studies. Studies that have used the Gaussian approach have considered uncertainties both in the measured input pair of carbonate system variables and in the equilibrium constants needed to make these calculations; conversely, with the full or simplified Monte Carlo approach, only Williams et al. (2017) have accounted for uncertainties from the equilibrium constants.

<footnote><sup>a</sup> Laboratoire des Sciences du Climat et de l'Environnement, LSCE/IPSL, CEA-CNRS-UVSQ, Université Paris-Saclay, Gif-sur-Yvette, France</footnote>

<footnote><sup>b</sup> Geoscientific Programming Services, Fanjeaux, France</footnote>

<footnote><sup>c</sup> University of California, San Diego (UCSD), La Jolla, CA, USA</footnote>

<footnote><sup>d</sup> Sorbonne Université, CNRS, Laboratoire d'Océanographie de Villefranche, Villefranche-sur-mer, France</footnote>

<footnote><sup>e</sup> Institute for Sustainable Development and International Relations, Sciences Po, Paris, France</footnote>

<footnote>* Corresponding author.</footnote>

<footnote>E-mail address: james.orr@lsce.ipsl.fr (J.C. Orr).</footnote>
<table><caption>**Table 1**
Various estimates of standard uncertainties in input variables and constants.</caption><thead><tr><th rowspan="2">Variable<sup>a</sup></th><th colspan="8">Previous work</th><th colspan="2">This study</th><th rowspan="2">Units</th></tr><tr><th>DR78<sup>b</sup></th><th>M95<sup>b</sup></th><th>M06<sup>b</sup></th><th>M07<sup>b</sup></th><th>D10a<sup>c</sup></th><th>D10b<sup>c</sup></th><th>D10c<sup>c</sup></th><th>McL15<sup>d</sup></th><th>Random</th><th>Total</th></tr></thead><tbody><tr><td>A<sub>T</sub></td><td>9, 2</td><td>4</td><td>3</td><td>3</td><td>1.2</td><td>2–3</td><td>4–10</td><td>5</td><td>2</td><td>2</td><td>µmol kg<sup>-1</sup></td></tr><tr><td>C<sub>T</sub></td><td>10, 4</td><td>2</td><td>2</td><td>2</td><td>1</td><td>2–3</td><td>4–10</td><td>5</td><td>2</td><td>2</td><td>µmol kg<sup>-1</sup></td></tr><tr><td>pCO<sub>2</sub></td><td>7</td><td>2</td><td>2</td><td>2</td><td>1</td><td>2</td><td>5–10</td><td>6</td><td>2</td><td>2</td><td>µatm</td></tr><tr><td>pH</td><td>0.01</td><td>0.002</td><td>0.002</td><td>0.002</td><td>0.003</td><td>0.005</td><td>0.01–0.03</td><td>0.01</td><td>0.003</td><td>0.01</td><td></td></tr><tr><td>pK<sub>0</sub></td><td>0.002</td><td></td><td></td><td></td><td>0.002</td><td>0.002</td><td>0.002</td><td>0.002</td><td>0.002</td><td>0.002</td><td></td></tr><tr><td>pK<sub>1</sub></td><td>0.01</td><td>0.002</td><td>0.006</td><td>0.002</td><td>0.01</td><td>0.01</td><td>0.01</td><td>0.0075</td><td>0.0055<sup>e</sup></td><td>0.0075</td><td></td></tr><tr><td>pK<sub>2</sub></td><td>0.02</td><td>0.005</td><td>0.011</td><td>0.005</td><td>0.02</td><td>0.02</td><td>0.02</td><td>0.015</td><td>0.01<sup>e</sup></td><td>0.015</td><td></td></tr><tr><td>pK<sub>8</sub></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td>0.01</td><td>0.01</td><td></td></tr><tr><td>pK<sub>W</sub></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td>0.01</td><td>0.01</td><td></td></tr><tr><td>pK<sub>A</sub></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td>0.02<sup>f</sup></td><td>0.02</td><td></td></tr><tr><td>pK<sub>c</sub></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td>0.02<sup>f</sup></td><td>0.02</td><td></td></tr><tr><td>B<sub>T</sub>/S</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td>0.02<sup>g</sup></td><td></td></tr></tbody></table>

<sup>a</sup> These variables are defined as in chapter 2 of Dickson et al. (2007) except that here we write pCO₂ rather than p(CO₂). The equilibrium constants K₀, K₁, K₂, K₃ etc. are in their conventional logarithmic form: pK = −log₁₀K, while S is practical salinity. Thus whenever we refer to [H⁺] we mean the amount content (mol kg⁻¹) of “total hydrogen ion” (i.e., implicitly accounting for interactions with sulfate). Thus here as in Dickson et al. (2007), pH is defined as −log₁₀[H⁺] rather than the more conventional log₁₀(a[H⁺] + 1) (Buck et al., 2002). This definition ensures compatibility with the recommended (Dickson et al., 2007) stoichiometric acid dissociation constants for seawater media.

<sup>b</sup> Input from Dickson and Riley (1978), Millero (1995), Millero et al. (2006), and Millero (2007), respectively. M06 also estimated u(pK₁) and u(pK₂) for earlier formulations. Numbers from DR78 and M07 are based on precision, not accuracy.

<sup>c</sup> Input from Dickson (2010), Table 1.4 for (a) reference methods, (b) state-of-the-art, and (c) other.

<sup>d</sup> Input from McLaughlin et al. (2015)

<sup>e</sup> Based on the sample standard deviation of the difference between observed and predicted values for the Mehrbach et al. (1973) constants refitted by Lueker et al. (2000).

<sup>f</sup> Based on a 5% relative uncertainty (precision) from Mucci (1983)

<sup>g</sup> Relative uncertainty (equivalent to 2%)

But to what extent do uncertainties in the equilibrium constants matter? Dickson and Riley (1978) found that these uncertainties were usually a minor contributor to the overall propagated uncertainty, unlike uncertainties in the input pair of carbonate system variables. But decades have elapsed, and measurement accuracy and precision have improved. Millero (1995) used lower uncertainty estimates for both input variables and constants, finding that contributions from the constants to the overall propagated uncertainty were no longer a minor contributor although they were generally outweighed by the contribution from the input pair of measurements. Dickson (2010) reevaluated propagated uncertainties, using three sets of estimates for measurement uncertainties, while maintaining the same set of uncertainty estimates for the dissociation constants and using the same sensitivities from Dickson and Riley (1978) (Table 1). In that update, uncertainties from the constants generally dominated the overall propagated uncertainty when present-day, state-of-the-art methods were used with reference materials; conversely, with other methods, uncertainties in measurements of the input pair contributed most to the overall propagated uncertainty. Yet propagated uncertainties will differ with choices for input uncertainties, with geographic variations in input variables, and with consequent changes to the sensitivities of computed variables to input variables. The original sensitivities from Dickson and Riley (1978) remain in use today (e.g., Dickson, 2010; Sutton et al., 2016). However those were provided for only one particular set of seawater conditions (temperature T = 25C, salinity S = 35, carbonate alkalinity A<sub>C</sub> = 2248 µmol kg<sup>-1</sup>, total dissolved inorganic carbon C<sub>T</sub> = 2017 µmol kg<sup>-1</sup>, and zero nutrients). Much less of a concern is that the first and second dissociation constants of carbonic acid, K₁ and K₂, used by Dickson and Riley (1978) have changed. When those are converted from the NBS to the total scale, they agree within about 1% (at 25C) of the values now recommended for best practices (Lueker et al., 2000). The most important factor for uncertainty propagation of the marine CO₂ system is the choice of input uncertainties themselves, particularly for the dissociation constants. Choices of standard uncertainties for pK₁ and pK₂ (u(pK₁) and u(pK₂)) have ranged from 0.2 to 0.75 times the values used by Dickson and Riley (1978) (Table 1). Thus the proposed uncertainties in K₁ and K₂ (u(K₁) and u(K₂)) have also ranged from about 20 to 75% of the

Dickson and Riley (1978) values because for small uncertainties u(Kᵢ)/Kᵢ = 2.3 u(pKᵢ).

In recent years, uncertainties of carbonate system variables have been discussed in the context of establishing a Global Ocean Acidification Observing Network (GOA-ON) (Newton et al., 2015). To ensure that measurements are of appropriate quality to address the relevant problems, GOA-ON has proposed two levels of uncertainty. The GOA-ON Weather goal is aimed at assessing spatial and short-term variations (e.g., diurnal variability), while its Climate goal is stricter, focusing on deciphering decadal trends. To frame its Weather goal, GOA-ON has proposed that the relative uncertainty in calculated carbonate ion concentration [CO₃²⁻] must be < 10%, while for the Climate goal, the uncertainty for a difference in [CO₃²⁻] over time must be < 1%. While the Weather goal focuses on individual measurements, the Climate goal emphasizes the need to identify trends. These thresholds have been used to back calculate the corresponding maximum permissible uncertainties in measured input variables. For the Weather goal, GOA-ON estimates that measurement uncertainties must be no larger than 0.02 for pH, 10 µmol kg⁻¹ for total alkalinity A<sub>T</sub> and dissolved inorganic carbon C<sub>T</sub>, while pCO₂ must have a relative uncertainty of no more than 2.5%. For the Climate goal, corresponding estimates are 0.003 for pH, 2 µmol kg⁻¹ for A<sub>T</sub> and C<sub>T</sub>, and 0.5% for pCO₂. A recent interlaboratory comparison suggests that most research groups currently measuring ocean pH, A<sub>T</sub>, and C<sub>T</sub> are able to make those measurements within the criteria of the Weather goal, whereas rather few were able to achieve the Climate goal (Bockmon and Dickson, 2015). The same study hints that many groups may underestimate the true uncertainty in their measurements. Falling into this category would be research groups whose uncertainty estimates are based only on repeatability (short-term precision), i.e., not considering other possible sources of uncertainty (De Bièvre, 2008).

Traditionally, uncertainty propagation has involved identifying sources of bias (systematic errors), correcting for those as best as possible, and then propagating remaining random uncertainties (Taylor, 1996). But that changed when the International Organization for Standardization (ISO) published the Guide to the Expression of Uncertainty in Measurement (GUM, 1993), followed two years later by Eurachem's guide, subsequently
updated (Ellison and Williams, 2012), that emphasized GUM's application to analytical chemistry. GUM emphasized that the overall uncertainty includes all sources of uncertainty, both random and systematic, that both those kinds of uncertainties should be expressed as standard deviations, and that systematic as well as random uncertainties should be propagated together. GUM further clarified that if its sign and magnitude are known, the systematic error should be removed before propagating uncertainties while noting that such correction is imperfect and leaves a residual bias that should itself be treated as an uncertainty contribution; otherwise systematic uncertainty should be treated in the same way as random uncertainty, whether it was determined statistically (Type A) or by other means (Type B), such as by comparison to reference materials or by professional judgment. Although revolutionary at the time, the GUM principles have been widely adopted by the community, particularly by metrology institutes such as the National Institute of Standards and Technology (NIST).

Following suit, we adopt GUM's technical terminology and approach. For example, error and uncertainty are not synonyms. The error is the difference between a measured value and the true value. Since the true value cannot be known, neither can the error. The error may be positive or negative, while the uncertainty is always non-negative. That is, the uncertainty is the half-width of the interval around a measured value that is expected to contain the true value given a certain probability. The uncertainty of a measurement when characterized by its standard deviation is denoted as the standard uncertainty u. When the standard uncertainties of the input variables are used to propagate the uncertainty of a calculated variable, as outlined in GUM, the result is referred to as the combined standard uncertainty u_c. Often, u_c is multiplied by a coverage factor k to obtain the expanded uncertainty U, the half-width of a confidence interval that is expected to include the true value. For example, with k = 1, 2, or 3, there is 68, 95, or 99.7% confidence that the true value falls within the ± ku_c interval centered around the calculated value, assuming the variable is normally distributed.

Given the basic scientific requirement to understand uncertainty, experts from the ocean acidification community have emphasized the need to enhance existing CO₂ system software packages to allow them to also propagate uncertainties (Martz et al., 2015). Here our main aim is to do just that, thus providing marine scientists with public tools to easily propagate uncertainties using (1) automatically calculated sensitivities, which vary from region to region, and (2) input standard uncertainties that may be user specified. Standard uncertainties should be estimates that include both random and systematic components. Standard uncertainties based only on repeatability of measurements (precision) will underestimate the overall propagated uncertainty. Input uncertainties should be given as standard uncertainties u, as adopted here, so that computed outputs are the combined standard uncertainties u_c. In addition, we aim to quantify the contributions to the combined standard uncertainties calculated by the enhanced software packages and offer the ability to assess the extent to which they are affected by using alternative measurement techniques with lower or higher standard uncertainties.

## 2. Methods

Below we detail the approaches taken to compute sensitivities and propagate uncertainties, explain our choice of two sets of input uncertainties, and illustrate how results from the new routines can be synthesized graphically to provide a wider perspective.

### 2.1. Approaches

The classic way to propagate uncertainties is to use a first-order Taylor series expansion. For example, for an equation with two independent variables, y = f(x₁, x₂), the combined standard uncertainty in the dependent variable y is

$$u_c^2(y) = \left(\frac{\partial y}{\partial x_1}\right)^2 u^2(x_1) + \left(\frac{\partial y}{\partial x_2}\right)^2 u^2(x_2) + 2\left(\frac{\partial y}{\partial x_1}\right)\left(\frac{\partial y}{\partial x_2}\right) u(x_1, x_2) \quad (1)$$

where squared uncertainties are treated as a variances, i.e., u_c²(y) is the squared combined standard uncertainty for y, u²(x₁) is the squared standard uncertainty for x₁, and u²(x₂) is that for x₂, while the terms in parentheses on the right-hand side are the sensitivities of the dependent variable to each independent variable and u(x₁, x₂) is the covariance between the standard uncertainties in x₁ and x₂. The covariance is the product of the individual standard uncertainties and the correlation coefficient, namely u(x₁, x₂) = u(x₁) u(x₂) r(x₁, x₂) where -1 ≤ r(x₁, x₂) ≤ 1. Many are familiar with Eq. (1) only when the final covariance term is neglected, a simplification that is justified when x₁ and x₂ are uncorrelated; sometimes though, that is a poor assumption (Taylor, 1996; Tellinghuisen, 2001). Generalizing for n dependent variables, the propagation-of-uncertainties equation becomes

$$u_c^2(y) = \sum_{i=1}^{n} \sum_{j=1}^{n} \left( \frac{\partial y}{\partial x_i} \right) \left( \frac{\partial y}{\partial x_j} \right) u(x_i, x_j) \quad (2)$$

which follows the classical form of Eq. (1) given that the covariance of a variable with itself is its variance u(xᵢ, xᵢ) = u²(xᵢ) and that covariances are also symmetric, i.e., u(xᵢ, xⱼ) = u(xⱼ, xᵢ). Dickson and Riley (1978) applied the classic approach of propagation of uncertainties to the ocean carbonate system while neglecting covariances and expressing the equation in terms of relative uncertainties. That is, they neglected terms in Eq. (2) where i ≠ j, divided both sides of that by y², and multiplied the right side by (xᵢ/xⱼ)² to come up with

$$\left(\frac{u_c(y)}{y}\right)^2 = \sum_{i=1}^{n} \left(\frac{\partial y/y}{\partial x_i/x_i}\right)^2 \left(\frac{u(x_i)}{x_i}\right)^2, \quad (3)$$

The left-hand side is the square of the relative combined standard uncertainty in y, a function of the right-hand side's squared relative standard uncertainties of each input variable (u(xᵢ)/xᵢ), multiplied by the square of the associated relative sensitivity term (∂y/y)/(∂xᵢ/xᵢ), which is the same as (∂ ln(y))/(∂ ln(xᵢ)). Thus it is simple to understand the percent relative combined standard uncertainty in y as a function of a 1% relative standard uncertainty in each xᵢ, neglecting other contributions. The general approach shown in Eq. (2) may also be expressed in matrix form (Appendix A).

Eq. (2) is the most common form of the Method of Moments because it uses a first-order approximation of y to estimate that dependent variable's second moment, i.e., its standard deviation u_c(y); when covariance terms are neglected, it simplifies to the Gaussian approach (Taylor, 1996; Kirchner, 2001). We implemented both approaches in four public software packages that are commonly used to compute marine carbonate chemistry: CO2SYS-Excel (Pierrot et al., 2006), CO2SYS-MATLAB (van Heuven et al., 2011), mocsy (Orr and Epitalon, 2015), and seacarb (Proye and Gattuso, 2003; Lavigne and Gattuso, 2011; Gattuso et al., 2018). In the current implementation, users may specify the covariance between the uncertainties of the input pair of carbonate system variables. Other covariances are assumed to be negligible. Users may also exploit the intermediate routines that compute the absolute sensitivities (partial derivatives) needed for uncertainty propagation. These sensitivities of output variables to input variables, also known as buffer factors (Frankignoulle, 1994), are useful by themselves to assess rates of change and help deconvolve responsible mechanisms using Taylor series. A third approach was included in one public package for comparison. Rather than calculating sensitivities, the Monte Carlo approach relies on random number generation. That is, it computes the uncertainty in the output variable as the standard deviation of results from many calculations, each of which randomly selects each input variable from an assumed probability distribution, typically Gaussian, characterized by its input value and uncertainty (mean and standard deviation). Conversely, the Gaussian and Method of Moments approaches are valid independent of the probability distribution, which is often unknown.

### 2.2. Sensitivities

To propagate uncertainties with both the Gaussian approach and the Method of Moments, we must first compute the sensitivity of each computed variable to each input variable. These sensitivities were computed
numerically, optimizing the numerical step size by comparing results to
analytical solutions and so-called automatic derivatives.

A routine *derivnum* was coded to compute all partial derivatives numerically with centered differences, a commonly used approach but one for which a reference should be used to optimize the numerical step size. Optimal step sizes differ between derivatives. In parallel, we coded a reference routine to provide analytical solutions to some of the partial derivatives of the carbonate system needed for uncertainty propagation. That code was built upon an existing routine in seacarb, *buffesm* developed by Orr (2011) using equations from Egleston et al. (2010) with corrections to typographical errors consistent with those from Álvarez et al. (2014). Yet Egleston et al. (2010) ignored contributions to the total alkalinity from phosphoric and silicic acid systems. Hence both systems were added here following the procedure outlined in Appendix B. The modified version of *buffesm* used to make these calculations is available in the seacarb package.

A second reference routine *derivauto* was developed in the mocsy package using an approach known as automatic differentiation to compute all relevant partial derivatives (buffer factors). In some languages such as Fortran 90, automatic differentiation requires minimal coding, taking advantage of operator overloading, to adapt an existing routine to compute derivatives of its output variables with respect to its input variables. More importantly, its results are known to be as accurate as those from analytical solutions. Hence we implemented Dual Number Automatic Differentiation (Yu and Blair, 2013) in the *mocsy* package, modifying its routines that compute equilibrium constants and carbonate system variables. By comparing results from *derivnum* to those from *derivauto*, we optimized the step size used to compute numerical derivatives in *derivnum*. Then *derivnum*, including its optimized step sizes for each derivative, was translated into three other languages, from Fortran 90 (for *mocsy*) to R (seacarb), MATLAB (CO2SYS-MATLAB), and Visual Basic for Applications (CO2SYS-Excel).

## 2.3. Uncertainty propagation

Other new routines were written to call the numerical derivative routines and use their computed partial derivatives to propagate uncertainties with either the Method of Moments or the Gaussian approach. Both of those methods were implemented in all four public packages. For further comparison, the Monte Carlo approach was implemented in one of the packages (seacarb).

## 2.4. Standard uncertainties

To propagate uncertainties, we considered standard uncertainties for each of the variables in Table 1. We propagated these standard uncertainties through the various calculations, in some cases distinguishing contributions from systematic and random components, unlike previous assessments. Earlier estimates of standard uncertainties (Table 1) were used here for three out of the four commonly measured carbonate system variables ($A_T$, $C_T$, $pCO_2$), while for the fourth (pH), we attempted to distinguish systematic and random contributions to the standard uncertainty because the former have often been neglected. Standard uncertainties in state-of-the-art oceanographic pH measurements are typically estimated as $u(pH) = 0.003$ or less. Yet there is additional uncertainty, usually neglected, because buffer solutions are used to calibrate pH measurements, including those made by indicator dye as well as those made potentiometrically (Buck et al., 2002; Marion et al., 2011). In addition, there is another uncertainty for pH due to its temperature correction when it is not measured at in situ pH. Hence our subsequent propagation-of-uncertainty calculations use 0.003 as the random component of standard uncertainty in pH and 0.01 as an approximation for the total standard uncertainty (random plus systematic).

The random component and combined random and systematic components (total) of selected standard uncertainties (Table 1) were propagated to assess the random and systematic contributions to combined standard uncertainties. Those components were distinguished when assessing contributions to the combined standard uncertainties of derived variables using stacked bar charts. Conversely, only the total contribution of the two

components was accounted for when using a new type of diagram to assess
how the combined standard uncertainty of a given derived variable changes
across the ranges of standard uncertainties of the input pair.

Following GUM, the standard uncertainty of an equilibrium constant may be pictured as having random and systematic components, which when added in quadrature give its overall standard uncertainty. For a particular constant, we assume that the random component of its standard uncertainty may be estimated from the internal consistency of the data set, e.g., from the goodness of fit to an appropriate interpolating expression in S and T. For that component, we rely on estimates from Lueker et al. (2000) who found that in their refit of the Mehrbach et al. (1973) data, the sample standard deviation of the difference between observed and predicted values is 0.0055 for $pK_1$ and 0.010 for $pK_2$.

The systematic component of the standard uncertainty of an equilibrium constant is more difficult to assess but can be approximated by comparing results from different formulations of the same constant. Thus we compared three sets of formulations for $K_1$ and $K_2$: (1) the set from Lueker et al. (2000) that is applicable for salinities between 19 and 43 and recommended for best practices (Dickson et al., 2007); (2) a commonly used set that is applicable over a wider salinity range (0 to 40) from (Dickson and Millero (1987), and (3) a newer set based on more measurements and applicable over the same salinity range (Waters et al., 2014), an update to previous sets from the same research group (Millero et al., 2006; Millero, 2010). The $K_1$ and $K_2$ computed from these formulations are provided either directly on the total hydrogen ion scale (Dickson et al., 2007; Waters et al., 2014) or on the seawater scale (Dickson and Millero, 1987); the latter results were converted to the total scale. The Waters et al. (2014) and Lueker et al. (2000) sets of formulations agree within ± 0.007 for $pK_1$ and within ± 0.01 for $pK_2$ over the valid salinity range for the latter set; however, below that range, differences are much larger (Fig. 1). Comparison of the two formulations covering the full salinity range Waters et al., 2014, Dickson and Millero, 1987) show that they agree within ± 0.015 for $pK_1$ and ± 0.03 for $pK_2$, values that we use as 2u estimates for overall standard uncertainties (dominated by the systematic component) in subsequent uncertainty propagation. Another systematic contribution to the standard uncertainty could, in principle, be assigned either to pH or to the $pK$ values, with our preference being the latter. Namely, it is the extent to which the exact implementation of a pH scale differs between the pH measurement and the $pK$ determination. However, it may be included implicitly already insofar as we approximate the overall standard uncertainty of the equilibrium constants. Standard uncertainties in other equilibrium constants matter much less for uncertainty propagation, as shown later. Hence no attempt is made here to separate their random and systematic components. Thus we consider only their overall standard uncertainties based on previous estimates (Table 1).

Finally, our uncertainty propagation also accounts for the systematic component of standard uncertainty in the total-boron-to-salinity ratio ($B_T/S$). The two available formulations (Uppström (1974), Lee et al. (2010) differ by about 4%, leading in some cases to substantial differences in computed variables, e.g., of 4–6 µatm in computed surface $pCO_2$ when using the $A_T-C_T$ input pair (Orr and Epitalon, 2015). We assign a 2% default relative standard uncertainty to $B_T/S$ assuming that the difference between the two formulations is 2u (Table 1).

## 2.5. Error-space diagram

To illustrate how changes in standard uncertainties of a chosen input pair affect a derived variable's combined standard uncertainty, we constructed a new type of diagram (Fig. 2). This error-space diagram plots contours of the combined standard uncertainty as a function of the range of standard uncertainties in each member of the input pair after prescribing the standard uncertainties for the equilibrium constants as given in Table 1. Thus, it is analogous to the more conventional diagram where contours of a derived variable (e.g., $pCO_2$) are shown as a function of the two members of an input pair (e.g., $A_T$ and $C_T$) except that the axes and contours are the uncertainties. Once drawn for a particular set of conditions, error-space
<figure><img src="image_5.png" /><figcaption>**Fig. 1.** Contour plots showing systematic differences between formulations for (left) pK₁, (center) pK₂, and (right) pK, where K = K₁/K₂: (top row) formulations of Waters et al. (2014) minus Lueker et al. (2000) and (bottom row) formulations of Waters et al. (2014) minus Dickson and Millero (1987), the latter of which combines the measured data of Hansson (1973) and Mehrbach et al. (1973) and extrapolates from salinities of ~20 to zero. Both comparisons are perilous below S~20, which is beyond the valid salinity range of the measured data whether measured by Mehrbach et al. (and used in Lueker et al.) or by Hansson. That low-salinity region is included though as a warning about the large associated uncertainties. In the bottom center and right panels, contour lines converge towards zero difference near S = 0, e.g., for K₂ from the largest differences located around S = 5. Red dashed lines are negative, while blue solid lines are positive or zero. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

diagrams can be used to estimate the combined standard uncertainties graphically for the same conditions without further recourse to the uncertainty propagation routines. In addition, they allow users to quickly assess how using other methods with different standard uncertainties, which may also differ in cost and convenience, will affect the combined standard uncertainty of a derived variable.

In addition to showing contour lines for the combined standard uncertainty of a derived variable (Fig. 2a), an error-space diagram also includes two other features. It indicates (1) the *pair-constants curve* (Fig. 2b), along which the constants and the input pair contribute equally and (2) the *pair line* (Fig. 2c), along which each member of the input pair contributes equally. Inside the pair-constants curve, most of the combined standard uncertainty is from the constants; outside that curve, most of the combined standard uncertainty is from the input pair. Below the pair line, most of contribution to the combined standard uncertainty from the input pair comes from the member indicated on the x-axis label; above the pair line, most of that contribution comes from the member indicated on the y-axis label. The pair line and the pair-constants curve are calculated as described in Appendix C. Generally it will be seen that working to reduce the

uncertainties of members of the input pair by improving methods has a diminishing return on improving the combined standard uncertainty as one crosses the pair-constants curve and moves towards the origin. At the origin of an error-space diagram, the combined standard uncertainty is that attributable only to the constants. The basic error-space diagram provides results for one set of input conditions (e.g., A_T, C_T, T, and S), but up to two sets of input conditions may be shown if they differ substantially, e.g., average conditions for the tropics and Southern Ocean (Fig. 2d).

As implemented here, error-space diagrams neglect the minute contributions from standard uncertainties in *T* and *S*, for which oceanographic measurements are usually highly accurate (e.g., *u(T)*0.01C and *u(S)*0.01), as well as standard uncertainties from *P*<sub>T</sub> and *S*<sub>T</sub>, whose contributions to the combined standard uncertainty are negligible even in high-nutrient regions. Although *P*<sub>T</sub> and *S*<sub>T</sub> significantly affect values of computed variables in high-nutrient waters when *A*<sub>T</sub> is a member of the input pair (Orr et al., 2015), the joint contribution of their standard uncertainties to the combined standard uncertainty does not usually exceed 0.1%. Smaller still are the contributions to the combined standard uncertainty from the typical standard uncertainties in *T* and *S* mentioned above, which always result in less
<figure><img src="image_6.png" /><figcaption>**Fig. 2.** Anatomy of an error-space diagram. As its base, this diagram includes (a) contours of the percent relative combined standard uncertainty of a calculated variable (blue), in this case ΩA, as a function of the standard uncertainties in the input pair, in this case C<sub>T</sub> (x-axis) and A<sub>T</sub> (y-axis). Then in panel (b) is added the *pair-constants curve*, along which the contribution to the combined standard uncertainty from the input pair is equal to that from the constants, the latter being indicated by the blue number below the origin. In addition, in panel (c) is added the *pair line* along which the two members of the input pair contribute equally to the combined standard uncertainty. Finally, panel (d) shows results not only for average surface water of the Southern Ocean (blue), as before, but also for the tropics (red), including the contribution to the combined standard uncertainty from only the constants (red number to the left of the origin). A symbol (x) is also added to indicate an arbitrary reference for the standard uncertainties for the input pair, in this case 2 µmol kg<sup>-1</sup> for each of C<sub>T</sub> and A<sub>T</sub> (Table 1). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

than a 0.02% relative uncertainty in derived variables. Thus these diagrams focus entirely on the main contributors, i.e., standard uncertainties from the equilibrium constants and from the chosen input pair of carbonate system variables.

To facilitate future production of error-space diagrams by others, we offer an interactive interface available at https://katirg.shinyapps.io/marineco2errorpropagation for which users may raise issues or contribute improvements via GitHub at https://github.com/KatiRG/marineCO2ErrorPropagation.

## 3. Results

Packages were evaluated by comparing their numerical derivatives to analytic solutions and automatic derivatives and by comparing their combined standard uncertainties computed with a common approach to those from two other methods as well as those from a previous study. Packages were also compared to one another. Subsequently, the new routines from one package were used to characterize combined standard uncertainties as a function of standard uncertainties of input variables for two contrasting ocean regions by means of error-space diagrams.

### 3.1. Agreement

Comparison first focused on the consistency of computed partial derivatives (sensitivities) and combined standard uncertainties among packages, among approaches, and in comparison to previous work. Comparison of these sensitivities among packages reveals that they agree to at least four

significant figures, i.e., within at least 0.01% (Tables 2 and 3).

Results from the seacarb package's numerical derivatives generally agree to at least five significant figures (within 0.001%) with the analytical solution from the same package. Results from mocsy's numerical derivatives and its automatic derivatives also agree within at least four significant figures. Overall, the buffer factors calculated from the four packages and from the analytical solutions always agree to better than 0.01%. Similar agreement is found for the combined standard uncertainties with packages differing by < 0.008% for computed concentrations and partial pressures and by < 0.03% for ΩA and ΩC (Table 4). Although still small, the differences for the Ω's are larger because the formulation for the [Ca²⁺]/S ratio differs slightly among packages (Orr et al., 2015).

For a historical perspective, Table 5 shows our sensitivities and combined standard uncertainties calculated with the same input conditions as those from Dickson and Riley (1978). Our K₁ and K₂ values are only 0.4% and 1.2% smaller (when converted to the same pH scale), being calculated following Lueker et al. (2000) as recommended for best practices (Dickson et al., 2007). For the pH-A<sub>T</sub> and pH-C<sub>T</sub> input pairs, calculated sensitivities generally differ by a few percent at most from those of Dickson and Riley (1978) except that our sensitivities of computed C<sub>T</sub> and A<sub>T</sub> to [H⁺] are up to 36% higher. For the pCO₂-A<sub>T</sub> and pCO₂-C<sub>T</sub> input pairs, the magnitude of our computed sensitivities are within 12% of those from Dickson and Riley's estimates, except for our sensitivities of computed A<sub>T</sub> or C<sub>T</sub> to input pCO₂, K₀, and K₁, which are up to 36% larger. There is also a sign inconsistency with the pCO₂-A<sub>T</sub> pair, for which our computed ∂ ln [CO₃²⁻]/∂ ln (K₂) is negative while Dickson and Riley's estimate is positive due to a typo in the latter. For the pCO₂-pH input pair, all sensitivities are within 5% of those
<table><caption>**Table 2**
Comparison of partial derivatives across packages calculated from surface-ocean global means<sup>a</sup> for A<sub>T</sub>, C<sub>T</sub>, T, and S.</caption><thead><tr><th></th><th>[H<sup>+</sup>]</th><th>pCO<sub>2</sub></th><th>[CO<sub>2</sub><sup>+</sup>]</th><th>[HCO<sub>3</sub><sup>-</sup>]</th><th>[CO<sub>3</sub><sup>2-</sup>]</th><th>Ω<sub>A</sub></th></tr></thead><tbody><tr><td colspan="7"><strong>∂/∂A<sub>T</sub></strong></td></tr><tr><td>co2sys<sup>b</sup></td><td>-0.02558196<sup>b</sup></td><td>-1.192136<sup>c</sup></td><td>-0.04073334</td><td>-0.6310411</td><td>0.6717745</td><td>0.01037201</td></tr><tr><td>Mocsy</td><td>-0.02558194</td><td>-1.192138</td><td>-0.04073333</td><td>-0.6310401</td><td>0.6717725</td><td>0.01037431</td></tr><tr><td>Ref.<sup>d</sup></td><td>-0.02558194</td><td>-1.192138</td><td>-0.04073333</td><td>-0.6310401</td><td>0.6717725</td><td>0.01037431</td></tr><tr><td>Ref.<sup>e</sup></td><td>-0.02558201</td><td>-1.192140</td><td>-0.04073345</td><td>-0.6310539</td><td>0.6717729</td><td>0.01037431</td></tr><tr><td>Seacarb</td><td>-0.02557970</td><td>-1.192033</td><td>-0.04072973</td><td>-0.6310761</td><td>0.6718058</td><td>0.01037482</td></tr><tr><td>Ref.<sup>f</sup></td><td>-0.02558003</td><td>-1.192046</td><td>-0.04073025</td><td>-0.6310842</td><td>0.6718144</td><td>0.01037496</td></tr><tr><td colspan="7"><strong>∂/∂C<sub>T</sub></strong></td></tr><tr><td>co2sys<sup>b</sup></td><td>0.02810982</td><td>1.461280</td><td>0.04992956</td><td>1.584240</td><td>-0.6341699</td><td>-0.00979141</td></tr><tr><td>Mocsy</td><td>0.02810980</td><td>1.461283</td><td>0.04992955</td><td>1.584239</td><td>-0.6341677</td><td>-0.00979357</td></tr><tr><td>Ref.<sup>d</sup></td><td>0.02810980</td><td>1.461283</td><td>0.04992955</td><td>1.584239</td><td>-0.6341677</td><td>-0.00979357</td></tr><tr><td>Ref.<sup>e</sup></td><td>0.02810988</td><td>1.461285</td><td>0.04992968</td><td>1.584255</td><td>-0.6341683</td><td>-0.00979358</td></tr><tr><td>Seacarb</td><td>0.02810754</td><td>1.461163</td><td>0.04992547</td><td>1.584277</td><td>-0.6342023</td><td>-0.00979411</td></tr><tr><td>Ref.<sup>f</sup></td><td>0.02810790</td><td>1.461178</td><td>0.04992604</td><td>1.584286</td><td>-0.6342118</td><td>-0.00979425</td></tr><tr><td colspan="7"><strong>∂/∂T</strong></td></tr><tr><td>co2sys<sup>b</sup></td><td>0.2525344</td><td>12.57638</td><td>0.1332509</td><td>-0.5848728</td><td>0.4516218</td><td>0.01716385</td></tr><tr><td>Mocsy</td><td>0.2525339</td><td>12.57640</td><td>0.1332499</td><td>-0.5848648</td><td>0.4516475</td><td>0.01716809</td></tr><tr><td>Ref.<sup>c</sup></td><td>0.2525339</td><td>12.57640</td><td>0.1332499</td><td>-0.5848648</td><td>0.4516475</td><td>0.01716809</td></tr><tr><td>Seacarb</td><td>0.2525121</td><td>12.57521</td><td>0.1332345</td><td>-0.5849710</td><td>0.4517364</td><td>0.01717017</td></tr><tr><td colspan="7"><strong>∂/∂S<sup>g</sup></strong></td></tr><tr><td>co2sys<sup>b</sup></td><td>0.2099016</td><td>8.299964</td><td>0.2268601</td><td>1.004573</td><td>-1.231433</td><td>-0.03713724</td></tr><tr><td>Mocsy</td><td>0.2099018</td><td>8.299996</td><td>0.2268604</td><td>1.004582</td><td>-1.231429</td><td>-0.03714548</td></tr><tr><td>Ref.<sup>d</sup></td><td>0.2099018</td><td>8.299996</td><td>0.2268604</td><td>1.004582</td><td>-1.231429</td><td>-0.03714548</td></tr><tr><td>Seacarb</td><td>0.2098768</td><td>8.298895</td><td>0.2268276</td><td>1.004441</td><td>-1.231269</td><td>-0.03714426</td></tr></tbody></table>

<sup>a</sup> Input conditions: AT = 2300 µmol kg<sup>-1</sup>, CT = 2000 µmol kg<sup>-1</sup>, T = 18°C, S = 35, PT = 0 µmol kg<sup>-1</sup>, and S iT = 0 µmol kg<sup>-1</sup>.

<sup>b</sup> CO2SYS-Matlab

<sup>c</sup> The first number in the table is for ∂[H<sup>+</sup>]/∂AT, the second number for ∂pCO2/∂AT, both from the CO2SYS-Matlab package.

<sup>d</sup> Automatic derivative in mocsy (derivauto.f90 routine)

<sup>e</sup> Analytic derivative in mocsy (buffesm.f90 routine)

<sup>f</sup> Analytic derivative in seacarb (buffesm.R function)

<sup>g</sup> Since salinity is on the practical salinity scale, no units are given for the denominator of ∂/∂S partial derivatives

from Dickson and Riley. For the A<sub>T</sub>-C<sub>T</sub> pair, most sensitivities are lower but remain within 30% of the Dickson and Riley estimates, except for sensitivities to K<sub>2</sub> of computed [HCO<sub>3</sub><sup>-</sup>] and [CO<sub>3</sub><sup>2-</sup>], which are about 4 and 7 times larger. These differences are partly explained by our use of A<sub>T</sub> versus Dickson and Riley's use of carbonate alkalinity A<sub>C</sub> (Sect. 4.1). Despite roughly consistent sensitivities between the two studies as well as generally comparable contributions to the combined standard uncertainties from the input pair, our overall combined standard uncertainties in Table 5 are usually smaller than those from Dickson and Riley because our default input uncertainties for K<sub>1</sub> and K<sub>2</sub> are 25% smaller.

By default, all packages use the Gaussian approach to propagate uncertainties, an approach that assumes no correlation between the standard uncertainties of the two members of the input pair. For comparison, we also implemented the Monte Carlo approach in the errors function of the seacarb package. It too was coded to include no correlation between uncertainties of the input variables. Results from the Monte Carlo approach agree with those from the Gaussian approach, our reference, within the statistics permitted by the number of random samples taken (Fig. 3), an additional input parameter that is specified only for the Monte Carlo approach. Random sample sizes must be as large as 10<sup>5</sup> for the Monte Carlo approach to be able

<table><caption>**Table 3**
Comparison of computed partial derivatives with respect to P<sub>T</sub> and Si<sub>T</sub> for typical surface conditions<sup>a</sup>.</caption><thead><tr><th></th><th>[H<sup>+</sup>]</th><th>pCO<sub>2</sub></th><th>[CO<sub>2</sub><sup>+</sup>]</th><th>[HCO<sub>3</sub><sup>-</sup>]</th><th>[CO<sub>3</sub><sup>2-</sup>]</th><th>Ω<sub>A</sub></th></tr></thead><tbody><tr><td colspan="7"><strong>∂/∂P<sub>T</sub></strong></td></tr><tr><td>co2sys<sup>b</sup></td><td>0.02935462</td><td>1.368130</td><td>0.04674678</td><td>0.7009685</td><td>-0.7477153</td><td>-0.01154451</td></tr><tr><td>Mocsy</td><td>0.02935621</td><td>1.368207</td><td>0.04674931</td><td>0.7010073</td><td>-0.7477555</td><td>-0.01154773</td></tr><tr><td>Ref.<sup>c</sup></td><td>0.02935621</td><td>1.368207</td><td>0.04674931</td><td>0.7010073</td><td>-0.7477555</td><td>-0.01154773</td></tr><tr><td>Seacarb</td><td>0.02935388</td><td>1.368098</td><td>0.04674557</td><td>0.7010528</td><td>-0.7477984</td><td>-0.01154839</td></tr><tr><td colspan="7"><strong>∂/∂Si<sub>T</sub></strong></td></tr><tr><td>co2sys<sup>b</sup></td><td>0.001069476</td><td>0.04984502</td><td>0.001703124</td><td>0.02553836</td><td>-0.02724148</td><td>-0.0004206009</td></tr><tr><td>Mocsy</td><td>0.001069695</td><td>0.04985535</td><td>0.001703473</td><td>0.02554362</td><td>-0.02724705</td><td>-0.0004207814</td></tr><tr><td>Ref.<sup>c</sup></td><td>0.001069695</td><td>0.04985535</td><td>0.001703473</td><td>0.02554362</td><td>-0.02724705</td><td>-0.0004207814</td></tr><tr><td>Seacarb</td><td>0.001069679</td><td>0.04985458</td><td>0.001703446</td><td>0.02554693</td><td>-0.02725038</td><td>-0.0004208327</td></tr></tbody></table>

<sup>a</sup> Input conditions as in Table 2 except that PT = 2.0 µmol kg<sup>-1</sup>, and S iT = 60 µmol kg<sup>-1</sup>.

<sup>b</sup> CO2SYS-Matlab

<sup>c</sup> Automatic derivative in mocsy (derivauto.f90 routine)
<table><caption>**Table 4**

Combined standard uncertainties in variables computed with the A<sub>T</sub>-C<sub>T</sub> pair at typical surface conditions<sup>a</sup>.</caption><thead><tr><th rowspan="2"></th><th>[H<sup>+</sup>]</th><th>[CO<sub>2</sub><sup>+</sup>]</th><th>fCO<sub>2</sub></th><th>pCO<sub>2</sub></th><th>[HCO<sub>3</sub><sup>-</sup>]</th><th>[CO<sub>3</sub><sup>2-</sup>]</th><th rowspan="2">Ω<sub>A</sub></th><th rowspan="2">Ω<sub>C</sub></th></tr><tr><th>(nmol kg<sup>-1</sup>)</th><th>(μmol kg<sup>-1</sup>)</th><th>(μatm)</th><th>(μatm)</th><th>(μmol kg<sup>-1</sup>)</th><th>(μmol kg<sup>-1</sup>)</th></tr></thead><tbody><tr><td colspan="9">Without standard uncertainties from the constants</td></tr><tr><td>CO2SYS-Excel</td><td>0.07668398</td><td>0.12995259</td><td>3.7900747</td><td>3.8032995</td><td>3.4135224</td><td>1.8539528</td><td>0.02862451</td><td>0.04426536</td></tr><tr><td>CO2SYS-Matlab</td><td>0.07668398</td><td>0.12995259</td><td>3.7900747</td><td>3.8032994</td><td>3.4135224</td><td>1.8539528</td><td>0.02862451</td><td>0.04426536</td></tr><tr><td>Mocsy</td><td>0.07669003</td><td>0.12996267</td><td>3.7903686</td><td>3.8036035</td><td>3.4135222</td><td>1.8539595</td><td>0.02863104</td><td>0.04427546</td></tr><tr><td>Seacarb</td><td>0.07668414</td><td>0.12995262</td><td>3.7900755</td><td>3.8033093</td><td>3.4136038</td><td>1.8540405</td><td>0.02863230</td><td>0.04427739</td></tr><tr><td colspan="9">With default standard uncertainties from the constants</td></tr><tr><td>CO2SYS-Excel</td><td>0.19429498</td><td>0.32911625</td><td>9.6994199</td><td>9.7332644</td><td>4.4247168</td><td>3.2612783</td><td>0.15578845</td><td>0.24091348</td></tr><tr><td>CO2SYS-Matlab</td><td>0.19429498</td><td>0.32911625</td><td>9.6994199</td><td>9.7332642</td><td>4.4247168</td><td>3.2612783</td><td>0.15578845</td><td>0.24091348</td></tr><tr><td>Mocsy</td><td>0.19430788</td><td>0.32913860</td><td>9.7000746</td><td>9.7339444</td><td>4.4247421</td><td>3.2613327</td><td>0.15581821</td><td>0.24095949</td></tr><tr><td>Seacarb</td><td>0.19429739</td><td>0.32911714</td><td>9.6994403</td><td>9.7333078</td><td>4.4247272</td><td>3.2612654</td><td>0.15582641</td><td>0.24097219</td></tr></tbody></table>

to routinely obtain combined standard uncertainties that agree with those from the default Gaussian approach within < 1%. With random samples of that size, the computation time for the Monte Carlo approach is 10 times longer than the Gaussian approach for the pCO₂-pH pair but 440 times longer for the A<sub>T</sub>-C<sub>T</sub> pair based on our implementation in seacarb.

The third approach, Method of Moments, is unlike the others in that its calculated combined standard uncertainties can be affected by a user-specified correlation (r) between the standard uncertainties of the two members of the input pair. With r = 0, the Method of Moments matches results from the Gaussian approach; with nonzero values of r (between -1 and +1), the combined standard uncertainties are affected variously, depending on the computed variable and the input pair as well as the value of r. For instance, for the pCO₂-pH pair, by assuming a perfect correlation (r = 1) between the uncertainties in pCO₂ and [H⁺], the overall propagated uncertainty in computed [CO₃²⁻] is 24% less than with no correlation (Table 6). There is a reduction because of the negative sign of the covariance term (last term in Eq. (1)), i.e., twice the product of the two sensitivities, the individual uncertainties, and the correlation coefficient (between the two uncertainties):

$$2 \left( \frac{\partial[\text{CO}_3^{2-}]}{\partial p\text{CO}_2} \right) \left( \frac{\partial[\text{CO}_3^{2-}]}{\partial [\text{H}^+]} \right) u(p\text{CO}_2) u([\text{H}^+]) r(p\text{CO}_2, [\text{H}^+]) \quad (4)$$

The covariance term is negative because the two standard uncertainties are positive (by definition), the two sensitivities are opposite in sign (e.g., Table 5c), and r(pCO₂, [H⁺]) is positive. Thus combined standard uncertainties can be affected by correlations between the standard uncertainties of the two members of the input pair. However, the correlation between the standard uncertainties is not the same as a correlation between the measurements themselves. For uncertainty propagation, the interest is in the correlation between standard uncertainties of input variables, not between the values of the input variables themselves. If there is a measurement uncertainty that is common to both members of the input pair, then that commonality should be accounted for explicitly. Correlation between standard uncertainties of other input variables such as K₁ and K₂ remains possible but requires further investigation. Hence adding room for those additional correlations in the uncertainty propagation routines is left to future work.

## 3.2. Error space

To quantify combined standard uncertainties more generally and assess the potential for improvement, we used our errors function from the seacarb package to provide data for error-space diagrams (Sect. 2.5), which show how combined standard uncertainties in derived variables are affected by the range of possible uncertainties in input variables. Uncertainties were propagated for five calculated variables [CO₃²⁻], Ω<sub>A</sub>, pCO₂, [H⁺], and [HCO₃⁻], along with the ratio of the latter two [HCO₃⁻]/[H⁺], all as a function of a range of standard uncertainties in each member of six input

pairs (pH-A<sub>T</sub>, pH-C<sub>T</sub>, pCO₂-A<sub>T</sub>, pCO₂-C<sub>T</sub>, pCO₂-pH, and A<sub>T</sub>-C<sub>T</sub>). The computed combined standard uncertainties were compared using error-space diagrams drawn for each input pair and each computed variable (Figs. 4–11).

### 3.2.1. [CO₃²⁻]

Error-space diagrams for [CO₃²⁻] emphasize that its combined standard uncertainty is lowest when calculated with the A<sub>T</sub>-C<sub>T</sub> input pair and highest when calculated with the pCO₂-pH pair (Fig. 4). Each member of the A<sub>T</sub>-C<sub>T</sub> pair contributes nearly equally to the combined standard uncertainty of [CO₃²⁻] when their standard uncertainties are identical because the sensitivities of [CO₃²⁻] to A<sub>T</sub> and C<sub>T</sub> have similar magnitudes. With identical sensitivities for each member of the input pair, an error-space diagram would be perfectly symmetrical (Eq. (1)), an ideal case that is nearly reached for [CO₃²⁻] as well as most other variables computed from the A<sub>T</sub>-C<sub>T</sub> pair. For that pair, the pair-constants curves occur at about 2% relative uncertainty in [CO₃²⁻] for surface waters of both the tropics and the Southern Ocean. For other pairs, relative combined standard uncertainties in [CO₃²⁻] are larger, particularly for the pCO₂-pH input pair, and contributions of the standard uncertainties from the two members of the input pair are less symmetric. Given state-of-the-art standard uncertainties (Table 1), if pH is a member of the input pair, its contribution to the combined standard uncertainty is substantially larger than that from the other member, as revealed by the position of the reference points (crosses) which sit well below the pair line. Likewise, when pCO₂ is a member of the input pair, its contribution is larger except when the other member is pH. In the error-space diagrams, the much heavier weight of pH or pCO₂ (x axis) when paired with A<sub>T</sub> or C<sub>T</sub> (y axis) is also manifested by near-vertical lines for overall propagated uncertainty. For all pairs, current state-of-the-art standard uncertainties in each of their members place the corresponding propagated uncertainty within the pair-constants curve, where the uncertainty from the constants outweighs the uncertainty from the input pair. Thus further improvements to methods to reduce their measurement uncertainties would do little to reduce the combined standard uncertainty, given the substantial propagated uncertainty from the constants. For example, a threefold reduction of the standard uncertainty in pH from 0.01 to 0.003 lowers the combined standard uncertainty of [CO₃²⁻] by only a factor of one-seventh with the pH-A<sub>T</sub> and pH-C<sub>T</sub> pairs and one-third with the pCO₂-pH pair.

### 3.2.2. Ω<sub>A</sub>

For Ω<sub>A</sub>, error-space diagrams are similar to those for [CO₃²⁻], except that relative combined standard uncertainties are typically a few percent larger and pair-constants curves are farther from the origin, simply due to the standard uncertainty in K<sub>A</sub> (u(pK<sub>A</sub>) = 0.02), a constant (solubility product) that only affects computed Ω<sub>A</sub> (Fig. 5). The pair-constants curve for the A<sub>T</sub>-C<sub>T</sub> pair occurs at about 7% relative combined standard uncertainty for both regions, i.e., at a value that is more than twice that seen for

<footnote><sup>a</sup> Input conditions (± u): AT = 2300 ± 2 µmol kg<sup>-1</sup>, CT = 2000 ± 2 µmol kg<sup>-1</sup>, T = 18 ± 0°C, S = 35 ± 0, PT = 2.0 ± 0.1 µmol kg<sup>-1</sup>, and S iT = 60 ± 4 µmol kg<sup>-1</sup>.</footnote>
<table><caption>**Table 5**
Relative sensitivities of output to input (dy/y)/(dx/x) and propagated percent relative combined standard uncertainties $u_c$ for six input pairs<sup>a</sup>.</caption><thead><tr><th colspan="6">(a) Input (pH-A<sub>T</sub> pair) sensitivities</th><th colspan="2">u<sub>c</sub> (%)<sup>c</sup></th></tr><tr><th>Output</th><th>[H<sup>+</sup>]</th><th>C<sub>T</sub></th><th>K<sub>0</sub></th><th>K<sub>1</sub></th><th>K<sub>2</sub></th><th>(i)</th><th>(ii)</th></tr></thead><tbody><tr><td>[CO<sub>2</sub><sup>+</sup>]</td><td><sup>b</sup>1.249</td><td><sup>b</sup>1.048</td><td></td><td>-1.000</td><td>-0.212</td><td>2.9</td><td>3.5</td></tr><tr><td>[HCO<sub>3</sub><sup>-</sup>]</td><td>0.249</td><td>1.048</td><td></td><td></td><td>-0.212</td><td>0.7</td><td>1.0</td></tr><tr><td>[CO<sub>3</sub><sup>2-</sup>]</td><td>-0.751</td><td>1.048</td><td></td><td></td><td>0.788</td><td>1.8</td><td>3.3</td></tr><tr><td>C<sub>T</sub></td><td>0.136</td><td>1.048</td><td></td><td>-0.005</td><td>-0.094</td><td>0.5</td><td>0.6</td></tr><tr><td>pCO<sub>2</sub></td><td>1.249</td><td>1.048</td><td>-1.000</td><td>-1.000</td><td>-0.212</td><td>2.9</td><td>3.5</td></tr><tr><td>Ω<sub>A</sub></td><td>-0.751</td><td>1.048</td><td></td><td></td><td>0.788</td><td>1.8</td><td>5.6</td></tr></tbody><thead><tr><th colspan="6">(b) Input (pH-C<sub>T</sub> pair) sensitivities</th><th colspan="2">σ (%)</th></tr><tr><th>Output</th><th>[H<sup>+</sup>]</th><th>C<sub>T</sub></th><th>K<sub>0</sub></th><th>K<sub>1</sub></th><th>K<sub>2</sub></th><th>(i)</th><th>(ii)</th></tr></thead><tbody><tr><td>[CO<sub>2</sub><sup>+</sup>]</td><td>1.113</td><td>1.000</td><td></td><td>-0.995</td><td>-0.118</td><td>2.6</td><td>3.2</td></tr><tr><td>[HCO<sub>3</sub><sup>-</sup>]</td><td>0.113</td><td>1.000</td><td></td><td>0.005</td><td>-0.118</td><td>0.6</td><td>0.7</td></tr><tr><td>[CO<sub>3</sub><sup>2-</sup>]</td><td>-0.887</td><td>1.000</td><td></td><td>0.005</td><td>0.882</td><td>2.1</td><td>3.7</td></tr><tr><td>A<sub>T</sub></td><td>-0.130</td><td>0.955</td><td></td><td>0.005</td><td>0.090</td><td>0.6</td><td>0.6</td></tr><tr><td>pCO<sub>2</sub></td><td>1.113</td><td>1.000</td><td>-1.000</td><td>-0.995</td><td>-0.118</td><td>2.6</td><td>3.2</td></tr><tr><td>Ω<sub>A</sub></td><td>-0.887</td><td>1.000</td><td></td><td>0.005</td><td>0.882</td><td>2.1</td><td>5.9</td></tr></tbody><thead><tr><th colspan="6">(c) Input (pCO<sub>2</sub>-pH pair) sensitivities</th><th colspan="2">σ (%)</th></tr><tr><th>Output</th><th>[H<sup>+</sup>]</th><th>pCO<sub>2</sub></th><th>K<sub>0</sub></th><th>K<sub>1</sub></th><th>K<sub>2</sub></th><th>(i)</th><th>(ii)</th></tr></thead><tbody><tr><td>[CO<sub>2</sub><sup>+</sup>]</td><td></td><td>1.000</td><td>1.000</td><td></td><td></td><td>2.0</td><td>2.1</td></tr><tr><td>[HCO<sub>3</sub><sup>-</sup>]</td><td>-1.000</td><td>1.000</td><td>1.000</td><td>1.000</td><td></td><td>3.0</td><td>3.5</td></tr><tr><td>[CO<sub>3</sub><sup>2-</sup>]</td><td>-2.000</td><td>1.000</td><td>1.000</td><td>1.000</td><td>1.000</td><td>5.0</td><td>6.4</td></tr><tr><td>A<sub>T</sub></td><td>-1.192</td><td>0.954</td><td>0.954</td><td>0.954</td><td>0.203</td><td>3.3</td><td>3.8</td></tr><tr><td>C<sub>T</sub></td><td>-1.113</td><td>1.000</td><td>1.000</td><td>0.995</td><td>0.118</td><td>3.3</td><td>3.7</td></tr><tr><td>Ω<sub>A</sub></td><td>-2.000</td><td>1.000</td><td>1.000</td><td>1.000</td><td>1.000</td><td>5.0</td><td>7.8</td></tr></tbody><thead><tr><th colspan="6">(d) Input (A<sub>T</sub>-C<sub>T</sub> pair) sensitivities</th><th colspan="2">u<sub>c</sub> (%)</th></tr><tr><th>Output</th><th>A<sub>T</sub></th><th>C<sub>T</sub></th><th>K<sub>0</sub></th><th>K<sub>1</sub></th><th>K<sub>2</sub></th><th>(i)</th><th>(ii)</th></tr></thead><tbody><tr><td>[CO<sub>2</sub><sup>+</sup>]</td><td>-8.536</td><td>9.146</td><td>0.000</td><td>-0.955</td><td>0.652</td><td>2.0</td><td>3.5</td></tr><tr><td>[HCO<sub>3</sub><sup>-</sup>]</td><td>-0.875</td><td>1.835</td><td>0.000</td><td>0.009</td><td>-0.040</td><td>0.4</td><td>0.4</td></tr><tr><td>[CO<sub>3</sub><sup>2-</sup>]</td><td>6.787</td><td>-5.476</td><td>0.000</td><td>-0.027</td><td>0.268</td><td>1.3</td><td>1.7</td></tr><tr><td>[H<sup>+</sup>]</td><td>-7.661</td><td>7.311</td><td>0.000</td><td>0.036</td><td>0.692</td><td>1.7</td><td>3.0</td></tr><tr><td>pCO<sub>2</sub></td><td>-8.536</td><td>9.146</td><td>-1.000</td><td>-0.955</td><td>0.652</td><td>2.0</td><td>3.5</td></tr><tr><td>Ω<sub>A</sub></td><td>6.787</td><td>-5.476</td><td>0.000</td><td>-0.027</td><td>0.268</td><td>1.3</td><td>4.9</td></tr></tbody><thead><tr><th colspan="6">(e) Input (pCO<sub>2</sub>-A<sub>T</sub> pair) sensitivities</th><th colspan="2">σ (%)</th></tr><tr><th>Output</th><th>A<sub>T</sub></th><th>pCO<sub>2</sub></th><th>K<sub>0</sub></th><th>K<sub>1</sub></th><th>K<sub>2</sub></th><th>(i)</th><th>(ii)</th></tr></thead><tbody><tr><td>[CO<sub>2</sub><sup>+</sup>]</td><td></td><td>1.000</td><td>1.000</td><td></td><td></td><td>2.0</td><td>2.1</td></tr><tr><td>[HCO<sub>3</sub><sup>-</sup>]</td><td>0.838</td><td>0.200</td><td>0.200</td><td>0.200</td><td>-0.171</td><td>0.5</td><td>0.9</td></tr><tr><td>[CO<sub>3</sub><sup>2-</sup>]</td><td>1.676</td><td>-0.599</td><td>-0.599</td><td>-0.599</td><td>0.659</td><td>1.4</td><td>2.9</td></tr><tr><td>[H<sup>+</sup>]</td><td>-0.838</td><td>0.800</td><td>0.800</td><td>0.800</td><td>0.171</td><td>1.6</td><td>2.3</td></tr><tr><td>C<sub>T</sub></td><td>0.933</td><td>0.109</td><td>0.109</td><td>0.104</td><td>-0.071</td><td>0.4</td><td>0.5</td></tr><tr><td>Ω<sub>A</sub></td><td>1.676</td><td>-0.599</td><td>-0.599</td><td>-0.599</td><td>0.659</td><td>1.4</td><td>5.4</td></tr></tbody><thead><tr><th colspan="6">(f) Input (pCO<sub>2</sub>-C<sub>T</sub> pair) sensitivities</th><th colspan="2">σ (%)</th></tr><tr><th>Output</th><th>C<sub>T</sub></th><th>pCO<sub>2</sub></th><th>K<sub>0</sub></th><th>K<sub>1</sub></th><th>K<sub>2</sub></th><th>(i)</th><th>(ii)</th></tr></thead><tbody><tr><td>[CO<sub>2</sub><sup>+</sup>]</td><td></td><td>1.000</td><td>1.000</td><td></td><td></td><td>2.0</td><td>2.1</td></tr><tr><td>[HCO<sub>3</sub><sup>-</sup>]</td><td>0.898</td><td>0.102</td><td>0.102</td><td>0.107</td><td>-0.107</td><td>0.5</td><td>0.6</td></tr><tr><td>[CO<sub>3</sub><sup>2-</sup>]</td><td>1.795</td><td>-0.795</td><td>-0.795</td><td>-0.786</td><td>0.786</td><td>1.8</td><td>3.6</td></tr><tr><td>[H<sup>+</sup>]</td><td>-0.898</td><td>0.898</td><td>0.898</td><td>0.893</td><td>0.107</td><td>1.9</td><td>2.5</td></tr><tr><td>A<sub>T</sub></td><td>1.071</td><td>-0.117</td><td>-0.117</td><td>-0.112</td><td>0.076</td><td>0.6</td><td>0.7</td></tr><tr><td>Ω<sub>A</sub></td><td>1.795</td><td>-0.795</td><td>-0.795</td><td>-0.786</td><td>0.786</td><td>1.8</td><td>5.8</td></tr></tbody></table>

<footnote><sup>a</sup> For comparison to Table II of Dickson and Riley (1978) calculated with the same input: A<sub>T</sub> = 2355 µmol kg<sup>-1</sup> computed from A<sub>C</sub> = 2248 µmol kg<sup>-1</sup>, C<sub>T</sub> = 2017 µmol kg<sup>-1</sup>, pH = 8.10 on the total scale computed from pH = 8.25 on the NBS scale using Eq. (15) from Lueker et al. (2000), pCO<sub>2</sub> = 350 µatm, T = 25°C, and S = 35, having zero standard uncertainty in T and S, a standard uncertainty of 0.01 in pH, a relative standard uncertainty of 2% for pCO<sub>2</sub>, and relative standard uncertainties of 0.4% for A<sub>T</sub> and 0.5% for C<sub>T</sub> except for the A<sub>T</sub>-C<sub>T</sub> pair where the latter two standard uncertainties were instead 0.1 and 0.2%. Combined standard uncertainty estimates from Dickson and Riley (1978) were based on precision at that time, not accuracy.</footnote>

<footnote><sup>b</sup> The first number in the table is for $\partial \ln[CO_2]/\partial \ln[H^+]$, the second for $\partial \ln[CO_2]/\partial \ln A_T$.</footnote>

<footnote><sup>c</sup> The final two columns list combined standard uncertainties computed considering standard uncertainties in (i) only the input pair and (ii) the input pair and the equilibrium constants (total uncertainty).</footnote>
<figure><img src="image_7.png" /><figcaption>**Fig. 3.** Percent relative combined standard uncertainty of carbonate ion concentration computed with the Gaussian approach (blue dotted line) vs. the Monte Carlo approach (black dots with error bars). The Monte Carlo approach requires an input parameter to specify the number of iterations (runs). That approach was used with the same input, altering only the number of iterations to make five separate estimates of the propagated uncertainty (runs = 10⁷, where n = 2, 3, 4, 5, 6). That exercise was then repeated 9 more times for each case to compute statistics, i.e., mean ± 1 standard deviation (shown as the filled symbols and error bars). In both approaches, [CO₃²⁻] is computed from the following conditions (± u): pCO₂ = 400 ± 2 µatm, pH = 8.100 ± 0.005 (total scale), S = 35 ± 0, T = 18 ± 0 C, P_T = 2.0 ± 0.1 µmol kg⁻¹, and SiT = 60 ± 4 µmol kg⁻¹. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

[CO₃²⁻]. Because pair-constants curves are more distant from the origin owing to the relatively larger contribution to the combined standard uncertainty from the constants (due to Kₐ), improving standard uncertainties of the input pair would have even less effect on the combined standard uncertainty. Yet that also means that reducing methodological uncertainty requirements for the input pair from the state-of-the-art constraints, which are well under the pair-constants curve, results in little worsening of overall propagated uncertainties. For instance, a fivefold increase from the currently smallest possible uncertainties for the input pair (crosses in Fig. 5) adds < 1% to the overall propagated relative uncertainty in Ωₐ for input pairs where the pair-constants curve is either near the edge or outside of the plot domain (pH-A_T, pH-C_T, pCO₂-A_T, and pCO₂-C_T pairs).

### 3.2.3. pCO₂ and [CO₂*]

Three of the six input pairs include pCO₂ as a member, while for the three others it is calculated. For the latter three, the two that include pH as a member of the input pair have combined standard uncertainties in pCO₂ that vary only as a function of u(pH) throughout most of the plot domain (Fig. 6). Those close ties are explained by the near-linear relationship between [CO₂*] and [H⁺] (Orr, 2011; Kwiatkowski and Orr, 2018) and the relative insensitivity to likely uncertainties in the opposing member of each input pair (A_T or C_T). In those two plot domains, curvature is obvious only at lower u(pH) and higher u(A_T) or u(C_T) (above 10 µmol kg⁻¹). For the remaining A_T-C_T pair, relative combined standard uncertainties of pCO₂ are nearly symmetrical with respect to each member as they are for [CO₃²⁻], but with pCO₂ they are larger and the pair-constants curves are about twice as far from the origin. Thus, changes in standard uncertainties of the input pair from the same state-of-the-art reference have even less effect on the combined standard uncertainty of pCO₂ relative to that for [CO₃²⁻]. That difference is driven by the relative contribution to the combined standard uncertainty from the constants being larger for pCO₂ than for [CO₃²⁻]. For pCO₂, increasing standard uncertainties of A_T and C_T from 2 to 10 µmol kg⁻¹ results in an increase in the percent relative combined standard uncertainty from about 3 to 6% in the tropics. For the same three pairs, relative combined standard uncertainties in [CO₂*] (Fig. 7) are visually indistinguishable from those for pCO₂ (Fig. 6). When pCO₂ is a member of the

<caption>**Table 6**</caption>

Relative change in combined standard uncertainty (%) that would result from correlation between the uncertainties of pCO₂ and pH.

<table><thead><tr><th>r</th><th>[HCO<sub>3</sub><sup>-</sup>]</th><th>[CO<sub>3</sub><sup>2-</sup>]</th><th>C<sub>T</sub></th><th>A<sub>T</sub></th><th>Ω<sub>A</sub></th></tr></thead><tbody><tr><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td></tr><tr><td>0.5</td><td>-20</td><td>-11</td><td>-19</td><td>-18</td><td>-11</td></tr><tr><td>01.0</td><td>-48</td><td>-24</td><td>-44</td><td>-40</td><td>-24</td></tr></tbody></table>

input pair, it may be used to directly compute [CO₂*], knowing only K₀ and that the relationship is linear. For example at u(pCO₂) = 10 µatm, the percent relative standard uncertainty in [CO₂*] is about 3% for surface waters of both the tropics and the Southern Ocean based on their average surface composition given in Fig. 4.

### 3.2.4. pH and [H⁺]

Error-space diagrams showing absolute combined standard uncertainties in pH (Fig. 8) actually represent a relative change in [H⁺]. They may be converted to percent relative combined standard uncertainties in [H⁺] simply by multiplying by c = 2.3 × 10⁷. The resulting error-space diagrams for [H⁺] (Fig. 9) are analogous to those for pCO₂ (Fig. 6), although input pairs with pH are excluded in the former and pairs with pCO₂ are excluded in the latter, as described above. Owing to the near-linear relationship between these two calculated variables, relative combined standard uncertainties and pair-constants curves for [H⁺] computed from the pCO₂-A_T and pCO₂-C_T input pairs closely resemble those for pCO₂ computed from the pH-A_T and pH-C_T input pairs, respectively. The corresponding error-space diagrams for [H⁺] and pCO₂ each computed from the A_T-C_T pair also closely resemble one another.

### 3.2.5. [HCO₃⁻]

Unlike previous error-space diagrams, those for [HCO₃⁻] differ in several ways (Fig. 10). First, the percent relative combined standard uncertainties are around an order of magnitude lower than those for other variables, e.g., < 1% over most of the plot domain at least for the Southern Ocean, except for the pCO₂-pH pair. Second, there is much less symmetry
<figure><img src="image_8.png" /><figcaption>**Fig. 4.** Error-space diagrams of the percent relative combined standard uncertainty in [CO₃²⁻] computed from six input pairs as a function of standard uncertainties in each pair's two members. Input pair members are indicated by the x- and y-axis labels. Combined standard uncertainties depend on the composition of seawater. Blue solid lines indicate average conditions for surface waters of the Southern Ocean (A_T = 2295 µmol kg⁻¹, C_T = 2155 µmol kg⁻¹, temperature T = -0.49 C, and salinity S = 33.96), while red dashed lines indicate average conditions for tropical surface waters (A_T = 2300 µmol kg⁻¹, C_T = 1960 µmol kg⁻¹, temperature T = 27.01 C, and salinity S = 34.92). Crosses indicate the state-of-the-art standard uncertainty for each member of the input pair (next to last column in Table 1) as well as the random component of the standard uncertainty for pH (previous column). Black lines indicate the pair line (thin) and pair-constants curve (thick) for both regions (Southern Ocean, solid; tropics, dashed); at their intersections lie the pair-constants midpoints (Southern Ocean, filled circle; tropics, open circle). Pair-constants midpoints fall outside of the plot domain for four of the pairs: pH-A_T, pH-C_T, pCO₂-A_T, and pCO₂-C_T. Colored numbers in the margin near each panel's origin indicate the contribution to the overall propagated uncertainty from the constants alone for the Southern Ocean (blue) and the tropics (red). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

than seen for previous error-space diagrams between members of the A<sub>T</sub>-C<sub>T</sub> pair, while there is greater symmetry for the pCO<sub>2</sub>-pH pair. The former is explained by the high sensitivity of computed [HCO<sub>3</sub><sup>-</sup>] to C<sub>T</sub> as also evident in the error-space diagrams for pH-C<sub>T</sub> and pCO<sub>2</sub>-C<sub>T</sub>. Third, there are larger differences in patterns of percent relative combined standard uncertainties between the two regions for the four pairs where either pH or pCO<sub>2</sub> is a member.

### 3.2.6. [HCO₃⁻]/[H⁺]

The [HCO₃⁻]/[H⁺] ratio (Fig. 11) is a derived variable that has been proposed as being potentially physiologically more relevant to determining the gross calcification rate than either Ω<sub>A</sub> or Ω<sub>C</sub> because [HCO₃⁻] stimulates calcification while [H⁺] inhibits it (Bach, 2015). This derived quantity is also known as the substrate-to-inhibitor ratio (SIR) (Fassbender et al., 2016b). With the pCO₂-pH pair, overall propagated relative uncertainties of SIR are much like those for [HCO₃⁻]; for other pairs, they resemble those for [H⁺]. The reason is that the percent relative combined standard uncertainties in [H⁺] are larger than those for [HCO₃⁻] with all pairs except for pCO₂-pH. Given state-of-the-art total standard uncertainties for the four measured variables (Table 1), A<sub>T</sub>-C<sub>T</sub> is the worst pair to use to compute SIR even though the pCO₂-pH pair is the poorest combination to compute

[HCO₃⁻]. The best combinations for SIR are the four remaining pairs: pH-A<sub>T</sub>, pH-C<sub>T</sub>, pCO₂-A<sub>T</sub>, and pCO₂-C<sub>T</sub>. In the first two of those pairs, where pH is measured not calculated, it might be expected that combined standard uncertainties for SIR would be lower because u(HCO₃⁻)/[HCO₃⁻] ≪ u([H⁺])/[H⁺]. That would be the result if one considered only our random component for the standard uncertainty in pH of 0.003 (Table 1). However, using the total standard uncertainty of 0.01 raises the combined standard uncertainty in SIR to the same level as for the two other pairs where pCO₂ is a member.

Although we provide error-space diagrams for SIR, doing so appears to offer little new insight. The patterns of SIR's combined standard uncertainties follow those for [CO₃²⁻] (Fig. 4) given that SIR = [CO₃²⁻]/K₂. Likewise they follow those for Ω<sub>A</sub> (Fig. 5), when [Ca²⁺] is not manipulated, because Ω<sub>A</sub> = [Ca²⁺][CO₃²⁻]/Kₐ. Given these proportionalities, SIR does not provide independent information.

## 4. Discussion

To place results in context, let us now look further into sensitivities, approach suitability, sources of uncertainties, impacts of covariance, and uncertainty criteria for large-scale measurement programs.
<figure><img src="image_9.png" /><figcaption>**Fig. 5.** Error-space diagrams of the percent relative combined standard uncertainty in ΩA computed from six input pairs as a function of standard uncertainties in each pair's two members. Pair members and line and point signatures are as indicated in Fig. 4. Pair-constants curves and midpoints are near the edge or outside the plot domain for the pCO₂-C_T, pCO₂-A_T, and pH-A_T input pairs.</figcaption></figure>

## 4.1. Sensitivities and approaches

Our choice to use A_T rather than A_C affects the calculated sensitivities. Dickson and Riley (1978) used A_C, the preference at that time (Park, 1969; Skirrow, 1975), with which the carbonate system equations can be solved analytically; conversely, with today's choice of A_T, one must resort to

numerical techniques. Another advantage of the A_C-A_T and pCO₂-A_C pairs is that when solving the required quadratic expressions, the equilibrium constants for carbonic acid do not appear separately but only as the K₁/K₂ ratio. Hence, the sensitivity to K₁ is generally equal and opposite to that for K₂. The same symmetry holds for the pCO₂-C_T input pair as shown in Table 5. This symmetry disappears when calculating [H⁺], which involves

<figure><img src="image_10.png" /><figcaption>**Fig. 6.** Error-space diagrams of the percent relative combined standard uncertainty in pCO₂ computed from three input pairs as a function of standard uncertainties in each pair's two members. Input pair members and line and point signatures are as indicated in Fig. 4. Plots are shown only for the three pairs where the uncertainty in pCO₂ is not an input variable.</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>Fig. 7. Error-space diagrams of the percent relative combined standard uncertainty in [CO₂⁻] computed from three input pairs as a function of uncertainties in each pair's two members. Pair members and line and point signatures are as indicated in Fig. 4. Plots are shown only for the three pairs where computed [CO₂⁻] depends on both members of the input pair.</figcaption></figure>

other dependencies, as does the [H⁺]-dependent calculation of A_T from A_C. Thus present-day studies, which refer to A_T, should avoid using sensitivities to K₁ and K₂ calculated in regards to A_C.

One could question though our large sensitivities of [HCO₃⁻] and [CO₃²⁻] to K₂ calculated with the A_T-C_T input pair rather than with the A_C-C_T input pair used by Dickson and Riley (1978). To determine if they are accurate, we compared our sensitivities computed numerically for use with the Gaussian approach to those back-calculated from the Monte Carlo approach. For example, the magnitude of the absolute derivative ∂[CO₃²⁻]/∂K₂ can be back calculated by propagating uncertainties for only the input variable of interest (K₂), i.e., setting all other input uncertainties to zero. Thus only one term is left on the right-hand side of Eq. (1), which can be rearranged to solve for the derivative. For [CO₃²⁻] then,

$$ \left( \frac{\partial[\text{CO}_3^{2-}]}{\partial K_2} \right)^2 = \left( \frac{u_c (|\text{CO}_3^{2-}|)}{u(K_2)} \right)^2 \quad (5) $$

i.e., we only need to divide its component of the combined standard uncertainty from K₂ by the standard uncertainty in K₂. That was done for the Monte Carlo approach as well as the Gaussian approach. Both yield the same result for the absolute value of ∂[CO₃²⁻]/∂K₂, and both agree with the estimate from our routine that calculates numerical derivatives (derivatum). When these absolute derivatives are converted to relative derivatives, they agree with our relative sensitivity estimate ∂ ln [CO₃²⁻]/∂ ln (K₂) in Table 5. The same holds for ∂ ln [HCO₃⁻]/∂ ln (K₂). This agreement is expected

when using results from the Gaussian approach, if the back calculation is done correctly. Further agreement with results from the Monte Carlo approach, a method that does not use derivatives to propagate uncertainties, implies that our a priori numerical estimates of these derivatives are accurate. Therefore our estimates of ∂ ln [CO₃²⁻]/∂ ln (K₂) and ∂ ln [HCO₃⁻]/∂ ln (K₂) for the A_T-C_T input pair (Table 5d) appear to be correct.

The consistent sensitivities for the Gaussian and Monte Carlo approaches detailed just above and the overall agreement between their propagated uncertainties (Fig. 3) both indicate that the linear approximations inherent in the Gaussian approach do not significantly skew propagated uncertainties of the marine carbonate system. As in any Taylor series, those nonlinearities due to higher order terms may be safely neglected when the uncertainty specified for each input variable, expressed as standard deviation u(xᵢ), is sufficiently small when compared to the corresponding input variable xᵢ. Therefore our findings do not support the assertion by Lauvset and Gruber (2014) that nonlinearities in the carbonate system would cause the Gaussian approach to underestimate propagated uncertainties. Rather than the approach, it is the magnitude of the uncertainties that matters.

## 4.2. Contributions to propagated uncertainties

To synthesize relative contributions, we follow the lead of GUM (1993), rewriting Eq. (1) as

<figure><img src="image_12.png" /><figcaption>Fig. 8. Error-space diagrams of the combined standard uncertainty in pH computed from three input pairs as a function of standard uncertainties in each pair's two members. Pair members and line and point signatures are as indicated in Fig. 4. Plots are shown only for the three pairs where pH uncertainty is not an input variable.</figcaption></figure>
<figure><img src="image_13.png" /><figcaption>Fig. 9. Error-space diagrams of the percent relative combined standard uncertainty in [H⁺] computed from three input pairs as a function of standard uncertainties in each pair’s two members. Pair members and line and point signatures are as indicated in Fig. 4. Plots are shown only for the three pairs where computed [H⁺] depends on both members of the input pair.</figcaption></figure>

$$u_c^2(y) = e^2(x_1) + e^2(x_2) + 2r(x_1, x_2) e(x_1) e(x_2) \quad (6)$$

where $e(x_i) = (\partial y/\partial x_i)u(x_i)$, a definition that combines each input sensitivity and standard uncertainty $u(x_i)$ into a term with the same units as $u_c(y)$. Thus the three terms on the right side of Eq. (6) have the same units as $u_c^2(y)$. With the same approach for Eq. (2), we assessed contributions to the total propagated uncertainty $u_c$ from uncertainties in each member of the input pair (var1 and var2) and all constants grouped together but separated into

systematic (sys) and random (ran) components, all of which may be added in quadrature,

$$u_c^2(y) = e^2(\text{var } 1) + e^2(\text{var } 2) + e^2(K_{\text{sys}}) + e^2(K_{\text{ran}}) \quad (7)$$

Given our standard uncertainties listed in Table 1, most of the $u_c$ in calculated variables usually derives from the standard uncertainties in the equilibrium constants, with much of that coming from their systematic

<figure><img src="image_14.png" /><figcaption>Fig. 10. Error-space diagrams of the percent relative combined standard uncertainty in [HCO₃⁻] computed from six input pairs as a function of standard uncertainties in each pair’s two members. Pair members and line and point signatures are as indicated in Fig. 4.</figcaption></figure>
<figure><img src="image_15.png" /><figcaption>Fig. 11. Error-space diagrams of the percent relative combined standard uncertainty in the [HCO₃⁻]/[H⁺] ratio computed from six input pairs as a function of standard uncertainties in each pair's two members. Pair members and line and point signatures are as indicated in Fig. 4.</figcaption></figure>

component (Fig. 12). There are exceptions when pCO₂ is computed from pH combined with either A_T or C_T and when [CO₃²⁻] or Ω_A is computed from the pCO₂-pH pair. In those cases, the largest contribution to overall uncertainty comes from our total standard uncertainty in pH (0.01); conversely, if we considered only the random component of standard uncertainty in pH (0.003), its contribution to the overall combined standard uncertainty would remain minor and propagated uncertainties from the equilibrium constants would again dominate. Thus with standard uncertainties for state-of-the art measurement techniques, in most cases much less than half of the combined standard uncertainty derives from the joint contribution from the members of the input pair. However, doubling those standard uncertainties would quadruple the size of the input pairs' bar segments, while segments for the constants would not change. In some cases, that would cause the contributions from the standard uncertainties of the input pair to dominate the combined standard uncertainty. That dominance would increase even more when using the largest standard uncertainties from Dickson (2010) for members of the input pair, which are three to five times higher than our default standard uncertainties. The contribution to combined standard uncertainty from the input pair comes largely from pCO₂ or pH if one of those is a member and from pH for the pCO₂-pH input pair. For the A_T-C_T input pair, the contributions are more balanced. These stacked histograms reveal the different contributions to the combined standard uncertainty but are less useful to study how changing the standard uncertainty from the input pair affects the combined standard uncertainty (square root of each y axis in Fig. 12), for which is better to refer back to the error-space diagrams (Figs. 4–11).

Second, we distinguished the contribution from the individual constants to the total from all constants

$$e^2(K) = e^2(K_0) + e^2(K_1) + e^2(K_2) + e^2(K_B) + e^2(K_A) + \dots \quad (8)$$

where the e²(K) is also the sum of the last 2 terms in Eq. (7) and thus has the same units as u_c²(y). Contributions from unlisted equilibrium constants (K_W, K_Si, K_1P, K_2P, K_3P) are insignificant. The contribution from K_A concerns only Ω_A, and that is replaced by the contribution from K_C when calculating Ω_C. The constants contribute in different proportions to e²(K), depending on the output variable and the input pair. Two constants generally dominate, two others are mostly negligible, and two others contribute about half of e²(K) but only for one computed variable each. Contributions to e²(K) from e²(K₁) or e²(K₂) usually dominate those from other constants (Fig. 13). The contribution of e²(K₁) dominates other contributions to u_c²(pCO₂) from e²(K) for input pairs having pH as one of its members; likewise, e²(K₁) dominates the contributions to u_c²(pH) from e²(K) for input pairs having pCO₂ as a member. But the contribution from e²(K₁) is second in importance to e²(K) for u_c²([CO₃²⁻]) when pCO₂ is a member of the input pair. In those cases, it is e²(K₂) that dominates other contributions to e²(K) when computing [CO₃²⁻] with any input pair. Likewise e²(K₂) contributes most to e²(K) for pCO₂ and [H⁺] when those two variables are computed from the A_T-C_T input pair.

Two other equilibrium constants generally contribute little. The contribution from e(K₀) is null except when pCO₂ is either an input variable or calculated variable, and even then its contribution is negligible. Slightly more important is e(K_B), which makes a small contribution to combined standard uncertainties of variables calculated with the A_T-C_T input pair. Although small, these contributions from e(K_B) outweigh those from e(K₁) when the same input pair is used to calculate [H⁺], [CO₃²⁻], and Ω_A; conversely for calculated pCO₂, the
<figure><img src="image_16.png" /><figcaption>**Fig. 12.** Stacked bar chart showing $u_c^2$ of (a) pCO₂, (b) [H⁺], (c) [CO₃²⁻], and (d) ΩA separated into the components derived from the standard uncertainties in each member of each input pair (var1 and var2) and the joint contribution from the standard uncertainties of the equilibrium constants, separated into random (ran) and systematic (sys) contributions. Only three input pairs are shown in the top row because pCO₂ or pH are members of the unshown input pairs. Input pair labels are given on the x axis in the order of var1-var2. Input conditions are for average Southern Ocean surface waters: A_T = 2295 µmol kg⁻¹, C_T = 2155 µmol kg⁻¹, T = -0.49 C, S = 33.96, P_T = 2.0 µmol kg⁻¹, and Si_T = 60 µmol kg⁻¹ as well as corresponding calculated pH (8.114) and pCO₂ (324.8 µatm). Standard uncertainties are from the Total column in Table 1.</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>**Fig. 13.** Stacked bar chart showing contributions from individual constants to the $u_c^2$ of (a) pCO₂, (b) [H⁺], (c) [CO₃²⁻], and (d) ΩA. Reds and blues indicate the total contribution (random + systematic) from each of K₀, K₁, K₂, K₃, and K₄ (the solubility product of aragonite). Green indicates the contribution from the 2% default standard uncertainty in the B_T/S ratio. Contributions from other constants are negligible. Input conditions are as in Fig. 12. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>
contribution from e(K₁) is larger. But what stands out most relative to contributions from e(K₁) and e(K₂) are those from e(Kₐ) and e(Kc), although only for uₑ(Ωₐ) and uₑ(Ωc), respectively. These large contributions make the percent relative combined standard uncertainty for Ωₐ substantially larger than that for [CO₃²⁻]. The e(Kₐ) contribution to uₑ(Ωₐ) essentially adds to the same general pattern of contributions to e(K) seen for uₑ([CO₃²⁻]), i.e., from e(K₂), e(K₁) and e(Kₐ). The contribution of e(Kₐ) to the e(K) component of uₑ(Ωₐ) dominates for all input pairs. The contribution of e(Kc) to uₑ(Ωc) is just like that for e(Kₐ) to uₑ(Ωₐ). Hence the standard uncertainties in Kₐ and Kc should not be neglected when propagating uncertainties in Ωₐ and Ωc.

In addition, Fig. 13 shows the contribution of the uncertainty for the total-boron-to-salinity ratio (e²(BT)), a term that is included but not explicitly shown on the right-hand side of Eq. (8). That contribution is negligible except when computing pCO₂, [H⁺], and [CO₃²⁻] from the Aₜ-Cₜ pair, and even then it is small, comparable in magnitude to the contribution from Kₐ. If the default 2% uncertainty in e(BT) is doubled, its contribution to e²(K) would quadruple, but it would still contribute little to the total and be outweighed by the default contributions from the other constants. However, the uncertainty in the total-boron-to-salinity ratio remains poorly constrained as does its effect on the Mehrbach et al. (1973) equilibrium constants.

### 4.3. Covariance

Previous efforts to propagate uncertainties for the marine carbonate system have neglected covariances between the uncertainties of the measured variables. With the new uncertainty-propagation software developed here, users may specify the correlation *r* between uncertainties in the two members of the input pair but not other correlations. The reason is that even that correlation is difficult to assess. By choice, we do not provide a general covariance matrix. That matrix would depend on many factors, being situation dependent. At best, measurements of each member of the input pair come with an estimate of standard uncertainty that is a constant, either absolute or relative, with no a priori concept that it varies and much less that it covaries with the uncertainty of the other member of the input pair. Correlations between uncertainties in the equilibrium constants would also need to be included in a complete covariance matrix. The equilibrium constants are all functions of *T* and *S* so that the constants themselves will be correlated as *T* and *S* vary. Yet when standard uncertainties in *T* and *S* are neglected, as adopted for our analysis, there is no covariance between constants generated by their relationship with *T* and *S*. If *T* and *S* input uncertainties were larger, one could in theory come up with a covariance matrix between variables (*T*, *S*, and the pertinent equilibrium constants), although that would be laborious. But such would not suffice because to propagate uncertainties, we need to specify covariances between the uncertainties in input variables, not between the input variables themselves.

We can think of a few cases where standard uncertainties might have been correlated between members of the input pair or between constants. For instance, Aₜ and Cₜ were measured using a single titration during the Geochemical Ocean Sections Study (GEOSECS) (Bradshaw et al., 1981) and the Transient Tracers in the Ocean (TTO) program's North Atlantic Study (Bradshaw and Brewer, 1988). Likewise, some studies have measured both K₁ and K₂ with a single titration (Hansson, 1973; Millero et al., 2006). Conversely, the measurements of K₁ and K₂ from Mehrbach et al. (1973) refitted by Lueker et al. (2000) and recommended for best practices (Dickson et al., 2007; Dickson, 2010) were made separately. Nonetheless Mehrbach et al. measured K₁ and the product K₁K₂, later separating out K₂, so there could still be some correlation. But in all cases it is difficult to quantify potential covariance. Nowadays any two carbonate system input variables (Aₜ, Cₜ, pCO₂, and pH) are seldom measured with the same instrument. We are aware of just two exceptions: (1) where pCO₂ and Cₜ are measured with same infra-red analyzer (see, e.g., Barton et al., 2012) in an instrument referred to as the Burke-o-Lator (https://www.hakai.org/blog/life-at-hakai/meet-burke-o-lator) and (2) where Aₜ, Cₜ, and pH are all measured during a single titration (https://www.locean-ipsl.upmc.fr/

SNAPO).

### 4.4. GOA-ON assessment

Our analysis illustrates how GOA-ON's Weather goal of a 10% percent relative combined standard uncertainty when [CO₃²⁻] is computed from point measurements is achievable even with the largest measurement uncertainties listed in Table 1 except when using the pCO₂-pH pair (Fig. 4). Likewise, a 10% relative uncertainty in Ωₐ is attainable with the same large input uncertainties for all input pairs except pCO₂-pH, despite overall propagated uncertainties for Ωₐ being at least 2% larger than those for [CO₃²⁻] owing to the uncertainty in the solubility product of aragonite Kₐ (Fig. 5).

In contrast, the more stringent GOA-ON Climate goal referring to a 1% relative combined standard uncertainty is defined not in terms of point measurements but rather for a difference in computed [CO₃²⁻] at two different times (a trend). It is framed in that way so that some of the contributing uncertainties will cancel. Thus the systematic components of the standard uncertainties of the constants are neglected, being assumed identical for both measurements. The random components of the standard uncertainties of the constants are also assumed to cancel, considering that differences in time are due only to a change in Cₜ while *T* and *S* remain unchanged. However, if *T* or *S* would change substantially between the two points in time, so would the constants and their absolute uncertainties. If uncertainties from the constants do not entirely cancel, the Climate goal may be too stringent. For example, using only the random uncertainties from the constants, the resulting propagated uncertainty in a single [CO₃²⁻] estimate from the best input pair (Aₜ-Cₜ) would be about 1% in the tropics and 1.3% in the Southern Ocean (Fig. 14). The corresponding propagated combined standard uncertainties for a difference (adding in quadrature) would be 1.4% and 1.8%, respectively. Those uncertainties would of course increase after accounting for measurement uncertainties from the input pair. Although the Weather goal is easily attained, the Climate goal may be ambiguous when climate change and variability are prominent. Therefore, the GOA-ON community should consider also setting a less ambitious intermediate goal, one that would focus on individual measurements to calculate Ωₐ as well as [CO₃²⁻]. A higher-end goal for individual measurements is needed, beyond the Weather goal, e.g., to accurately establish when waters cross chemical thresholds such as Ωₐ = 1, a concern that is not directly addressed by GOA-ON's Climate goal.

### 4.5. Neglected uncertainties

When propagating uncertainties, the four packages rely on the same total alkalinity equation,

$$A_T = [\text{HCO}_3^-] + 2[\text{CO}_3^{2-}] + [\text{B(OH)}_4^-] + [\text{OH}^-] - [\text{H}^+]_F - [\text{HSO}_4^-] - [\text{HF}] + [\text{HPO}_4^{2-}] + 2[\text{PO}_4^{3-}] - [\text{H}_3\text{PO}_4] + [\text{SiO(OH)}_3^-] \quad (9)$$

which in some cases does not include all components that make a significant contribution to the measured Aₜ. For example, seawater may contain non-negligible contributions to total alkalinity from organic acid systems (e.g., Cai et al., 1998; Yang et al., 2015), a missing component in all packages. Likewise, anoxic and hypoxic zones contain NH₃ and HS⁻, both of which contribute to measured Aₜ but are neglected in Eq. (9). Hence this low-oxygen alkalinity is not considered in our uncertainty propagation add-ons, which rely on the standard routines to compute derived variables in the four packages, each of which uses Eq. (9).

Although lacking these components, the uncertainty propagation software developed here can still be used to propagate the associated uncertainties. For example, considering well oxygenated seawater that contains organic-acid alkalinity (one additional alkalinity component Aᵢ), its quantity can be estimated from

$$A_i = A'_i - A_{\text{calc}}^i \quad (10)$$

where A'ᵢ is the measured total alkalinity and Aᵢ³calc is calculated from
<figure><img src="image_18.png" /><figcaption>Fig. 14. Percent relative combined standard uncertainty in [CO₃²⁻] as a function of standard uncertainties in C_T and A_T, which are assumed equal. Results are shown for surface waters of (left) the Southern Ocean and (right) the tropics, using their average conditions specified in Fig. 4. The total combined standard uncertainty (red line) accounts for uncertainties from the measurements (black line) and systematic and random uncertainties from the constants (blue line). Also shown is the combined standard uncertainty component calculated from only the random uncertainties in the constants (Table 1, cyan line). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

simultaneous measurements of an input pair such as pH-C_T. The associated standard uncertainty u(A_i) could be estimated from

$$u(A_i) = \sqrt{u^2(A'_T) + u^2(A_T^{calc})} \qquad (11)$$

where u(A'_T) is the estimated standard uncertainty of a total alkalinity measurement (probably as in Table 1) and u(A_T^{calc}) is the combined standard uncertainty in total alkalinity calculated from pH and C_T.

If the goal is to estimate a more complete uncertainty for a total alkalinity measurement that is being treated as though Eq. (9) holds true and the organic-acid alkalinity is being systematically ignored, then it would be appropriate to propagate the following uncertainty:

$$u(A_T) = \sqrt{u^2(A'_T) + u^2(A_i)} = \sqrt{2u^2(A'_T) + u^2(A_T^{calc})} \qquad (12)$$

Finally, the resulting u(A_T) could be used as input to the uncertainty propagation software introduced here. In low oxygen systems, the approach could be expanded assuming that concentrations of NH₃ and HS⁻ are known along with estimates of their uncertainties. Ultimately, it would be better to add these various missing components directly to the uncertainty propagation software's total alkalinity equation, e.g., by building on recent improvements to account for low-oxygen alkalinity when computing derived variables in CO2SYS-Excel (Xu et al., 2017) and seacarb (carbfull function from M. Hagens). This additional effort is left for future work.

These neglected components in the A_T equation only affect propagated uncertainties when A_T is calculated as a derived variable or when it is a member of the input pair. One may estimate the effect of these systematic uncertainties from the error-space diagrams presented in Sect. 3.2. As illustrated, with the A_T-C_T input pair, combined standard uncertainties in calculated pCO₂, pH, [CO₃²⁻], and Ω_A are much more sensitive to uncertainties in A_T than when they are computed with the pH-A_T and pCO₂-A_T input pairs.

## 5. Conclusions

We have written software to propagate uncertainties and provided it as add-ons to four commonly used public packages that compute marine carbonate chemistry. These packages include CO2SYS-Excel (Visual Basic), CO2SYS-MATLAB (MATLAB), seacarb (R), and mocsy (Fortran). The last three can be used in interactive Jupyter notebooks, while the last two can be used directly in Python (Appendix D). Results agree across packages and with analytical solutions (for the intermediate sensitivities) to better than 0.01%. Given the consistent implementation across packages, users may use any of the four updated packages to compute accurate numerical buffer factors and propagate

uncertainties. Uncertainties can now be propagated routinely as easily as derived variables can be computed, whether working with pre-existing or new data sets. There remains no good reason for marine scientists not to compute and report propagated uncertainties along with their computed carbonate system variables.

Along with this software, we provide here a new kind of diagram to graphically estimate combined standard uncertainties quickly. These error-space diagrams are useful to compare combined standard uncertainties among input pairs across the spectrum of standard uncertainties that may be used as input and to estimate the implications of changes in these standard uncertainties. These diagrams further demonstrate that often, reducing the standard uncertainties of the measurements beyond the current state of the art would not substantially reduce the propagated combined standard uncertainties. Conversely, relaxing these measurement constraints by a factor of two or more often does not substantially increase combined standard uncertainty. Substantial reduction in combined standard uncertainties below values obtained with today's best measurement uncertainties for the input pair can only be realized by reducing the standard uncertainties of the equilibrium constants, particularly K₁ and K₂. The standard uncertainties for the equilibrium constants and hence their contribution to the propagated combined standard uncertainties remain poorly understood. Lack of consideration of the systematic components of these uncertainties, or even neglect of uncertainties from the equilibrium constants altogether, has led some previous studies to underestimate combined standard uncertainty. Future work should focus on improving estimates of the standard uncertainties in the equilibrium constants, particularly in regard to how they change with pressure. Generally, it has not been appreciated that the relative combined standard uncertainties for Ω_A and Ω_C are substantially larger than that for [CO₃²⁻] because of the large standard uncertainties in the solubility products for aragonite and calcite.

It is hoped that the new software and the new type of diagram proposed here will be used widely and expanded upon so that propagating uncertainties will become standard practice when calculating marine CO₂ system variables. Although the convenience of these tools offers a step forward, more is still required. The uncertainties from the various constants need to be evaluated more rigorously as emphasized above. Scientists will also need to assess their measurement uncertainties more carefully than has usually been the case. To assist with this fundamental effort, there are two valuable resources: the Eurachem Guide on quantifying uncertainty in analytical measurement (Ellison and Williams, 2012) and a software package known as the GUM Tree Calculator (Hall, 2013; Lovell-Smith et al., 2017), which allows users to propagate uncertainties for any arbitrary function. Widespread use of
all of these resources will allow the community to transition from one that seldom propagates uncertainties to one where uncertainties are routinely propagated following best practices from metrology.

## Acknowledgments

This effort was proposed by the SOLAS-IMBER Ocean Acidification Working group (SIOA) and the Advisory Board for the Ocean Acidification International Coordination Centre (OA-ICC) as an effort to promote best practices. It was funded by OA-ICC, a project of the Peaceful Uses Initiative

of the International Atomic Energy Agency (IAEA), which supported the consultancy of JME. We thank Lina Hansson for managing this project within the OA-ICC and Cathy Nanini for developing the interactive interface that produces error-space diagrams. JCO was supported by the EU H2020 Project CRESCENDO (grant #641846) and the French ANR project SOBUMS. AGD was supported, in part, by the US National Science Foundation under grants OCE-1233648 and OCE-1657799. JPG was supported by the EU-H2020 project INTAROS (grant #727890) and IPEV, the French Polar Institute.

## Appendix A. Matrix form of uncertainty propagation

To streamline the uncertainty propagation equation when dealing with multiple input and output variables, it may be written in terms of matrices. When only one variable y is computed from n input variables (x₁, x₂, ..., xₙ), its combined variance uₑ²(y) is given by multiplication of 2 vectors (1-D matrices) and 1 matrix

$$u_c^2(y) \approx \gamma C_{xx} \gamma^T \quad (A.1)$$

where γ is the row vector of partial derivatives

$$\gamma = \left[ \frac{\partial y}{\partial x_1} \quad \frac{\partial y}{\partial x_2} \quad \frac{\partial y}{\partial x_3} \quad \dots \quad \frac{\partial y}{\partial x_n} \right] \quad (A.2)$$

γT is the corresponding transpose (a column vector), and Cxx is the covariance matrix of the input variables

$$C_{xx} = \begin{bmatrix} u^2(x_1) & u(x_1 x_2) & u(x_1 x_3) & \dots & u(x_1 x_n) \\ u(x_2, x_1) & u^2(x_2) & u(x_2, x_3) & \dots & u(x_2, x_n) \\ u(x_3, x_1) & u(x_3, x_2) & u^2(x_3) & \dots & u(x_3, x_n) \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ u(x_n, x_1) & u(x_n, x_2) & u(x_n, x_3) & \dots & u^2(x_n) \end{bmatrix} \quad (A.3)$$

Eq. (A.1) is thus the matrix form of Eq. (2). More generally with m dependent variables (y₁, y₂, ..., yₘ) computed from the same set of n independent variables (x₁, x₂, ..., xₙ), the combined variances of those dependent variables (uₑ²(y₁), uₑ²(y₂), ..., uₑ²(yₘ)) are simply the diagonal terms of covariance matrix of the computed results

$$C_{yy} = J_{yx} C_{xx} J_{yx}^T \quad (A.4)$$

where Cxx is defined in Eq. (A.3), Cyy is the analogous covariance matrix except that it is for the output variables, and Jyx is the Jacobian matrix, i.e., the m × n array of first derivatives

$$J_{yx} = \begin{bmatrix} \frac{\partial y_1}{\partial x_1} & \dots & \frac{\partial y_1}{\partial x_n} \\ \vdots & \ddots & \vdots \\ \frac{\partial y_m}{\partial x_1} & \dots & \frac{\partial y_m}{\partial x_n} \end{bmatrix} \quad (A.5)$$

A more detailed development is provided by Arras (1998).

## Appendix B. Analytical solutions for partial derivatives

To evaluate the numerical solutions of the partial derivatives, they were compared to available analytical solutions. For the analytical solutions, we used the six buffer factors proposed by Egleston et al. (2010), Sabine and Morel], using known corrections and after modifying them to also account for alkalinity from phosphoric and silicic acid systems. These buffer factors are

$$\gamma_{AT} = \left( \frac{\partial \ln[CO_2^*]}{\partial A_T} \right)^{-1} \quad (B.1a)$$

$$\gamma_{CT} = \left( \frac{\partial \ln[CO_2^*]}{\partial C_T} \right)^{-1} \quad (B.1b)$$

$$\beta_{AT} = \left( \frac{\partial \ln[H^+]}{\partial A_T} \right)^{-1} \quad (B.1c)$$

$$\beta_{CT} = \left( \frac{\partial \ln[H^+]}{\partial C_T} \right)^{-1} \quad (B.1d)$$

$$\omega_{AT} = \left( \frac{\partial \ln\Omega}{\partial A_T} \right)^{-1} \quad (B.1e)$$
$$ \omega_{C_T} = \left( \frac{\partial \ln \Omega}{\partial C_T} \right)^{-1} \qquad (B.1f) $$

Each of them are calculated with formulas that involve the $S_e$ term, either directly or indirectly, where

$$ S_e = [\text{HCO}_3^-] + 4[\text{CO}_3^{2-}] + \frac{[\text{H}^+] [\text{B(OH)}_4^-]}{K_B + [\text{H}^+]} + [\text{H}^+] + [\text{OH}^-] \qquad (B.2) $$

and the final term on the right-hand side has been corrected for the sign error in (Egleston et al.'s Eq. A27) [Álvarez et al., 2014; Orr, 2011). The key here is to realize that the full right-hand side of Eq. (B.2) may also be obtained from the alkalinity equation, written in terms of only equilibrium constants, $[H^+]$, and $[CO_2^*]$ (Egleston et al., 2010, Eq.10), i.e., by taking the derivative of both of its sides with respect to $[H^+]$ when $[CO_2^*]$ is constant. Thus,

$$ S_e = - \left( \frac{\partial A_T}{\partial \ln[H^+]} \right)_{CO_2^*} = -[H^+] \left( \frac{\partial A_T}{\partial [H^+]} \right)_{CO_2^*} \qquad (B.3) $$

where for $A_T$ Egleston et al. considered contributions only from carbonic and boric acid systems and the dissociation of water. Thus we can now expand the Egleston et al. buffer factors to also take into account the alkalinity from silicic and phosphoric acid systems. To do so, let us first identify the pertinent terms to be added to the alkalinity equation, namely $[\text{Si(OH)}_3^-] + [\text{HPO}_4^{2-}] + 2[\text{PO}_4^{3-}] - [\text{H}_3\text{PO}_4]$. Then we need to take their derivatives with respect to $[H^+]$, adding those to the right-hand side of Eq. (B.2). Thus

$$ S_e = [\text{HCO}_3^-] + 4[\text{CO}_3^{2-}] + \frac{[\text{H}^+] [\text{B(OH)}_4^-]}{K_B + [\text{H}^+]} + [\text{H}^+] + [\text{OH}^-] \\ + ([\text{H}_3\text{PO}_4] ([\text{H}_2\text{PO}_4^-] + 2[\text{HPO}_4^{2-}] + 3[\text{PO}_4^{3-}]) \\ + [\text{HPO}_4^{2-}] (2[\text{H}_3\text{PO}_4] + [\text{H}_2\text{PO}_4^-] - [\text{PO}_4^{3-}]) \\ + 2[\text{PO}_4^{3-}] (3[\text{H}_3\text{PO}_4] + 2[\text{H}_2\text{PO}_4^-] + [\text{HPO}_4^{2-}]))/P_T \\ + \frac{[\text{H}^+] [\text{SiO(OH)}_3^-]}{K_{\text{Si}} + [\text{H}^+]} \qquad (B.4) $$

where the first line is Eq. (B.2), lines 2–4 account for the added contributions from the phosphoric acid system, using the analytical solution from Hagens and Middelburg (2016), and the final line is for the addition from the simplified silicic acid system (Dickson et al., 2007), whose contribution has the same form as the term for the boric acid system in line 1.

Besides $S_e$, other published formulas from Egleston et al. (2010) had typographical errors. They too were previously corrected by comparing the published formulas (in Egleston et al.'s Table 1) to those in the authors' spreadsheet (C. Sabine, personal communication). Two subsequent studies detailed the same corrections independently: Álvarez et al. (2014) published the corrected equations in their Table 3; Orr (2011) used the corrected equations while providing them in a subroutine that computes these buffer factors (buffesm in the seacarb package). The additional corrections concern two of the six buffer factors. The corrected equations are

$$ \omega_{C_T} = C_T - A_C \frac{P_e}{Q_e} \qquad (B.5a) $$

$$ \omega_{A_T} = A_C - C_T \frac{Q_e}{P_e} \qquad (B.5b) $$

where

$$ P_e = 2[\text{CO}_2^*] + [\text{HCO}_3^-] \qquad (B.6a) $$

$$ Q_e = [\text{HCO}_3^-] + \frac{[\text{H}^+] [\text{B(OH)}_4^-]}{K_B + [\text{H}^+]} - [\text{H}^+] - [\text{OH}^-] \qquad (B.6b) $$

For both terms we use the same nomenclature as in the spreadsheet from Egleston et al. (2010); those terms were partially confused in the published formulas. In addition, we note that Eq. (B.6b) simplifies to

$$ Q_e = 2A_C - S_e \qquad (B.7) $$

Hence $Q_e$ also changes when accounting for alkalinity contributions from phosphoric and silicic acid systems. Thus all six buffer factors depend on $S_e$, which has been modified here to depend on $P_T$ and $Si_T$.

These buffer factors now only need to be converted to absolute partial derivatives to be used in the uncertainty propagation equation. For example, knowing that $\partial \ln x = \partial x/x$, it follows that

$$ \frac{\partial[\text{CO}_2^*]}{\partial A_T} = [\text{CO}_2^*] \frac{\partial \ln[\text{CO}_2^*]}{\partial A_T} = \frac{[\text{CO}_2^*]}{\gamma_{A_T}} \qquad (B.8) $$

Likewise, the same procedure is used to obtain all other partial derivatives with respect to $A_T$

$$ \frac{\partial p\text{CO}_2^*}{\partial A_T} = \frac{p\text{CO}_2^*}{\gamma_{A_T}} \qquad (B.9a) $$

$$ \frac{\partial [\text{H}^+]}{\partial A_T} = \frac{[\text{H}^+]}{\beta_{A_T}} \qquad (B.9b) $$
$$ \frac{\partial \Omega_A}{\partial A_T} = \frac{\Omega_A}{\omega_{A_T}} \qquad (B.9c) $$

$$ \frac{\partial \Omega_C}{\partial A_T} = \frac{\Omega_C}{\omega_{A_T}} \qquad (B.9d) $$

$$ \frac{\partial[\text{CO}_3^{2-}]}{\partial A_T} = \frac{[\text{CO}_3^{2-}]}{\omega_{A_T}} \qquad (B.9e) $$

For the analogous partial derivatives with respect to $C_T$, one only needs to change the subscripts of $\gamma$, $\beta$, and $\omega$ from $A_T$ to $C_T$. Partial derivatives of $[\text{HCO}_3^-]$ are not as readily accessible but may be calculated from others. Knowing that $[\text{HCO}_3^-] = [\text{CO}_3^{2-}][\text{H}^+]/K_2$, let us modify each side of that equation by taking the logarithm and then taking the derivative with respect to $A_T$, which yields

$$ \frac{\partial \ln[\text{HCO}_3^-]}{\partial A_T} = \frac{\partial \ln[\text{CO}_3^{2-}]}{\partial A_T} + \frac{\partial \ln[\text{H}^+]}{\partial A_T} \qquad (B.10) $$

That is the same as

$$ \frac{1}{\delta_{A_T}} = \frac{1}{\omega_{A_T}} + \frac{1}{\beta_{A_T}} \qquad (B.11) $$

where

$$ \delta A_T = \left( \frac{\partial \ln[\text{HCO}_3^-]}{\partial A_T} \right)^{-1} \qquad (B.12) $$

Hence as for Eq. (B.8), the corresponding partial derivative is

$$ \frac{\partial[\text{HCO}_3^-]}{\partial A_T} = \frac{[\text{HCO}_3^-]}{\delta_{A_T}} \qquad (B.13) $$

The same approach is taken to derive $\partial[\text{HCO}_3^-]/\partial C_T$.

## Appendix C. Error-space diagram

As described in Sect. 2.5, an error-space diagram combines three principal objects: (1) a basic contour plot of the combined standard uncertainty as a function of standard uncertainties from the two members of the input pair, (2) a *pair-constants curve* indicating where total propagated uncertainties from the constants are equal to total propagated uncertainties from the input pair, (3) a *pair line* indicating where the contribution from the uncertainty of each member of the input pair contributes equally to the combined standard uncertainty. Here we detail how the second and third objects are computed.

### C.1. *Pair-constants curve*

Let us define a *pair-constants curve* along which the contribution to the combined variance from both members of input pair is equal to that from all the constants:

$$ e^2(x_1) + e^2(x_2) = \sum_i e^2(K_i) \qquad (C.1) $$

We can simultaneously generate pairs of coordinates $u(x_1)$ and $u(x_2)$ along this *pair-constants curve* by rearranging Eq. (C.1) to

$$ \frac{e^2(x_1)}{\sum_i e^2(K_i)} + \frac{e^2(x_2)}{\sum_i e^2(K_i)} = 1 \qquad (C.2) $$

which may then be considered in terms of the Pythagorean identity, $\sin^2\theta + \cos^2\theta = 1$. That is, for angle $\theta$ of a right triangle, $e(x_1)$ is the adjacent side (x), $e(x_2)$ is the opposite side (y), and $(\sum e^2(K_i))^{0.5}$ is the hypotenuse (r). Given the coordinates of the Pythagorean identity ($x = r \cos \theta$ and $y = r \sin \theta$), we can find the corresponding coordinates for the *pair-constants curve* by plugging in the three uncertainty components and expanding $e(x_1)$ and $e(x_2)$, yielding

$$ u(x_1) = \left( \sum_i e^2(K_i) \right)^{1/2} \cos \theta \left| \frac{\partial y}{\partial x_1} \right|^{-1} \qquad (C.3a) $$

$$ u(x_2) = \left( \sum_i e^2(K_i) \right)^{1/2} \sin \theta \left| \frac{\partial y}{\partial x_2} \right|^{-1} \qquad (C.3b) $$

For a plot of $u(x_1)$ vs. $u(x_2)$, the *pair-constants curve* is computed by varying $\theta$ from 0 to $\pi/2$ radians. Along that curve, there are three points of particular interest. First is its *pair-constants midpoint* where uncertainties from each member of the input pair are equivalent

$$ e^2(x_1) = e^2(x_2) = \frac{1}{2} \sum_i e^2(K_i) \qquad (C.4) $$

occurring at $\sin \theta = \cos \theta$, i.e., at $\pi/4$. Second and third are its x and y intercepts, occurring when the total uncertainty from the input pair comes solely from the first member (at $\theta = 0$, where $u(x_2) = 0$) and when it comes solely from the second member (at $\theta = \pi/2$, where $u(x_1) = 0$).
## C.2. Pair line

More generally, let us determine the exact point on any isoline (for the propagated uncertainty in the computed variable) where $e(x_1) = e(x_2)$. This midpoint can be defined from Eq. (7) after combining terms for the constants, rearranging, and simplifying names, to be

$$e^2(x_1) = e^2(x_2) = \frac{1}{2} \left( u_c^2(y) - \sum_i e^2(K_i) \right) \qquad (C.5)$$

Solving separately for $u(x_1)$ and $u(x_2)$, we have

$$u_1^2(x_1) = \frac{1}{2} \left( u_c^2(y) - \sum_i e^2(K_i) \right) \left( \frac{\partial y}{\partial x_1} \right)^{-2} \qquad (C.6a)$$

$$u_2^2(x_2) = \frac{1}{2} \left( u_c^2(y) - \sum_i e^2(K_i) \right) \left( \frac{\partial y}{\partial x_2} \right)^{-2} \qquad (C.6b)$$

In Eq. (C.6a,b), $u(x_1)$ and $u(x_2)$ each have a single value per isoline, the pair of which defines the isoline's midpoint. That midpoint on the pair-constants curve is the same as the pair-constants midpoint, i.e., where

$$u_c^2(y) = 2 \sum_i e^2(K_i) \qquad (C.7)$$

Thus on the pair-constants curve, Eq. (C.6a,b) becomes Eq. (C.3a,b) for $\theta = \pi/4$. The pair line is determined by varying $u_c(y)$ to compute multiple midpoints.

# Appendix D. Code availability and documentation

The add-on routines for calculating sensitivities (derivnum) and propagating uncertainties (errors) can be downloaded from CRAN for seacarb and from GitHub for CO2SYS-Excel, CO2SYS-MATLAB, and mocsy, as detailed below. These sites include not only the add-ons but the complete packages, including documentation and examples. The new add-ons follow the initial design of each package's approach to calculate carbonate chemistry variables.

## D.1. CO2SYS-Excel

For CO2SYS-Excel, the modified software may be downloaded from https://github.com/jamesorr/CO2SYS-Excel. In CO2SYS-Excel, a fourth sheet (tab) was added to specify input uncertainties and display output uncertainties once the “Start” button is clicked on third sheet. Calculated variables are, as before, displayed on the third sheet, while their combined standard uncertainties are provided on the new fourth sheet, which maintains the same layout and units as the third sheet. Additional documentation has been added to the preexisting first tab and to the new fourth tab in order to allow users to quickly make the transition to routinely computing and reporting propagated uncertainties.

## D.2. CO2SYS-MATLAB

For CO2SYS-MATLAB users, the new add-ons and the slightly modified CO2SYS.m routine can be downloaded from https://github.com/jamesorr/CO2SYS-MATLAB. In that repository, the src, examples, and notebooks subdirectories include the source code, examples as MATLAB scripts, and interactive examples as Jupyter notebooks. The errors routine is called as follows:

```
errors (par1, par2, par1type, par2type, sal, tempin, tempout, presin, presout, si, po4, ...
        epar1, epar2, esal, etemp, esi, epo4, epk, ebt, r, ...
        pHSCALEIN, K1K2CONSTANTS, KSO4CONSTANTS);
```

The first and third lines contain the arguments that are identical to those used with the CO2SYS base routine CO2SYS.m. Conversely, the second line is used to enter the standard uncertainties for each member of the input pair, $S, T, Si_T$, and $P_T$ (same units as in line 1), the standard uncertainties in the $pK$ values for the equilibrium constants (a vector of 7 values), the fractional relative standard uncertainty in $B_T$, and the correlation between $ePAR1$ and $ePAR2$. The derivnum routine is called as follows:

```
derivnum (varid, par1, par2, par1type, par2type, sal, ...
        tempin, tempout, presin, presout, si, po4, ...
        pHSCALEIN, K1K2CONSTANTS, KSO4CONSTANTS);
```

where one argument varid has been added to the beginning of the same list of arguments as used for CO2SYS.m. That new argument is simply a character string that specifies the input variable for which partial derivatives of computed variables will be taken with respect to. More detailed documentation can be obtained with the traditional help command in Matlab. All routines also run under Octave, the GNU clone of MATLAB.

## D.3. Seacarb

For seacarb users, the new add-ons for uncertainty propagation and sensitivity calculations are included since version 3.2.4 of the complete package, which can be downloaded from https://cran.r-project.org/package/seacarb. To propagate uncertainties, seacarb's errors function is called as follows:

```
errors (flag, var1, var2, S, T, Patm, P, Pt, Sit,
        evar1, evar2, eS, eT, ePt, eSit, epK, eBt,
        k1k2, kf, ks, pHscale, b, gas)
```

Its arguments are identical to those in seacarb's carb function except a second line has been added for the user to specify the standard uncertainties in the input
variables. Optional arguments include epK (a 7-member vector) to change the default standard uncertainties for 7 equilibrium constants, eBt (a scalar) to change the default fractional relative standard uncertainty in total boron-to-salinity ratio, and r to specify the correlation between the standard uncertainties of the two members of the input pair evar1 and evar2. To calculate sensitivities, seacarb's *derivnum* function is called as follows:

```
derivnum(varid, flag, var1, var2, S, T, Patm, P, Pt, Sit,
         k1k2=, kf=, ks, pHscale, b, gas)
```

where only one argument varid has been added to the beginning of the same list of arguments as used for seacarb's carb function to define the input variable for which partial derivatives of computed variables will be taken with respect to. All arguments are detailed and examples given in seacarb's documentation at the CRAN site mentioned above. Further documentation is provided as interactive Jupyter notebooks at https://github.com/jpgattuso/seacarb.git. Those can be downloaded and used locally, after installing the software that is indicated in the HTML version of the notebook file (after clicking on the name of the file on the GitHub site). Alternatively, these notebooks can be used interactively and modified without installing local software by directing one's browser to https://mybinder.org/v2/gh/jamesorr/seacarb.git/master?filepath=notebooks%2FIRkernel and waiting for the binder application to load (usually about 1 min), an effort that is more straightforward for first-time users of Jupyter notebook.

### D.4. Mocsy

The mocsy routines including the new add-ons for uncertainty propagation and sensitivity calculations can be downloaded from https://github.com/jamesorr/mocsy/. The errors routine may be called from within Fortran or Python. In Fortran, it may be called as

```
call errors(eh, epco2, efco2, eco2, ehco3, eco3, eOmegaA, eOmegaC,
           &
           temp, sal, alk, dic, sil, phos, Patm, depth, lat, N,
           &
           temp_e, sal_e, ALK_e, DIC_e, sil_e, phos_e,
           &
           optCON, optT, optP, optB=optB, optK1K2=optK1K2, optKf=optKf )
```

where the first line of arguments is the calculated combined standard uncertainties (output), the second line is the standard input (as in mocsy's *vars* routine), the third line is for the standard uncertainties (input), and the final line is for the options (also the same as in the *vars* routine). Similarly, mocsy's *derivnum* routine is called as

```
call derivnum (dh Dx, dpcO2 Dx, dfco2 Dx, dco2 Dx, dhco3 Dx,
             &
             dco3 Dx, dOmegaA Dx, dOmegaC Dx,
             &
             temp, sal, alk, dic, sil, phos, Patm, depth, lat, N,
             &
             derivar,
             &
             optCON, optT, optP, optB=optB, optK1K2=optK1K2, optKf=optKf )
```

where the first two lines are the calculated partial derivatives (output), the fourth line is the input variable for which partial derivatives of computed variables will be taken with respect to, and the third and fifth lines are the input arguments and options as in the *errors* routine.

In Python, the mocsy *errors* routine may be called as follows:

```
[eH, epCO2, efCO2, eCO2, eHCO3, eCO3, eOmegaA, eOmegaC] = \
    mocsy.merrors(
        temp, sal, alk, dic, sil, phos, Patm, depth, lat,
        etemp, esal, ealk, edic, esil, ephos,
        optCON, optT, optP, optb, optk1k2, optkf)
```

where the first line contains the combined standard uncertainties (output), the second line is the routine specification, the third line is the standard input, the fourth line includes the standard uncertainties (input), and the fifth line includes the standard options. As with the Fortran version, users may add three optional arguments:

```
r=0.0,
    epk= np.array([0.002, 0.0075, 0.015, 0.01, 0.01, 0.02, 0.02]),
    ebt= 0.02
```

where r is for the correlation between ealk and edic, epk is for the vector of uncertainties for 7 equilibrium constants (in pK), and ebt, is for the uncertainty in total boron-to-salinity ratio (fractional units). The call in Python for mocsy's *derivnum* routine takes the following form:

```
dh Dx, dpcO2 Dx, dfco2 Dx, dco2 Dx, dhco3 Dx, dco3 Dx, domegaA Dx, domegac Dx = \
    mocsy.mderivnum(
        temp, sal, alk, dic, sil, phos, patm, depth, lat,
        derivar, optCON, optp, optb, optk1k2, optkf)
```

The mocsy archive site mentioned above also provides more detailed documentation and examples, including Fortran routines, Python scripts, and interactive Jupyter notebooks.

## References

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

Arras, K.O., 1998. An Introduction To Error Propagation: Derivation, Meaning and Examples of Equation C<sub>Y</sub> = F<sub>X</sub>C<sub>X</sub>F<sub>X</sub><sup>T</sup>. Tech. Report EPFL-ASL-TR-98-01 R3. Autonomous Systems Lab, Institute of Robotic Systems, Swiss Federal Institute of Technology Lausanne (EPFL). http://srl.informatik.uni-freiburg.de/papers/arrasTR98.pdf

Bach, L.T., 2015. Reconsidering the role of carbonate ion concentration in calcification by marine organisms. Biogeosciences 12, 4939–4951. https://doi.org/10.5194/bg-12-4939-2015.

Barton, A., Hales, B., Walbusser, G.G., Langdon, C., Feely, R.A., 2012. The pacific oyster, Crassostrea gigas, shows negative correlation to naturally elevated carbon dioxide levels: Implications for near-term ocean acidification effects. Limnol. Oceanogr. 57, 698–710. https://doi.org/10.4319/lo.2012.57.3.0698.

Bockman, E.E., Dickson, A.G., 2015. An inter-laboratory comparison assessing the quality of seawater carbon dioxide measurements. Mar. Chem. 171, 36–43.
Bradshaw, A.L., Brewer, P.G., 1988. High precision measurements of alkalinity and total carbon dioxide in seawater by potentiometric titration — 1. presence of unknown protolyte(s)? Mar. Chem. 23, 69–86.

Bradshaw, A.L., Brewer, P.G., Shafer, D.K., Williams, R.T., 1981. Measurements of total carbon dioxide and alkalinity by potentiometric titration in the GEOSECS program. Earth Planet. Sci. Lett. 55, 99–115.

Buck, R.P., Rondinini, S., Covington, A.K., Baucke, F.G.K., Brett, C.M.A., Camoes, M.F., Milton, M.J.T., Mussini, T., Naumann, R., Pratt, K.W., Spitzer, P., Wilson, G.S., 2002. Measurement of pH. Definition, standards, and procedures (IUPAC recommendations 2002). Pure Appl. Chem. 74, 2169–2200. https://doi.org/10.1351/pac200274112169.

Cai, W.J., Wang, Y., Hodson, R.E., 1998. Acid-base properties of dissolved organic matter in the estuarine waters of Georgia, USA. Geochim. Cosmochim. Acta 62, 473–483.

De Bièvre, P., 2008. Measurement uncertainty is not synonym of measurement repeatability or measurement reproducibility. Accred. Qual. Assur. 13, 61–62. https://doi.org/10.1007/s00769-008-0371-0.

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

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

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

Dickson, A.G., Sabine, C.L., Christian, J.R., 2007. Guide to Best Practices for Ocean CO₂ measurements. Technical Report. PICES Special Publication 3. pp. 191. http://aquaticcommons.org/1443/ (ISBN 1-897176-07-4).

Egleston, E., Sabine, C.L., Morel, F.M.M., 2010. Revelle revisited: Buffer factors that quantify the response of ocean chemistry to changes in DIC and alkalinity. Glob. Biogeochem. Cycles 24, GB1002. https://doi.org/10.1029/2008GB00340.

Ellison, S., Williams, A., 2012. Eurachem/CITAC Guide: Quantifying Uncertainty in Analytical Measurement, Third edn. www.eurachem.org/index.php/publications/guides/quam (ISBN 978-0-948926-30-3).

Fassbender, A.J., Alin, S.R., Feely, R.A., Sutton, A.J., Newton, J.A., Byrne, R.H., 2016a. Estimating total alkalinity in the Washington State coastal zone: complexities and surprising utility for ocean acidification research. Estuaries Coast 40, 404–418. https://doi.org/10.1007/s12237-016-0168-z.

Fassbender, A.J., Sabine, C.L., Feilé, K.M., 2016b. Consideration of coastal carbonate chemistry in understanding biological calcification. Geophys. Res. Lett. 43, 4467–4476. https://doi.org/10.1002/2016GL068860.

Frankignoule, M., 1994. A complete set of buffer factors for acid/base CO₂ system in seawater. J. Mar. Syst. 5, 111–118.

Gattuso, J.P., Epitalon, J.M., Lavigne, H., Orr, J.C., 2018. seacarb: Seawater Carbonate Chemistry with R. R Package Version 3.2.8. http://CRAN.R-project.org/package/seacarb (accessed June 26, 2018).

GUM, 1993. Guide to the Expression of Uncertainty in Measurement. ISO, Geneva. (ISBN 92-67-10188-9. (Reprinted 1995: Reissued as ISO Guide 98-3 (2008), also available from http://www.bipm.org as JCGM 100:2008)) (1993).

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

Hall, B.D., 2013. Object-oriented software for evaluating measurement uncertainty. Meas. Sci. Technol. 24, 055004. https://doi.org/10.1088/0957-0233/24/5/055004.

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

Kirchner, J., 2001. Data Analysis Toolkit #5: Uncertainty Analysis and Error Propagation. http://seismo.berkeley.edu/kirchner/eps_120/Toolkits/Toolkit_05.pdf (accessed September 3, 2016). Univ. Cal. Berkeley).

Kwiatkowski, L., Orr, J.C., 2018. Diverging seasonal extremes for ocean acidification during the twenty-first century. Nat. Clim. Chang. 8, 141–145. https://doi.org/10.1038/s41558-017-0054-0.

Lavset, S.K., Gruber, N., 2014. Long-term trends in surface ocean pH in the North Atlantic. Mar. Chem. 162, 71–76.

Lavigne, H., Gattuso, J.P., 2011. seacarb: Seawater Carbonate Chemistry with R. R package version 2.4.2. In: The Comprehensive R Archive Network. http://cran.r-project.org/package/seacarb.

Lee, K., Kim, T.W., Byrne, R.H., Millero, F.J., Feely, R.A., Liu, Y.M., 2010. 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.

Lovell-Smith, J.W., Saunders, P., Feistel, R., 2017. Unleashing empirical equations with “nonlinear fitting” and “GUM Tree Calculator”. Int. J. Thermophys. 38. https://doi.org/10.1007/s10765-017-2282-y.

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

Marion, G., Millero, F.J., Camoes, M., Spitzer, P., Feistel, R., Chen, C.T., 2011. pH of seawater. Mar. Chem. 126, 89–96.

Martz, T., Daly, K., Byrne, R., Stillman, J., Turk, D., 2015. Technology for ocean acidification research: needs and availability. Oceanography 28, 40–47. https://doi.org/10.5670/oceanog.2015.30.

McLaughlin, K., Weisberg, S.B., Dickson, A.G., Hofmann, G.E., Newton, J.A., Aseltine-Neilson, D., Barton, A., Cudd, S., Feely, R., Jefferds, I., Jewett, E., King, T., Langdon, C., McAfee, S., Pleschner-Steele, D., Steele, B., 2015. Core principles of the California current acidification network: linking chemistry, physics, and ecological effects. Oceanography 28https://doi.org/10.5670/oceanog.2015.39. 160aë~169.

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

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

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

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

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

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

Newton, J., Feely, R.A., Jewett, E.B., Williamson, P., Mathis, J., 2015. Global Ocean Acidification Observing Network: Requirements and Governance Plan, Second Edition. http://goa-on.org/documents/general/GOA-ON_2nd_edition_final.pdf (accessed June 15, 2018).

Nisumaa, A.M., Pesant, S., Bellerby, R.G.J., Delille, B., Middelburg, J.J., Orr, J.C., Riebesell, U., Tyrrell, T., Wolf-Gladrow, D., Gattuso, J.P., 2010. EPOCA/EUROCEANS data compilation on the biological and biogeochemical responses to ocean acidification. Earth Sys. Sci. Data 2, 167–175. https://doi.org/10.5194/essd-2-167-2010.

Orr, J.C., 2011. Recent and future changes in ocean carbonate chemistry. In: Gattuso, J.P., Middelburg, J. (Eds.), Ocean Acidification. Oxford Univ. Press, pp. 41–66 chapter 3.

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

Orr, J.C., Gattuso, J.P., Epitalon, J.M., 2015. Comparison of ten packages that compute ocean carbonate chemistry. Biogeosciences 12, 1483–1510. https://doi.org/10.5194/bg-12-1483-2015.

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

Pierrot, D., Lewis, E., Wallace, D.W.R., 2006. MS Excel Program Developed for CO₂ System Calculations. Technical Report. Carbon Dioxide Inf. Anal. Cent. Oak Ridge Natl. Lab., U.S. DOE, Oak Ridge, Tenn.

Proye, A., Gattuso, J.P., 2003. Seacarb, an R package to Calculate Parameters of the Seawater Carbonate System.

Skirrow, G., 1975. The dissolved gases-carbon dioxide. In: Riley, J.P., Skirrow, G. (Eds.), Chemical Oceanography, vol. 2. Academic Press, London, pp. 1–192.

Sutton, A.J., Sabine, C.L., Feely, R.A., Cai, W.J., Cronin, M.F., McPhaden, M.J., Morell, J.M., Newton, J.A., Noh, J.H., Ólafsdóttir, S.R., Salisbury, J.E., Send, U., Vandemark, D.C., Weller, R.A., 2016. Using present-day observations to detect when anthropogenic change forces surface ocean carbonate chemistry outside preindustrial bounds. Biogeosciences 13, 5065–5083. https://doi.org/10.5194/bg-13-5065-2016.

Taylor, J.R., 1996. An Introduction to Error Analysis: The Study of Uncertainties in Physical Measurements, 2nd ed. University Science Books, Mill Valley, CA.

Tellinghuisen, J., 2001. Statistical error propagation. J. Phys. Chem. A 105, 3917–3921.

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

van Heuven, S., Pierrot, D., Rae, J.W.B., Lewis, E., Wallace, D.W.R., 2011. MATLAB Program Developed for CO₂ System Calculations. ORNL/CDIAC-105b. Carbon Dioxide Inf. Anal. Cent. Oak Ridge Natl. Lab., Oak Ridge, Tenn.

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

Williams, N.L., Juranek, L.W., Feely, R.A., Johnson, K.S., Sarmiento, J.L., Talley, L.D., Dickson, A.G., Gray, A.R., Wanninkhof, R., Russell, J.L., Riser, S.C., Takeshita, Y., 2017. Calculating surface ocean pCO₂ from biogeochemical ARGO floats equipped with pH: an uncertainty analysis. Glob. Biogeochem. Cycles 31, 591–604. https://doi.org/10.1002/2016GB005541.

Xu, Y.Y., Pierrot, D., Cai, W.J., 2017. Ocean carbonate system computation for anoxic waters using an updated CO₂SYS program. Mar. Chem. 195, 90–93. https://doi.org/10.1016/j.marchem.2017.07.002.

Yang, B., Byrne, R.H., Lindemuth, M., 2015. Contributions of organic alkalinity to total alkalinity in coastal waters: a spectrophotometric approach. Mar. Chem. 176, 199–207. https://doi.org/10.1016/j.marchem.2015.09.008.

Yang, Y., Hansson, L., Gattuso, J.P., 2016. Data compilation on the biological response to ocean acidification: an update. Earth Sys. Sci. Data. 8, 79–87. https://doi.org/10.5194/essd-8-79-2016.

Yu, W., Blair, M., 2013. DNAD, a simple tool for automatic differentiation of Fortran codes using dual numbers. Comput. Phys. Commun. 184, 1446–1452.