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

HAL
open science

# Improved routines to model the ocean carbonate system: mocsy 2.0

J. C. Orr, J.-M. Epitalon

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

## To cite this version:

J. C. Orr, J.-M. Epitalon. Improved routines to model the ocean carbonate system: mocsy 2.0. Geoscientific Model Development, 2015, 8 (3), pp.485 - 499. 10.5194/gmd-8-485-2015. hal-01806242

HAL Id: hal-01806242

https://hal.science/hal-01806242

Submitted on 5 May 2021

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_3.png" /></figure>

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

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

# Improved routines to model the ocean carbonate system: mocsy 2.0

J. C. Orr¹ and J.-M. Epitalon¹, *

¹LSCE/IPSL, Laboratoire des Sciences du Climat et de l'Environnement, CEA-CNRS-UVSQ, Gif-sur-Yvette, France

* now at: Geoscientific Programming Services, Fanjeaux, France

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

Received: 24 March 2014 – Published in Geosci. Model Dev. Discuss.: 5 May 2014

Revised: 18 December 2014 – Accepted: 19 December 2014 – Published: 9 March 2015

**Abstract.** Modelers compute ocean carbonate chemistry often based on code from the Ocean Carbon Cycle Model Intercomparison Project (OCMIP), last revised in 2005. Here we offer improved publicly available Fortran 95 routines to model the ocean carbonate system (mocsy 2.0). Both codes take as input dissolved inorganic carbon $C_T$ and total alkalinity $A_T$, tracers that are conservative with respect to mixing and changes in temperature and salinity. Both use the same thermodynamic equilibria to compute surface-ocean $pCO_2$ and simulate air-sea $CO_2$ fluxes, but mocsy 2.0 uses a faster and safer algorithm (SolveSAPHE) to solve the alkalinity-pH equation, applicable even under extreme conditions. The OCMIP code computes only surface $pCO_2$, while mocsy computes all other carbonate system variables throughout the water column. It also avoids three common model approximations: that density is constant, that modeled potential temperature is equal to in situ temperature, and that depth is equal to pressure. Errors from these approximations grow with depth, e.g., reaching 3% or more for $pCO_2$, $H^+$, and $\Omega_A$ at 5000 m. The mocsy package uses the equilibrium constants recommended for best practices. It also offers two new options: (1) a recently reassessed total boron concentration $B_T$ that is 4% larger and (2) new $K_1$ and $K_2$ formulations designed to include low-salinity waters. Although these options enhance surface $pCO_2$ by up to 7 µatm, individually, they should be avoided until (1) best-practice equations for $K_1$ and $K_2$ are reevaluated with the new $B_T$ and (2) formulations of $K_1$ and $K_2$ for low salinities are adjusted to be consistent among pH scales. The common modeling practice of neglecting alkalinity contributions from inorganic P and Si leads to substantial biases that could easily be avoided. With standard options for best practices, mocsy agrees with results from the CO2SYS package within 0.005% for the three inorganic carbon species (concentrations differ by less than

0.01 µmol kg⁻¹). Yet by default, mocsy's deep-water $fCO_2$ and $pCO_2$ are many times larger than those from older packages, because they include pressure corrections for $K_0$ and the fugacity coefficient.

## 1 Introduction

To compute air-sea $CO_2$ fluxes, ocean carbon cycle models compute the partial pressure of carbon dioxide $pCO_2$ from two passive tracers, namely dissolved inorganic carbon $C_T$ and total alkalinity $A_T$. In many models, that thermodynamic calculation is based on documented code from the Ocean Carbon Cycle Model Intercomparison Project (OCMIP), which is publicly available at http://ocmip5.ipsl.jussieu.fr/OCMIP/. Although modified versions of that code are used widely (e.g., Müller et al., 2008; Doney et al., 2006; Aumont and Bopp, 2006), it has not been updated since 2005. Meanwhile, there have been developments in recommended community standards for equilibrium constants (Dickson et al., 2007; Dickson, 2010).

Models require computationally efficient routines that are compatible with other model components, typically written in Fortran. Hence, model simulations are not made with other public software packages used widely to compute remaining carbonate system variables from any input pair, given corresponding temperature $T$, salinity $S$, and pressure $P$, as well as concentrations of total dissolved inorganic phosphorus $P_T$ and total dissolved inorganic silicon $Si_T$ (Dickson et al., 2007). Models also differ because they typically carry potential temperature $\theta$, use concentration units of moles per cubic meter ($mol m^{-3}$), and are referenced to depth (m); conversely, equations for carbonate system thermodynamics require in situ $T$, concentrations in moles per kilogram
(mol kg⁻¹), and in situ P. Unit conversion is straightforward, but for simplicity modelers often make three approximations: (1) that θ is equivalent to T, (2) that ocean density is constant (e.g., 1028 kg m⁻³), and (3) that depth (m) is equivalent to in situ pressure (dbar).

Errors associated with these simplifications are considered to be negligible, while the reasons behind them are largely historical. Most studies with ocean carbon cycle models have focused on large-scale patterns of air-sea CO₂ fluxes and related near-surface changes in the open ocean. More recently, with growing concern for ocean acidification (Caldeira and Wickett, 2003; Orr et al., 2005), attention has also turned to the deep ocean, the high latitudes, and local impacts in the coastal zone. Errors associated with the three approximations will be larger at depth, where θ diverges from T, where densities are greater than average surface values, and where there are larger absolute differences between pressure and depth. One may also question use of the constant density approximation in waters affected by excess evaporation (e.g., in the equatorial Pacific, Arabian Sea, and Mediterranean Sea) or by large freshwater input (e.g., in the Arctic and coastal zones with heavy river influence).

To fill these gaps, we provide here an improved set of routines to model the ocean carbonate system (mocsy 2.0). This new package uses the classic approach, taking simulated $A_T$ and $C_T$ and computing all other carbonate system variables, while adding refinements. Relative to its predecessor (OCMIP model code), mocsy offers several ameliorations: (1) it no longer makes the three approximations mentioned above; (2) it computes all carbonate system variables, not only $pCO_2$ and pH; (3) it provides these variables at all model levels, not only at the surface; and (4) it uses, as a default, the constants and the pH scale recommended for best practices (Dickson et al., 2007). All constants except $K_0$ (the $CO_2$ solubility) are corrected for pressure effects (Millero, 1995) with modified coefficients adopted from CO2SYS (Lewis and Wallace, 1998); both $K_0$ and the fugacity coefficient are corrected for pressure effects following Weiss (1974). Options are also provided to replace formulations for $K_1$ and $K_2$ and for total boron that are recommended for best practices (Dickson et al., 2007) with more recent formulations, choices that will be shown to affect results significantly.

# 2 Methods

## 2.1 Code description

### 2.1.1 Basics

