Improved Routines to Model the Ocean Carbonate System: Mocsy 2.0

Modelers compute ocean carbonate chemistry often based on code from the Ocean Carbon Cycle Model In-tercomparison 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 alka-linity 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 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 −1). Yet by default, mocsy's deep-water f CO 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.


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 θ , 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 Published by Copernicus Publications on behalf of the European Geosciences Union.
(mol kg −1 ), 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 −3 ), 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 2 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.

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 2− 3 concentrations divided by the apparent solubility product (either K A or K C ) yields the saturation state (i.e., for aragonite A or for calcite C , respectively).

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 highresolution, 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 Ob estuary, 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 2 S acid systems.

From aqueous CO 2 to pCO 2
In mocsy, we use K 0 from Weiss (1974) to compute f CO 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 f CO 2 : For this calculation of oceanic pCO 2 , mocsy does not use the combined coefficient F from Weiss and Price (1980), where 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): where P is the total pressure (atmospheric + hydrostatic, both in units of atm), vCO 2 is the partial molal volume of CO 2 (32.3 cm 3 mol −1 ), R is the gas constant (82.05736 cm 3 atm mol −1 K −1 ), and T is the absolute temperature (K).Thus the computed f CO 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): 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 1), and δ 12 = 57.7 − 0.118 T .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 f CO 2 and pCO 2 may be considered as pseudo-potential values that a water parcel would 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 f CO 2 and pCO 2 also at surface pressure but using the potential temperature θ in place of T in Eqs. ( 3) and (4).That is, these true potential f CO 2 and pCO 2 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 f CO 2 and pCO 2 with the first approach, at in situ temperature and total pressure.
The effects of these pressure corrections on K 0 and C f are shown in terms of their multiplicative effect on f CO 2 and pCO 2 .For K 0 , 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 f CO 2 and 1.05 for pCO 2 .At 5000 m, the total multiplicative effect reaches more than a factor of 2 for f CO 2 and 18 for pCO 2 .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 f CO 2 and pCO 2 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 0 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).

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 .
To do so, mocsy 2.0 replaces the analytical formula for R f from Frankignoulle (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 f CO 2 in place of pCO 2 .Secondly, in CO2SYS the step size h for the numerical approximation of the derivative (df/dx ≈ f (x 0 + h) − f (x 0 − h) /2h, where x is C T and f (x) is pCO 2 ) is 10 times larger than ours (h = 0.1 µmol kg −1 ), 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.

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 1 and K 2 , 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 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.

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.

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 3. Global-mean vertical profiles of carbonate system variables computed with mocsy 2.0 from the data described in Fig. 2. salinity (Antonov et al., 2010), and concentrations of P T and Si T (Garcia et al., 2010).

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 CO 2 and pCO 2 continue to increase with depth, particularly the latter, because of pressure corrections to K 0 and C f , respectively.Conversely, the 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 CO 2− 3 profile mirrors that for CO * 2 , as expected from the zonal-mean distributions.Vertical distributions of the saturation states A and C generally follow the CO 2− 3 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 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.

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 CO * 2 , HCO − 3 , CO 2− 3 , 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 −1 for each of the three inorganic carbon species, and within 0.0002 nmol kg −1 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 A and 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 remains 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.Differences in surface zonal means between mocsy and CO2SYS for computed variables shown as (top) percent relative differences (100 (V mocsy − V co2sys )/V co2sys ) and as (bottom) absolute differences (V mocsy − V co2sys ).The absolute differences are given in µmol kg −1 for the three inorganic carbon species, in nmol kg −1 for H + , and in µatm for f CO 2 and pCO 2 ; other computed variables are unitless.The two remaining computed variables, surface f CO 2 and pCO 2 , both track CO * 2 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 f CO 2 is twice as large as that in CO2SYS, while its pCO 2 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.

Model approximations
Errors from all three model approximations increase with depth (Fig. 6).Maximum total errors reach 3 % or more for pCO 2 , H + , and A at 5000 m.Yet the causes differ.For 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.
pCO 2 , the largest error comes from the pressure approximation, while errors from the other two approximations are smaller and nearly compensate one another.For 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 CO 2 in the deep ocean is about half that for pCO 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 A have somewhat similar patterns to that of CO 2− 3 .Yet at 5000 m, the magnitude of the total relative error of A reaches 4 %, more than 5 times greater than for CO 2− 3 .That relative enhancement of the former over the latter is partly due to the sharp decline in A , e.g., below 1000 m, whereas CO 2− 3 concentrations remain relatively stable (Fig. 3).The dominance of the pressure-induced increase in K A largely outweighs its sensitivity to temperature, unlike for CO 2− 3 .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 pCO 2 , H + , and f CO 2 ).Yet total absolute errors generally remain small, even at 5000 m, e.g., less than 0.6 µmol kg −1 for the three inorganic carbon species, less than 0.02 for 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 pCO 2 it reaches 350 µatm.

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- face CO 2− 3 everywhere, from −1.5 µmol kg −1 in the Southern Ocean to −3 µmol kg −1 in the tropics (Fig. 8).The corresponding decline in surface A is between 0.02 and 0.04.Simultaneously, computed surface pCO 2 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 R f decreases because pCO 2 increases (C T 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 2− 3 of −1 ± 0.1 µmol kg −1 .There is also an exponential increase in pCO 2 up to 110 µatm at 5000 m (Fig. 9).The latter represents a 1 % increase in the baseline pCO 2 (Fig. 3).
The new option where the formulations of K 1 and K 2 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 T .Changes from the new option for K 1 and K 2 are larger in the Southern Ocean: pCO 2 increases by up to +7 µatm, corresponding to a 2 % increase in CO * 2 , while pH declines by up to 0.006.In contrast, in the lowest latitudes corresponding changes in pCO 2 and CO * 2 are negligible.In the Southern Ocean, changes in pCO 2 , CO * 2 , and pH from using the new K 1 and K 2 option reinforce those from the new B T option.But in the same region, changes from the two options partly compensate one another for CO 2− 3 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 2 (+12 µatm), CO * 2 (+4 %), pH (−0.012), and R f (−0.30), whereas they are largest in the lower latitudes for CO 2− 3 (−3.5 µmol kg −1 ) and A (−0.06).With depth, mean absolute changes due to both options grow for pCO 2 and CO * 2 , change little for pH and R f , and decline for CO 2− 3 and A .These differences appear greatly affected by the structure of the baseline profiles (Fig. 3).

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 T and Si T (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 noncarbonate alkalinity and hence total alkalinity,  2010) formulation for the total boron-to-salinity ratio (blue B), (2) the Millero (2010) formulations for K 1 and K 2 (red, K 1 K 2 ), and (3) the Dickson and Riley (1979) formulation for K F (green K F ). 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). where 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 W provides an indicator of how a given seawater sample differs from acid-base neutrality.
As an input variable, A T 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 T is unaffected by our calculations.Yet A C is affected, being computed by difference.As A B , A P , or A Si increases, computed A C must decrease.In the surface ocean, nutrient alkalinity (A P +A Si ) 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 2 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 2− 3 changes by −1.6 µmol kg −1 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 2 differ between some models by only 0.1 Petagram of carbon per year (Pg C yr −1 ), equivalent to less than 1 µatm in the air-sea difference in pCO 2 (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 2 .
These biases concern models that neglect contributions of phosphoric and silicic acids to alkalinity (by assuming A T = A C + A B + A W ). 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 ( A < 1) that is too late, and a simulated interhemispheric north-tosouth 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  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 in the surface ocean (Figs. 8 and 10).In the deep ocean, the average combined effect shifts CO 2− 3 by −3 µmol kg −1 , pH by −0.018, and pCO 2 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 T nor that for K 1 and K 2 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 1 and K 2 that predate the Lee et al. (2010) study and depend themselves on B T (Mehrbach et al., 1973, Eq. 8).Likewise, the option to use new formulations of K 1 and K 2 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).