The ocean carbonate system is well constrained. Any pair of carbonate system variables can be used to compute all others. Because only two carbonate system variables, $C_T$ and $A_T$, are carried as passive tracers by all ocean carbon models, mocsy

offers only that input pair. Thus it is unlike other public packages designed largely for observationalists. To calculate the remaining carbonate system variables, one only needs to provide additional input for temperature, salinity, $P_T$, and $Si_T$ as well as pressure or depth. A precursor to mocsy was developed in 2004 and used to project future ocean acidification from simulated $C_T$ and $A_T$ in the OCMIP2 models (Orr et al., 2005). That precursor code was never released publicly and should not be confused with the preexisting OCMIP2 or OCMIP3 code, which only computes surface $pCO_2$. Its development began by combining the Fortran code for equilibrium constants from OCMIP2 and an iterative algorithm to solve for pH (Maier-Reimer, 1993; Aumont and Bopp, 2006). The precursor code was then modified to compute all carbonate system variables throughout the water column. Thus it included pressure corrections for equilibrium constants (Millero, 1995) with pressure-correction coefficients (some now known to be erroneous) taken from version 0.95 of seacarb (Proye and Gattuso, 2003; Lavigne and Gattuso, 2011), which itself adopted code from csys (Zeebe and Wolf-Gladrow, 2001). From the seacarb code, the precursor code also included the analytical formula for the Revelle factor from Frankignoulle (1994). Our feedback from this early development led to bug corrections that were later implemented in seacarb.

Since 2005, this precursor code has continued to be improved. The first public version was denoted as mocsy 1.0 (Orr and Epitalon, 2014). Here along with a revised version of that discussion paper, we provide an improved version of the code, mocsy 2.0. The equilibrium constants and pH scale adopted in mocsy are those recommended by the Guide to Best Practices for Ocean $CO_2$ Measurements (Dickson et al., 2007). The equilibrium constants are generally on the total pH scale. There are four exceptions: the $CO_2$ solubility $K_0$, the equilibrium constant for bisulfate $K_S$, the solubility product for aragonite $K_A$, and the analogous solubility product for calcite $K_C$ (Mucci, 1983). By definition, $K_0$, $K_A$, and $K_C$ are independent of the pH scale: they do not involve $[H^+]$. Although $K_S$ includes $[H^+]$, it is maintained on the free scale, as needed to convert between the free and total scales. Conversely, the equilibrium constant that is used to convert between total and seawater scales, i.e., $K_F$ for the dissociation of hydrogen fluoride HF, is maintained on the total scale along with all other constants.

For its basic calculations, mocsy adopts the recommendations of Dickson et al. (2007). These include (1) the Weiss (1974) formulation describing the solubility of $CO_2$ in seawater ($K_0$); (2) the Lueker et al. (2000) formulations for first and second dissociation constants of carbonic acid ($K_1$ and $K_2$), refits of measurements from Mehrbach et al. (1973) on the NBS scale to the total pH scale; (3) the Millero (1995) formulations for equilibrium constants of boric acid ($K_B$), phosphoric acid ($K_{1P}$, $K_{2P}$, $K_{3P}$), silicic acid ($K_{Si}$), and water ($K_W$), which mocsy converts from the seawater scale to the total scale; (4) the Dickson (1990) formulation for the
equilibrium constant for the dissociation of bisulfate ($K_S$) on the free scale (see above); (5) the Perez and Fraga (1987) formulation for the equilibrium constant for hydrogen fluoride HF ($K_F$) on the total scale; and (6) the Mucci (1983) formulations for the $CaCO_3$ solubility products for aragonite and calcite ($K_A$ and $K_C$). These equilibrium constants use concentrations, not activities. All of them except for $K_0$ are further adjusted for pressure using the approach of Millero (1995), with corrected coefficients from Lewis and Wallace (1998) (see Orr et al., 2015, Table 7). Constant ratios relative to salinity are used to compute concentrations of total inorganic boron (Uppström, 1974), fluoride (Riley, 1965), sulfur (Morris and Riley, 1966), and calcium (Riley and Tongudai, 1967). The product of the $Ca^{2+}$ and $CO_3^{2-}$ concentrations divided by the apparent solubility product (either $K_A$ or $K_C$) yields the saturation state (i.e., for aragonite $\Omega_A$ or for calcite $\Omega_C$, respectively).

### 2.1.2 Solver of the total alkalinity-pH equation

In mocsy 2.0, to solve the total alkalinity-pH equation, we have replaced the classic fixed-point iterative carbonate alkalinity scheme (ICAC) in mocsy 1.0 with a new, universally convergent algorithm from Munhoven (2013). We now call a cubic initialization routine followed by the standard safeguarded solver routine `solve_at_general`, both from Munhoven's SolveSAPHE package v1.0.1. Our tests using typical open-ocean conditions indicate that results from both are identical to at least the sixth digit after the decimal in terms of pH, but that SolveSAPHE is about 5 times faster than our former ICAC routine.

Moreover, the SolveSAPHE routines allow mocsy 2.0 to avoid any convergence problem under extreme conditions. The ICAC methods do not always converge, for example under low $C_T$ when $A_T/C_T > 1$ (Munhoven, 2013, Fig. 3c). Our experience is that non-convergence may occur in high-resolution, global biogeochemical models under present-day conditions, because local effects from freshwater river input are heightened. An example comes from recent test simulations with the ICAC scheme embedded in 2°, 0.5°, and 0.25° versions of the NEMO-PISCES model at our laboratory. At the highest resolution, the ICAC scheme generates negative $[H^+]$ in the model's Ocean Eustasy, where simulated $A_T$ and $C_T$ are both very low, as are salinities. We expect that these negative $[H^+]$ will disappear once the SolveSAPHE algorithms are implemented in the model.

For use with mocsy 2.0, we made minor modification to the SolveSAPHE routines (1) to use mocsy's equilibrium constants and (2) to remove arguments and equations for $NH_4^+$ and $H_2S$ acid systems.

### 2.1.3 From aqueous CO₂ to pCO₂

In mocsy, we use $K_0$ from Weiss (1974) to compute $fCO_2$ from $CO_2^*$ and we use the fugacity coefficient $C_f$ (Weiss,

1974; Dickson and Goyet, 1994; Dickson et al., 2007) to compute $pCO_2$ from $fCO_2$:

$$[CO_2^*] = K_0 fCO_2 = K_0 C_f pCO_2. \quad (1)$$

For this calculation of oceanic $pCO_2$, mocsy does not use the combined coefficient $F$ from Weiss and Price (1980), where

$$F = K_0 C_f (1 - pH_2O), \quad (2)$$

because it includes a wet-to-dry-air conversion (term in parentheses), which is inappropriate for the conversions in Eq. (1). Rather, that combined term is commonly used when converting between atmospheric $pCO_2$ and $xCO_2$. Although those atmospheric conversions are now also offered with new functions in mocsy 2.0, we still do not use $F$ in order to make calculations in a stepwise fashion and to avoid potential confusion.

The mocsy 2.0 package also differs substantially from other packages because it adjusts for effects of subsurface pressure on $K_0$ as well as $C_f$. All packages compared by Orr et al. (2014) compute $K_0$ with the same standard equation (Weiss, 1974, Eq. 12), but none of them make the exponential pressure correction (Weiss, 1974, Eq. 5):

$$[CO_2^*] = K_0 fCO_2 \exp[(1-P)\bar{v}_{CO_2/RT}], \quad (3)$$

where $P$ is the total pressure (atmospheric + hydrostatic, both in units of atm), $\bar{v}_{CO_2}$ is the partial molal volume of $CO_2$ (32.3 $cm^3mol^{-1}$), $R$ is the gas constant (82.05736 $cm^3atm mol^{-1}K^{-1}$), and $T$ is the absolute temperature (K). Thus the computed $fCO_2$ refers only to potential values considering the total pressure as that at the surface. The exponential pressure-correction term can be neglected but only when the total pressure is near 1 atm (Weiss, 1974). It is accounted for in mocsy 2.0.

Likewise the effects of subsurface pressure on the fugacity coefficient $C_f$ were not considered in the public packages compared by Orr et al. (2014), including mocsy 1.0 (Orr and Epitalon, 2014). Those packages considered that total pressure was the same as atmospheric pressure, always equal to 1 atm. Conversely, mocsy 2.0 also accounts for subsurface pressure effects on $C_f$ (Weiss, 1974, Eq. 9):

$$fCO_2 = pCO_2C_f = pCO_2 \exp \left[ \left( B + 2x_2^2 \delta_{12} \right) P/RT \right], \quad (4)$$

where $B$ is the virial coefficient of $CO_2$ (Weiss, 1974, Eq. 6), $x_2$ is the sum of the mole fractions of all remaining gases (1 - $xCO_2$, when $xCO_2 \ll 1$), and $\delta_{12} = 57.7 - 0.118T$. Once again, $P$ is the total pressure (atmospheric + hydrostatic) in atmospheres (atm), and $R$ and $T$ are as described for Eq. (3).

In Eqs. (3) and (4), other public packages have commonly used atmospheric pressure for $P$, typically fixed at 1 atm. Their computed subsurface $fCO_2$ and $pCO_2$ may be considered as pseudo-potential values that a water parcel would
<figure><img src="image_6.png" /><figcaption>Figure 1. Effect of pressure on K₀ and C_f given as factors that would need to be multiplied by (left) fCO₂ and (right) pCO₂ computed without pressure corrections. The reference lines (black dashed) are the pseudo-potential values computed from surface P and in situ T, the potential values (red) are computed from surface P and potential T, and the in situ values (blue) are computed from total P and in situ T. These global-mean vertical profiles were computed with mocsy 2.0 and the GLODAP-WOA2009 gridded input data.</figcaption></figure>

have if brought back to the surface. By pseudo, we mean that this common approach does not correct for the effect of pressure on temperature. Hence we offer a third approach to compute the true potential fCO₂ and pCO₂ also at surface pressure but using the potential temperature θ in place of T in Eqs. (3) and (4). That is, these true potential fCO₂ and pCO₂ also consider that subsurface water is brought adiabatically back to the surface. In mocsy 2.0, all three approaches are available, but the default is to compute fCO₂ and pCO₂ with the first approach, at in situ temperature and total pressure.

The effects of these pressure corrections on K₀ and C_f are shown in terms of their multiplicative effect on fCO₂ and pCO₂. For K₀, the effect is shown as components of the exponential term in Eq. (3) (Fig. 1, left). For C_f, it is shown as the ratio between the exponential term in Eq. (4), where P is total pressure over the same term with P equal to atmospheric pressure (Fig. 1, right). At 100 m, factors already reach 1.014 for fCO₂ and 1.05 for pCO₂. At 5000 m, the total multiplicative effect reaches more than a factor of 2 for fCO₂ and 18 for pCO₂. The effect of pressure on temperature (T) is a small effect compared with the total correction (T + P). Hence the second and third approaches above provide similar results. Because mocsy 2.0's default is to include these large subsurface pressure corrections, its subsurface fCO₂ and pCO₂ will be much larger than those from other packages that make no such corrections. Conversely, calculations of other variables will not be affected by pressure adjustments of K₀ and C_f. The same pressure corrections were recently included in a second package, seacarb 3.0.6, and its results compare well with mocsy 2.0 (Orr et al., 2015, Tables 9 and 10).

### 2.1.4 Revelle factor

Another new feature in mocsy 2.0 is that it accounts for the effects of P_T and Si_T on the Revelle factor R_f.

$$R_f = \frac{\partial pCO_2 / pCO_2}{\partial C_T / C_T} \qquad (5)$$

$$= \frac{\partial pCO_2}{\partial C_T} \left( \frac{C_T}{pCO_2} \right). \qquad (6)$$

To do so, mocsy 2.0 replaces the analytical formula for R_f from Frankignoule (1994) that was used in mocsy 1.0 with a numerical centered-finite-difference solution to compute the derivative. This second-order accurate numerical approach is like that used in CO2SYS, with two differences. The most significant of those is that CO2SYS uses fCO₂ in place of pCO₂. Secondly, in CO2SYS the step size h for the numerical approximation of the derivative (df/dx ≈ [f(x₀ + h) - f(x₀ - h)]/2h, where x is C_T and f(x) is pCO₂) is 10 times larger than ours (h = 0.1 µmol kg⁻¹), an optimal value that produced the minimum difference between the numerical and analytical solutions, when nutrient concentrations are set to 0. However, both step sizes yield identical computed R_f values up to the fourth decimal place.

## 2.2 Options

Since the publication of the last best-practices guide (Dickson et al., 2007), there have been developments that merit close attention given their potential impacts on computed carbonate chemistry variables. First, Lee et al. (2010) estimate that the total boron concentration in the ocean, i.e., its linear relationship with salinity, is 4 % greater than estimated previously (Uppström, 1974). Lee et al. (2010) used a more precise measurement technique on more samples (n = 139) collected from a wider geographic distribution than did Uppström (1974), whose total boron: salinity ratio is based on 20 samples from the deep Pacific.