Conclusions
Although f CO 2 and pCO 2 are typically measured or calculated at the surface, e.g., to estimate air-sea CO 2 fluxes, a limited number of studies also focus on those quantities below 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 f CO 2 and pCO 2 after making pressure corrections to K 0 and the fugacity coefficient following Weiss (1974) using total in situ pressure and in situ temperature or (2) calculate potential f CO 2 and pCO 2 at 1 atm pressure while using potential instead of in situ temperature.Results should then be clearly labeled as either in situ or potential f CO 2 and pCO 2 .We make this recommendation to avoid ambiguity, especially because pressure effects are large, e.g., multiplying f CO 2 by 2 and pCO 2 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 f CO 2 and pCO 2 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 f CO 2 and pCO 2 (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 equivalent to in situ temperature, and (3) that model depth is equivalent to pressure.None of the three approximations produces significant errors at the surface.Yet errors in computed variables grow with depth.At 5000 m, the sum of the relative errors from the three approximations reach more than 3% for pCO 2 , A , and H + (pH shift of up to 0.015).Yet they remain at less than 1 % for the three inorganic carbon species (<1 µmol kg −1 ).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 become available since the publication of the best-practices guide (Dickson et al., 2007).Those options concern an assessment that seawater contains 4 % more total boron than thought previously (Lee et al., 2010) and new formulations for K 1 and K 2 designed to include low-salinity waters (Millero, 2010).The new boron option leads to substantial shifts in computed surface variables, e.g., +4 to +6 µatm in pCO 2 , −1.2 to −2.5 µmol kg −1 in CO 2− 3 , 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 1 and K 2 leads to a shift of +7 µatm in surface pCO 2 in the Southern Ocean.When both options are combined, the Southern Ocean's surface pCO 2 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 T formulation, the community needs to assess how K 1 and K 2 will change because those measurements depend on B T .And before using the latest K 1 and K 2 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 2 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.
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/.

Figure 1 .
Figure1.Effect of pressure on K 0 and C f given as factors that would need to be multiplied by (left) f CO 2 and (right) pCO 2 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.

Figure 2 .
Figure2.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 .

Figure 4 .
Figure 4. Differences in surface zonal means between mocsy and CO2SYS for computed variables shown as (top) percent relative differences (100 (V mocsy − V co2sys )/V co2sys ) and as (bottom) absolute differences (V mocsy − V co2sys ).The absolute differences are given in µmol kg −1 for the three inorganic carbon species, in nmol kg −1 for H + , and in µatm for f CO 2 and pCO 2 ; other computed variables are unitless.

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

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

Figure 9 .
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.

Figure 10 .
Figure10.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.

Figure 11 .
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.

Figure 12 .
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 A .