Second for K₁ and K₂, Millero (2010) combined the same set of measurements used by Lueker et al. (2000) from Mehrbach et al. (1973) along with his own (Millero et al., 2006) to fit new formulations that are applicable over larger ranges of salinity (1 to 50) and temperature (0 to 50 °C) than are those recommended for best practices (Lueker et al., 2000). The latter are intended to be used only when 19 < S < 43 and 2 < T < 35 °C. Salinities and temperatures below these thresholds do occur even in coarse-resolution global models in areas such as the Arctic, which routinely experiences subzero temperatures and intense freshwater input from rivers as well as land- and sea-ice melt. Low salinities near rivers are also common in regional models and will become more prevalent in global models as resolution increases. To model such conditions properly, we may need to go beyond the best-practices recommendation (Dickson et al., 2007). Thus in mocsy, we have implemented options
<figure><img src="image_7.png" /><figcaption>**Figure 2.** Surface zonal means of carbonate system variables computed with mocsy 2.0 from the gridded GLODAP data for $A_T$ and $C_T$ combined with corresponding data from the 2009 World Ocean Atlas for in situ T, S, $P_T$, and $Si_T$.</figcaption></figure>

for the user to choose to use the newer formulations mentioned above for $K_1$, $K_2$, and total boron in place of those recommended for best practices. Likewise for $K_F$, we allow the user to choose either the Dickson and Riley (1979) formulation recommended by Dickson and Goyet (1994) or the Perez and Fraga (1987) formulation recommended by Dickson et al. (2007), but which is intended to be limited to waters where $10 < S < 40$ (practical salinity scale) and $9 < T < 33$ °C.

## 2.3 Exceptions to best practices

The best-practice formulations for $K_W$, $K_{1P}$, $K_{2P}$, and $K_{3P}$ proposed by Dickson et al. (2007) are those from Millero (1995) with 0.015 subtracted from the constant term as a simple means to convert from the seawater to the total hydrogen ion scale (Dickson and Goyet, 1994; Dickson et al., 2007, Chap. 5, footnote 5). In mocsy, we do not impose this constant correction, preferring instead to use the classic approach to convert equilibrium constants between the two pH scales (e.g., Millero, 2010, Eq. 6). That results in a pH-scale correction that varies with [HF]. The variable correction is also used in CO2SYS.

Dickson et al. (2007) do not discuss pressure corrections of equilibrium constants. Unlike other public packages, mocsy 2.0 makes pressure adjustments for $K_0$ and the related $C_f$ (Eqs. 3 and 4), thus affecting computed $f\,CO_2$ and $pCO_2$ (see Sect. 2.1.3). For the remaining equilibrium constants, mocsy follows the lead of CO2SYS using Millero's equations, quadratic functions of pressure and temperature

(Millero, 1995, Eqs. 90–92) with corrections to associated coefficients from Lewis and Wallace (1998) as detailed elsewhere (Orr et al., 2015, Table 7).

To compute the $Ca^{2+}$ concentration from salinity, mocsy 2.0 uses a Ca-to-chlorinity ratio of 0.02128 from Riley and Tongudai (1967) instead of the slightly different value of 0.02127 that appears in the best-practices guide (Dickson et al., 2007, Chapt. 5, Table 2). The latter states it uses the ratio from the former reference.

## 2.4 Evaluation

There is no absolute reference for computed carbonate system variables. To validate mocsy, its computed variables were compared to those from CO2SYS-MATLAB (van Heuven et al., 2011) run with identical input data. Although both packages could in principle be wrong even if both agree, we compare mocsy to CO2SYS because the latter is the first public package made available to compute these variables, it is used widely, and it was developed with great care (Lewis and Wallace, 1998). Moreover, it has served as a base for other public packages to build on. A companion paper (Orr et al., 2015) further details the reasons for our arbitrary choice of CO2SYS as the reference. The input data with which we compared the two packages include $A_T$ and $C_T$ from the three-dimensional, global gridded data product on a 1° × 1° grid known as GLODAP (Key et al., 2004). Other necessary input data on the same grid were taken from the 2009 World Ocean Atlas (WOA2009) gridded data product, which includes in situ temperature (Locarnini et al., 2010),
<figure><img src="image_8.png" /><figcaption>**Figure 3.** Global-mean vertical profiles of carbonate system variables computed with mocsy 2.0 from the data described in Fig. 2.</figcaption></figure>

salinity (Antonov et al., 2010), and concentrations of $P_T$ and $Si_T$ (Garcia et al., 2010).

## 3 Results

As a baseline reference for later evaluation and sensitivity tests, we present carbonate system variables computed with mocsy 2.0 from the gridded GLODAP-WOA2009 data. Magnitudes and patterns of the surface zonal-mean distributions (Fig. 2) are as expected based on previous studies (Orr et al., 2005; Orr, 2011). Corresponding global-mean vertical profiles (Fig. 3) have not been detailed previously. Both $f\text{CO}_2$ and $p\text{CO}_2$ continue to increase with depth, particularly the latter, because of pressure corrections to $K_0$ and $C_f$, respectively. Conversely, the $\text{CO}_2^*$ concentration increases from the surface to about 1000 m, where it reaches 3 times the surface level. Below, it declines slowly to values at 5000 m that are just over twice those at the surface. The $\text{CO}_3^{2-}$ profile mirrors that for $\text{CO}_2^*$, as expected from the zonal-mean distributions. Vertical distributions of the saturation states $\Omega_A$ and $\Omega_C$ generally follow the $\text{CO}_3^{2-}$ profile but are influenced by pressure effects on $K_A$ and $K_C$. The shape of the $H^+$ concentration profile is similar to that for $\text{CO}_2^*$ but shows two maxima at 1000 and 5000 m, both at about twice the surface concentration. The Revelle factor also reaches a maximum at 1000 m, but it declines only slightly below, with values remaining at nearly double the surface level.

### 3.1 Evaluation

Our evaluation reveals the extent to which mocsy's variables computed from the GLODAP-WOA2009 data dif-

fer from those computed when the same data is used with the CO2SYS-MATLAB package. Relative differences remain within ±0.005% for computed concentrations of $\text{CO}_2^*, \text{HCO}_3^-, \text{CO}_3^{2-}$, and $H^+$, both for surface zonal means (Fig. 4) and global average vertical profiles (Fig. 5). Corresponding absolute differences are within 0.01 µmol kg⁻¹ for each of the three inorganic carbon species, and within 0.0002 nmol kg⁻¹ for $H^+$. The latter translates into pH differences of less than 0.00002 between the two packages. Tighter agreement cannot be expected because CO2SYS computes pH with a Newton-Raphson method that stops iterating once the change in pH from the previous iteration is less than 0.0001. The slight difference in $H^+$ between packages explains the slight offsets in the three inorganic carbon species, because both packages yield identical results for the computed $K_1$ and $K_2$ equilibrium constants (Orr et al., 2015).

Relative differences between packages are slightly larger for $\Omega_A$ and $\Omega_C$, with a constant offset of 0.018% throughout the water column. Comparison of the source code of the two packages revealed that the offset is due to differences in the atomic weight of calcium: CO2SYS uses 40.087, whereas mocsy uses 40.078, the recommended value (Dickson et al., 2007, Chapt. 5, Table 1). Nonetheless, these relative differences remain small. They correspond with differences in computed saturation states of less than 0.001 unit between the two packages. Of similar magnitude are the relative differences for the Revelle factor $R_f$ of about 0.02% at the surface. Relative differences in $R_f$ grow to 0.2% at 5000 m, but the corresponding absolute differences remain less than 0.03. With the new numerical solution of $R_f$ in mocsy, differences relative to CO2SYS are 6 times smaller than with the analytical formula used in mocsy 1.0.
<figure><img src="image_9.png" /></figure>

**Figure 4.** Differences in surface zonal means between mocsy and CO2SYS for computed variables shown as (top) percent relative differences ($100(V_{\text{mocsy}} - V_{\text{co2sys}})/V_{\text{co2sys}}$) and as (bottom) absolute differences ($V_{\text{mocsy}} - V_{\text{co2sys}}$). The absolute differences are given in µmol kg⁻¹ for the three inorganic carbon species, in nmol kg⁻¹ for H⁺, and in µatm for fCO₂ and pCO₂; other computed variables are unitless.

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

**Figure 5.** Differences in global-mean profiles between mocsy and CO2SYS for computed variables, as in Fig. 4. The large differences in fCO₂ and pCO₂ (top right) are shown as ratios (mocsy/CO2SYS).

The two remaining computed variables, surface fCO₂ and pCO₂, both track CO₂* in both packages. Hence differences between packages for these variables remain small (0.005 %) at the surface. Yet discrepancies grow to extraordinary levels at greater depths. At 5000 m, mocsy's fCO₂ is twice as large as that in CO2SYS, while its pCO₂ is 18 times greater. The cause of these large differences is the pressure corrections to

$K_0$ and $C_f$ in mocsy 2.0 (Eqs. 3 and 4) as shown in Fig. 1, which are not accounted for in CO2SYS.

## 3.2 Model approximations

Errors from all three model approximations increase with depth (Fig. 6). Maximum total errors reach 3 % or more for pCO₂, H⁺, and ΩA at 5000 m. Yet the causes differ. For
<figure><img src="image_11.png" /><figcaption>**Figure 6.** Relative changes in global-mean profiles that result from avoiding each of the three ocean-model approximations: (1) that density is constant (blue ρ), (2) that potential temperature is equivalent to in situ temperature (red T), and (3) that depth (m) is equivalent to pressure (dbar) (green Z). Also shown is the sum of all three effects (black). Associated errors are equal but opposite in sign to the changes shown.</figcaption></figure>

$p\text{CO}_2$, the largest error comes from the pressure approximation, while errors from the other two approximations are smaller and nearly compensate one another. For $\Omega_A$ there is little compensation, while the error from the constant density approximation dominates; it alone reaches 3% at 5000 m. For $H^+$, there is no compensation and errors from the temperature approximation dominate. The total relative error for computed $f\text{CO}_2$ in the deep ocean is about half that for $p\text{CO}_2$ because the error associated with the pressure approximation is more than 4 times smaller while errors from the other two approximations are similar. Total errors from the three approximations remain less than 1% even in the deep ocean for the three inorganic carbon species as well as for the Revelle factor.

Errors in $\Omega_A$ have somewhat similar patterns to that of $\text{CO}_3^{2-}$. Yet at 5000 m, the magnitude of the total relative error of $\Omega_A$ reaches 4%, more than 5 times greater than for $\text{CO}_3^{2-}$. That relative enhancement of the former over the latter is partly due to the sharp decline in $\Omega_A$, e.g., below 1000 m, whereas $\text{CO}_3^{2-}$ concentrations remain relatively stable (Fig. 3). The dominance of the pressure-induced increase in $K_A$ largely outweighs its sensitivity to temperature, unlike for $\text{CO}_3^{2-}$.

None of the three approximations is without error, but it is the constant density approximation which leads to the largest errors in most computed variables (except for $p\text{CO}_2$, $H^+$, and $f\text{CO}_2$). Yet total absolute errors generally remain small, even at 5000 m, e.g., less than 0.6 µmol kg⁻¹ for the three inorganic carbon species, less than 0.02 for $\Omega_A$, and less than

0.07 for $R_f$ (Fig. 7). Conversely, at the same depth the absolute error in pH reaches 0.015, while for $p\text{CO}_2$ it reaches 350 µatm.

### 3.3 Options

The effects of the three options on computed variables differ by region and depth. The choice for the $K_F$ option has virtually no effect on computed quantities and can be neglected. That is, $K_F$ appears in the total alkalinity equation as detailed later, but its relative importance may be thought of simply in terms of pH, as it is used to convert between the seawater and total pH scales. The difference between the two scales is small (~0.01 units of pH), and the two formulations alter that very slightly. Practically then, computed carbonate system variables are insensitive to different formulations of $K_F$. More details about the sensitivity of computed variables to each of the constants can be found elsewhere (Orr et al., 2015, Table 11). Larger differences result from the choice of the two other new options, for total boron and for formulations for $K_1$ and $K_2$.

The 4% increase in total boron with the new formulation from Lee et al. (2010) increases borate alkalinity $A_B$, e.g., by about 3% in typical surface waters. The latter increases less because $A_B = B_T/(1+[H^+]/K_B)$ and $[H^+]$ also increases by 1%. Yet modeled total alkalinity $A_T$ is unaffected, being an input variable (along with $C_T$). Because $A_T$ is unchanged and $A_B$ is higher, modeled carbonate alkalinity $A_C$ must be lower. This decline in $A_C$ reduces sur-
<figure><img src="image_12.png" /><figcaption>**Figure 7.** Absolute changes in global-mean profiles of computed variables when each of the three ocean-model approximations is avoided, with line colors and patterns as in Fig. 6.</figcaption></figure>

face CO₃²⁻ everywhere, from -1.5 µmol kg⁻¹ in the Southern Ocean to -3 µmol kg⁻¹ in the tropics (Fig. 8). The corresponding decline in surface ΩA is between 0.02 and 0.04. Simultaneously, computed surface pCO₂ increases by 4 to 6 µatm and surface pH declines by 0.006 units. Yet the Revelle factor declines by 0.04 (the buffer capacity increases). That decline can be understood by studying Eq. (5). The partial derivative on the right-hand side increases by roughly 1% when the new boron formulation (Lee et al., 2010) is used in place of the standard (Uppström, 1974). However the adjacent concentration ratio (in parentheses) decreases by relatively more, about 1.4%. Hence Rf decreases because pCO₂ increases (Cₜ remains constant). In the deep ocean, the new total boron formulation leads to nearly uniform changes in pH of -0.005 ± 0.001 units, and in CO₃²⁻ of -1 ± 0.1 µmol kg⁻¹. There is also an exponential increase in pCO₂ up to 110 µatm at 5000 m (Fig. 9). The latter represents a 1% increase in the baseline pCO₂ (Fig. 3).

The new option where the formulations of K₁ and K₂ from Millero (2010) are used instead of those from Lueker et al. (2000) yields changes that are less uniform and in some places larger than those from using the new option for Bₜ. Changes from the new option for K₁ and K₂ are larger in the Southern Ocean: pCO₂ increases by up to +7 µatm, corresponding to a 2% increase in CO₂*, while pH declines by up to 0.006. In contrast, in the lowest latitudes corresponding changes in pCO₂ and CO₂* are negligible. In the Southern Ocean, changes in pCO₂, CO₂*, and pH from using the new K₁ and K₂ option reinforce those from the new Bₜ option. But in the same region, changes from the two options

partly compensate one another for CO₃²⁻ and ΩA, while they reinforce one another in the lower latitudes. Due to these regional differences in compensation, total absolute changes are largest in the Southern Ocean for pCO₂ (+12 µatm), CO₂* (+4%), pH (-0.012), and Rf (-0.30), whereas they are largest in the lower latitudes for CO₃²⁻ (-3.5 µmol kg⁻¹) and ΩA (-0.06). With depth, mean absolute changes due to both options grow for pCO₂ and CO₂*, change little for pH and Rf, and decline for CO₃²⁻ and ΩA. These differences appear greatly affected by the structure of the baseline profiles (Fig. 3).

## 4 Discussion

Our focus has been on quantifying errors in computed variables from the three model approximations and in assessing how variables are affected by the three user options. A more general concern is how computed variables are affected by the frequent practice of neglecting nutrient concentrations in carbonate system calculations. In high-nutrient regions, the changes in computed variables due to alkalinity from Pₜ and Siₜ (Figs. 10 and 11) are similar in magnitude and have the same sign as do those from changing to the new formulation for total boron (Figs. 8 and 9). All three contribute to non-carbonate alkalinity and hence total alkalinity,

$$A_T = A_C + A_B + A_W + A_P + A_{Si} + A_O, \quad (7)$$
<figure><img src="image_13.png" /><figcaption>**Figure 8.** Changes in zonal means of computed variables due to replacing best-practice recommendations with other options: (1) the Lee et al. (2010) formulation for the total boron-to-salinity ratio (blue B), (2) the Millero (2010) formulations for K₁ and K₂ (red, K₁K₂), and (3) the Dickson and Riley (1979) formulation for K<sub>F</sub> (green K<sub>F</sub>). Absolute differences are shown relative to (1) Uppström (1974), (2) Lueker et al. (2000), and (3) Perez and Fraga (1987). Also indicated is the sum of the first two changes (black).</figcaption></figure>

where

$$A_C = [\text{HCO}_3^-] + 2[\text{CO}_3^{2-}], \qquad (8)$$

$$A_B = [\text{B(OH)}_4^-], \qquad (9)$$

$$A_W = [\text{OH}^-] - [\text{H}^+]_F - [\text{HSO}_4^-] - [\text{HF}], \qquad (10)$$

$$A_P = [\text{HPO}_4^{2-}] + 2[\text{PO}_4^{3-}] - [\text{H}_3\text{PO}_4], \qquad (11)$$

$$A_{\text{Si}} = [\text{SiO(OH)}_3^-], \qquad (12)$$

$$A_O = [\text{NH}_3] + [\text{HS}^-] + \dots \qquad (13)$$

In brief, contributions on the right side of Eq. (7) come from components of carbonic acid, boric acid, water, phosphoric acid, silicic acid, and other species, respectively. The sum of last three terms in Eq. (10) is often represented simply as [H⁺], namely the hydrogen ion concentration on the seawater scale. Hence, water alkalinity A<sub>W</sub> provides an indicator of how a given seawater sample differs from acid-base neutrality.

As an input variable, A<sub>T</sub> is not affected by using a more simplified alkalinity equation, by neglecting nutrient concentrations, or by choosing a different formulation for total boron. Indeed, measured A<sub>T</sub> is unaffected by our calculations. Yet A<sub>C</sub> is affected, being computed by difference. As A<sub>B</sub>, A<sub>P</sub>, or A<sub>Si</sub> increases, computed A<sub>C</sub> must decrease. In the surface ocean, nutrient alkalinity (A<sub>P</sub> + A<sub>Si</sub>) substantially alters computed carbonate system variables where nutrient concentrations are largest, e.g., in the equatorial Pacific and

in the high latitudes (poleward of 40°). The largest surface effects occur in the Southern Ocean, where computed pCO₂ changes by +6 µatm, once nutrients are accounted for, which is 6 times more than estimated by Follows et al. (2006). At the same time, CO₃²⁻ changes by -1.6 µmol kg⁻¹ and pH by -0.007 (Fig. 10). Differences between some models may be larger in some regions, but that is far from a general rule. In the Southern Ocean, simulated air-sea fluxes of natural CO₂ differ between some models by only 0.1 Petagram of carbon per year (Pg C yr⁻¹), equivalent to less than 1 µatm in the air-sea difference in pCO₂ (Dufour et al., 2013, Table 1). The shift mentioned above would be more than enough to switch some of them from net sinks to net sources of natural CO₂.

These biases concern models that neglect contributions of phosphoric and silicic acids to alkalinity (by assuming A<sub>T</sub> = A<sub>C</sub> + A<sub>B</sub> + A<sub>W</sub>). Equivalent biases occur when nutrient concentrations are assumed to be 0 in offline calculations with output from models that include these nutrients in the alkalinity equation. Without other compensating biases, this simplification would lead to a simulated aragonite saturation horizon (ASH) that is too deep (Fig. 12), a simulated onset of aragonite undersaturation in polar surface waters (Ω<sub>A</sub> < 1) that is too late, and a simulated interhemispheric north-to-south oceanic transport of carbon (Sarmiento et al., 2000) that is too weak. Our results illustrate where and by how much models err when they neglect nutrient alkalinity, either because of a simplified alkalinity equation or equivalently
<figure><img src="image_14.png" /><figcaption>**Figure 9.** Absolute changes in global-mean vertical profiles of computed variables due to replacing best-practice recommendations with the three options mentioned in Fig. 8.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>**Figure 10.** Absolute changes in zonal means of computed variables due to alkalinity contributions from $P_T$ (blue solid), $Si_T$ (red dashed), and their sum (black dotted) relative to the case where alkalinity from phosphoric and silicic acids is neglected.</figcaption></figure>

by assuming null nutrient concentrations (e.g., in offline calculations). Remedying these errors requires little extra coding and does not add significantly to a model’s computation time. For models that do not carry $P_T$ and $Si_T$ as tracers, the bias in computed carbonate system variables would be reduced by calculating the alkalinity from those absent tracers as that from observed nutrient climatologies (e.g., Garcia

et al., 2010). Alternatively, models that carry inorganic nitrate but not inorganic phosphate as a tracer could approximate the latter by multiplying the former by a fixed Redfield ratio. Yet alkalinity associated with $Si_T$ would still be missing.

The combined effect from nutrient alkalinity and the new $B_T$ formulation are often greater than the individual effects
<figure><img src="image_16.png" /><figcaption>**Figure 11.** Absolute changes in global-mean profiles of computed variables attributable to nutrient alkalinity, i.e., for the same components shown in Fig. 10.</figcaption></figure>

in the surface ocean (Figs. 8 and 10). In the deep ocean, the average combined effect shifts CO₃²⁻ by -3 µmol kg⁻¹, pH by -0.018, and pCO₂ by 18 µatm (Fig. 9 + Fig. 11). The combined effect also shallows the computed modern ASH by about 100 m in the North Atlantic and up to 300 m in the Southern Ocean (Fig. 12).

Unfortunately, neither the new option for B<sub>T</sub> nor that for K<sub>1</sub> and K<sub>2</sub> can be recommended at present. The new boron option leads to substantial changes in computed variables, but those results rely on best-practice formulations of K<sub>1</sub> and K<sub>2</sub> that predate the Lee et al. (2010) study and depend them- selves on B<sub>T</sub> (Mehrbach et al., 1973, Eq. 8). Likewise, the option to use new formulations of K<sub>1</sub> and K<sub>2</sub> from Millero (2010) leads to large differences in computed variables, but those constants are not consistent when calculated from sets of coefficients made available on different pH scales (Orr et al., 2015).

## 5 Conclusions

Although fCO₂ and pCO₂ are typically measured or calcu- lated at the surface, e.g., to estimate air–sea CO₂ fluxes, a limited number of studies also focus on those quantities be- low the surface (e.g., Brewer and Peltzer, 2009; Cocco et al., 2013; Bates et al., 2013). Here we recommend that future studies concerned with subsurface values should choose one of two options: (1) compute in situ fCO₂ and pCO₂ af- ter making pressure corrections to K₀ and the fugacity co- efficient following Weiss (1974) using total in situ pressure and in situ temperature or (2) calculate potential fCO₂ and pCO₂ at 1 atm pressure while using potential instead of in

situ temperature. Results should then be clearly labeled as ei- ther in situ or potential fCO₂ and pCO₂. We make this rec- ommendation to avoid ambiguity, especially because pres- sure effects are large, e.g., multiplying fCO₂ by 2 and pCO₂ by 18 when at 5000 m. The mocsy 2.0 package offers the first choice as a default and the second as an option. Other public packages compared by Orr et al. (2014) offer neither option, computing fCO₂ and pCO₂ at 1 atm total pressure with in situ temperature. A second package, seacarb 3.0.6, also now offers all three approaches to compute subsurface fCO₂ and pCO₂ (Gattuso et al., 2015; Orr et al., 2015).

For simplicity, modelers typically make ocean carbonate chemistry calculations assuming (1) that model density is constant, (2) that simulated potential temperature is equiva- lent to in situ temperature, and (3) that model depth is equiv- alent to pressure. None of the three approximations produces significant errors at the surface. Yet errors in computed vari- ables grow with depth. At 5000 m, the sum of the relative errors from the three approximations reach more than 3% for pCO₂, Ω<sub>A</sub>, and H⁺ (pH shift of up to 0.015). Yet they re- main at less than 1% for the three inorganic carbon species (<1 µmol kg⁻¹). The mocsy modeling routines avoid these errors with little additional coding and trivial increases in computational time.

The same code also offers two new options that have be- come available since the publication of the best-practices guide (Dickson et al., 2007). Those options concern an as- sessment that seawater contains 4% more total boron than thought previously (Lee et al., 2010) and new formula- tions for K₁ and K₂ designed to include low-salinity waters (Millero, 2010). The new boron option leads to substantial
<figure><img src="image_17.png" /><figcaption>**Figure 12.** Changes in the aragonite saturation horizon (ASH) due to (top) using the new formulation of total boron (Lee et al. minus Uppström), (middle) accounting for nutrient alkalinity (relative to the common simplification where it is neglected), and (bottom) summing both corrections. Changes in ASH (m) are smoothed over 5° bands of latitude and shown as zonal means over the Atlantic, Pacific, and Indian Oceans. The peaks at 30 and 12° S in the Atlantic are caused by subtle shifts in horizontal gradients of Ω<sub>A</sub>.</figcaption></figure>

shifts in computed surface variables, e.g., +4 to +6 µatm in pCO₂, -1.2 to -2.5 µmol kg⁻¹ in CO₃²⁻, and -0.006 in pH. Comparable shifts at depth lead to a shallower computed ASH by 50 m in the North Atlantic and by 90 m in the Southern Ocean. The new option for K₁ and K₂ leads to a shift of +7 µatm in surface pCO₂ in the Southern Ocean. When both options are combined, the Southern Ocean's surface pCO₂ becomes 3 µatm higher than in the high northern latitudes and 6 µatm higher than in the tropics. Yet despite their large effects, we cannot yet recommend either new option. Before using the new B<sub>T</sub> formulation, the community needs to assess how K₁ and K₂ will change because those measurements depend on B<sub>T</sub>. And before using the latest K₁ and K₂ formulations designed to include low-salinity waters, their coefficients should be refit and given with enough precision to produce consistent results among pH scales.

For simplicity, models often neglect contributions of total phosphorus and silicon to total alkalinity. Doing so biases computed variables, e.g., shifting Southern Ocean surface pCO₂ by +6 µatm and the ASH below by up to +200 m. Because biases are not uniform, they also affect lateral and vertical gradients. For greater consistency, we invite all mod-

elers to account for nutrient alkalinity in carbonate system calculations and to avoid the three common model approximations for subsurface calculations. The mocsy 2.0 package removes these limitations while following best practices. It also adopts the fastest and safest method to solve the total alkalinity-pH equation (Munhoven, 2013). That solver always converges, even under extreme conditions, e.g., in estuaries subject to intense freshwater fluxes where other solvers may fail.

## Code availability

The mocsy package is distributed under the MIT license and is available from GitHub. It can be retrieved in the conventional manner with Git, namely by issuing the command

git clone git://github.com/jamesorr/mocsy.git

in an X terminal on Linux, Mac, or PC operating systems. If that fails, install Git and try again, or go to the main web page on GitHub (https://github.com/jamesorr/mocsy) and click on the link "Download ZIP". Once downloaded, mocsy can be compiled under Linux or Unix by typing

make.

More details on the code, its compilation, and examples of its use in Fortran or when called from Python can be found in the mocsy manual at http://ocmip5.ipsl.jussieu.fr/mocsy/.

**Acknowledgements.** We thank the original OCMIP code developers (C. L. Sabine, R. M. Key, and R. G. Najjar) and those who provided subsequent corrections or comments (A. Mouchet, K. Lindsay, J. Christian, and G. McKinley). We are grateful to J.-P. Gattuso, E. Lewis, and D. Wallace, as well as S. van Heuven, D. Pierrot, and J. Rae for coding and making public their packages (seacarb, CO2SYS, and CO2SYS-MATLAB), which have been no small inspiration for mocsy. We also thank P. Brockmann and Y. Meurdesoif for professional advice on coding and distribution strategies. The original manuscript and code have been improved because of three dedicated referees: G. Munhoven shared scientific and technical savvy while making his stable, efficient, and well-documented pH solver (SolveSAPHE) publicly available, allowing us to easily integrate it in mocsy; A. Dickson offered much critical expertise on carbonate chemistry, for instance inspiring us to include pressure corrections that dramatically alter computed subsurface fCO₂ and pCO₂; and an anonymous referee provided a wider modeler's perspective. We are grateful to F. Perez for pointing out that formulations of K₁ and K₂ depend upon B<sub>T</sub>. The OCMIP code was first modified with support from the EU FP5 project NOCES (grant EVK2-CT-2001-00134) and the FP6 project CARBOOCEAN (grant 511176(GOCE)). More recent efforts to refine, document, and make the code publicly available were supported by the French ANR project MACROES (MACRoscope for Oceanic Earth System ANR-09-CEP-003) and two EU FP7 projects, CARBOCHANGE (grant 264879) and MedSeA (grant 265103).
Edited by: D. Roche

References

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

Aumont, O. and Bopp, L.: Globalizing results from ocean in situ iron fertilization studies, Global Biogeochem. Cy., 20, GB2017, doi:10.1029/2005GB002591, 2006.

Bates, N. R., Orchowska, M. I., Garley, R., and Mathis, J. T.: Summertime calcium carbonate undersaturation in shelf waters of the western Arctic Ocean – how biological processes exacerbate the impact of ocean acidification, Biogeosciences, 10, 5281–5309, doi:10.5194/bg-10-5281-2013, 2013.

Brewer, P. and Peltzer, E.: Limits to marine life, Science, 324, 347–348, doi:10.1126/science.1170756, 2009.

Caldeira, K. and Wickett, M. E.: Anthropogenic carbon and ocean pH, Science, 425, p. 365, 2003.

Cocco, V., Joos, F., Steinacher, M., Frölicher, T. L., Bopp, L., Dunne, J., Gehlen, M., Heinze, C., Orr, J., Oschlies, A., Schneider, B., Segschneider, J., and Tjiputra, J.: Oxygen and indicators of stress for marine life in multi-model global warming projections, Biogeosciences, 10, 1849–1868, doi:10.5194/bg-10-1849-2013, 2013.

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

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

Dickson, A. G. and Goyet, C.: Handbook of Methods for the Analysis of the Various Parameters of the Carbon Dioxide System in Seawater, version 2, CDIAC-74, ORNL, 187 pp., 1994.

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

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

Doney, S. C., Lindsay, K., Fung, I., and John, J.: Natural variability in a stable, 1000-yr global coupled climate-carbon cycle simulation, J. Climate, 19, 3033–3054, doi:10.1175/JCLI3783.1, 2006.

Dufour, C. O., Le Sommer, J., Gehlen, M., Orr, J. C., Molines, J.-M., Simeon, J., and Barnier, B.: Eddy compensation and controls of the enhanced sea-to-air CO₂ flux during positive phases of the Southern Annular Mode, Global Biogeochem. Cy., 27, 950–961, doi:10.1002/gbc.20090, 2013.

Follows, M. J., Ito, T., and Dutkiewicz, S.: On the solution of the carbonate chemistry system in ocean biogeochemistry models, Ocean Model., 12, 290–301, doi:10.1016/j.ocemod.2005.05.004, 2006.

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

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

Gattuso, J.-P., Epitalon, J.-M., and Lavigne, H.: seacarb: seawater carbonate chemistry with R, R package version 3.0.6, The Comprehensive R Archive Network, available at: http://CRAN.R-project.org/package=seacarb, 2015.

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

Lavigne, H. and Gattuso, J.-P.: seacarb: seawater carbonate chemistry with R. R package version 2.4.2, in: The Comprehensive R Archive Network, available at: http://wwwobs-vlfr.fr/~gattuso/files/Lavigne_Gattuso_2011_seacarb_2.4.2.pdf (last access: 24 April 2014), 2011.

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

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

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

Lueker, T. J., Dickson, A. G., and Keeling, C. D.: Ocean pCO₂ 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, 2000.

Maier-Reimer, E.: Geochemical cycles in an ocean general circulation model: preindustrial tracer distributions, Global Biogeochem. Cy., 7, 645–677, 1993.

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

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

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

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

Morris, A. and Riley, J.: The bromide/chlorinity and sulphate/chlorinity ratio in sea water, in: Deep-Sea Research and Oceanographic Abstracts, 13, 699–705, 1966.
Mucci, A.: The solubility of calcite and aragonite in seawater at various salinities, temperatures, and one atmosphere total pressure, Am. J. Sci., 283, 780–799, 1983.

Müller, S. A., Joos, F., Plattner, G.-K., Edwards, N. R., and Stocker, T. F.: Modeled natural and excess radiocarbon: sensitivities to the gas exchange formulation and ocean transport strength, GBC, 22, GB3011, doi:10.1029/2007GB003065, doi:10.1029/2007GB003065, 2008.

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

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

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

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

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

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

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

Proye, A. and Gattuso, J.-P.: Seacarb, an R package to calculate parameters of the seawater carbonate system, available at: http://www obs-vlfr fr / ~ gattuso / files / Proye_Gattuso_2003 _seacarb_0.93.pdf (last access: 24 April 2014), 2003.

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

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

Sarmiento, J. L., Monfray, P., Maier-Reimer, E., Aumont, O., Murane, R., and Orr, J. C.: Sea–air CO₂ fluxes and carbon transport: a comparison of three ocean general circulation models, Global Biogeochem. Cy., 14, 1267–1281, 2000.

Uppström, L. R.: The boron/chlorinity ratio of deep-sea water from the Pacific Ocean, in: Deep-Sea Research and Oceanographic Abstracts, vol. 21, Elsevier, 161–162, 1974.

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

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

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