Temperature effect on seawater fCO₂ revisited: theoretical basis, uncertainty analysis and implications for parameterising carbonic acid equilibrium constants

The sensitivity of the fugacity of carbon dioxide in seawater (fCO₂) to temperature (denoted υ, reported in % °C⁻¹) is critical for the accurate fCO₂ measurements needed to build global carbon budgets and for understanding the drivers of air–sea CO₂ flux variability across the ocean. However, understanding and computing υ have been restricted to either using empirical functions fitted to experimental data or determining it as an emergent property of a fully resolved marine carbonate system, and these two approaches are not consistent with each other. The lack of a theoretical basis and an uncertainty estimate for υ has hindered resolving this discrepancy. Here, we develop a new approach for calculating the temperature sensitivity of fCO₂ based on the equations governing the marine carbonate system and the van 't Hoff equation. This shows that, to first order, ln (fCO₂) should be proportional to $\mathrm{1}/t_{\mathrm{K}}$ (where t_K is temperature in kelvin), rather than to temperature, as has previously been assumed. This new approach is, to first order, consistent with calculations from a fully resolved marine carbonate system, which we have incorporated into the PyCO2SYS software. Agreement with experimental data is less convincing but remains inconclusive due to the scarcity of direct measurements of υ, particularly above 25 °C. However, the new approach is consistent with field data, performing better than any other approach for adjusting fCO₂ by up to 10 °C if spatiotemporal variability in its single fitted coefficient is accounted for. The uncertainty in υ arising from only measurement uncertainty in the main experimental dataset where υ has been directly measured is in the order of 0.04 % °C⁻¹, which corresponds to a 0.04 % uncertainty in fCO₂ adjusted by +1 °C. However, spatiotemporal variability in υ is several times greater than this, so the true uncertainty due to the temperature adjustment in fCO₂ adjusted by +1 °C using the most widely used constant υ value is around 0.24 %. This can be reduced to around 0.06 % using the new approach proposed here, and this could be further reduced by more measurements. The spatiotemporal variability in υ arises mostly from the equilibrium constants for CO₂ solubility and carbonic acid dissociation ($K_{\mathrm{1}}^{\ast }$ and $K_{\mathrm{2}}^{\ast }$), and its magnitude varies significantly depending on which parameterisation is used for $K_{\mathrm{1}}^{\ast }$ and $K_{\mathrm{2}}^{\ast }$. Seawater fCO₂ can be measured accurately enough that additional experiments should be able to detect spatiotemporal variability in υ and distinguish between different parameterisations for $K_{\mathrm{1}}^{\ast }$ and $K_{\mathrm{2}}^{\ast }$. Because the most widely used constant υ was coincidentally measured from seawater with roughly global average υ, our results are unlikely to significantly affect global air–sea CO₂ flux budgets, but they may have more important implications for regional budgets and studies that adjust by larger temperature differences.

1 Introduction

The surface ocean continuously exchanges carbon dioxide (CO₂) with the atmosphere. The net direction of exchange varies through space and time, with some ocean areas acting as CO₂ sources and others as sinks (Takahashi et al., 2009). Whether an area acts as a CO₂ source or sink depends on the difference between the fugacity of CO₂ (fCO₂) in the sea and in the air (denoted $\mathrm{\Delta }f{\mathrm{CO}}_{\mathrm{2}}=f{\mathrm{CO}}_{\mathrm{2}}\left(\mathrm{sea}\right)-f{\mathrm{CO}}_{\mathrm{2}}\left(\mathrm{air}\right)$), with net CO₂ transfer from high to low fCO₂. The fugacity of CO₂ is similar to its partial pressure (pCO₂) but with a small correction for the non-ideal behaviour of CO₂ in air; fCO₂ is typically about 99.7 % of pCO₂. The gas exchange rate is modulated by a set of processes that are commonly parameterised as a function of wind speed (Wanninkhof, 2014). Most of the surface ocean is either a CO₂ source or sink, rather than being in equilibrium with the atmosphere, because many of the processes that modify seawater fCO₂ (e.g. seasonal cycles of temperature and primary productivity) operate on faster timescales than air–sea CO₂ equilibration, which takes from months to a year for the surface mixed layer (Jones et al., 2014).

On top of the natural background variability, the global ocean is a net CO₂ sink due to rising atmospheric fCO₂. About a quarter of the anthropogenic CO₂ emitted each year, some 2.8 ± 0.4 Gt C in 2022, is taken up by the ocean, mitigating the greenhouse effect of CO₂ on Earth's climate (Friedlingstein et al., 2023). Together, uptake of CO₂ and warming drive changes in the marine carbonate system, altering biologically important variables like seawater pH and the saturation states with respect to carbonate minerals (Humphreys, 2017). Quantifying marine CO₂ uptake requires highly accurate fCO₂ measurements: it is driven by a global mean ΔfCO₂ of −4.1 ± 1.6 µatm (Fay et al., 2024), which is around 2 % of mean fCO₂. Consequently, small systematic biases in seawater fCO₂ propagate into large errors in ΔfCO₂ and, thus, the total ocean CO₂ sink. The required measurement accuracy has been quantified by the Global Ocean Acidification Observing Network (GOA-ON), with the network suggesting a minimum accuracy of 0.5 % (i.e. around 2 µatm) for “climate quality” seawater fCO₂ measurements that can be used to detect decadal trends or 2.5 % (∼10 µatm) for “weather quality” fCO₂ data to investigate changes of greater magnitude, such as seasonal cycles (Newton et al., 2015).

The observational constraint on marine CO₂ uptake is based on measurements of surface ocean fCO₂ from ships, compiled into data products such as the Surface Ocean CO₂ Atlas (SOCAT; Bakker et al., 2016). Surface seawater is equilibrated with an enclosed headspace, which flows through to a measuring device, and the raw measurements are passed through a series of processing steps to obtain seawater fCO₂. One important processing step is the temperature adjustment. While it is being pumped from the surface ocean to the equilibration chamber, seawater may warm up or cool down. As fCO₂ is highly sensitive to temperature (Weiss, 1974), each measurement must be adjusted back to its in situ temperature in order to obtain environmentally relevant values of fCO₂ (Pierrot et al., 2009; Pfeil et al., 2013; Goddijn-Murphy et al., 2015). For the data in SOCAT with the most stringent quality control flags (A and B), a maximum change in temperature (Δt) of 1 °C between intake and measurement is permitted (although adjustments are usually around 0.2–0.3 °C in practice; Wanninkhof et al., 2022a), in order to minimise uncertainty in fCO₂ due to the temperature adjustment. Adjustments with greater Δt are often used in other applications. Some groups measure fCO₂ in discrete samples at a constant temperature, which must then be adjusted by up to 20 °C to return to in situ conditions (Wanninkhof et al., 2022a). Several studies have normalised seawater fCO₂ observations to a constant temperature in order to investigate the non-thermal (i.e. chemical) drivers of fCO₂ change on various timescales (e.g. Takahashi et al., 2002, 2003, 2006; Feely et al., 2006). Others have computed the temperature sensitivity in order to isolate the thermal component of observed fCO₂ change and unravel the drivers of fCO₂ variability (e.g. Metzl et al., 2010; Fassbender et al., 2022; Rodgers et al., 2023). Data from SOCAT have been adjusted from the intake temperature, typically measured a few metres below the sea surface, to the temperature directly at the sea surface in order to assess how the “cool skin” effect affects computed air–sea CO₂ fluxes (e.g. Watson et al., 2020; Dong et al., 2022). Finally, climatological fCO₂ datasets have been adjusted to match different sea surface temperature fields for specific analyses (e.g. Kettle et al., 2009; Land et al., 2013), and this functionality is part of published data processing toolkits (e.g. Shutler et al., 2016; Holding et al., 2019).

The most widely used approach to adjust fCO₂ data to different temperatures (t), used in all of the examples above, is based on data from Takahashi et al. (1993), who measured the variation in fCO₂ with temperature in one seawater sample from the North Atlantic Ocean (strictly speaking, they reported pCO₂, but the results apply equally to fCO₂). Their measurements indicated that the temperature sensitivity of fCO₂ (denoted υ; Eq. 1) was essentially independent of fCO₂ and temperature (i.e. temperature and ln (fCO₂) have a linear relationship). This apparent near-linearity was consistent with earlier measurements by Kanwisher (1960) on a single seawater sample collected from the vicinity of Woods Hole (Massachusetts, USA) and with later measurements by Lee and Millero (1995). Takahashi et al. (1993) also claimed that the value of υ was no different for higher-salinity waters, citing unpublished preliminary data.

However, it has long been known that υ is not constant across the ocean. For example, Gordon and Jones (1973) computed υ using a finite-difference approach across a range of different temperature, salinity, alkalinity and pH values and found that the variability in υ could be approximately parameterised as a function of pCO₂. Using a similar approach but with different parameterisations of the equilibrium constants in their marine carbonate system calculations, Weiss et al. (1982) and Copin-Montegut (1988) each fitted υ to functions of fCO₂ and temperature. It has been noted that υ may vary with the dissolved inorganic carbon (T_C) content of seawater (e.g. Takahashi et al., 2002) or, more generally, with the ratio between T_C and total alkalinity (A_T) (e.g. Goyet et al., 1993; Wanninkhof et al., 2022a).

There are also several lines of evidence that the relationship between temperature and ln (fCO₂) is not linear (i.e. υ is not constant for a single, isochemical sample). Takahashi et al. (1993) acknowledged this possibility by fitting ln (fCO₂) to a quadratic function of temperature as well as their linear function, although they did not give any reason to choose one over the other. The fits of Weiss et al. (1982), Copin-Montegut (1988) and Goyet et al. (1993) included various degrees of polynomial temperature dependence. Calculating how fCO₂ changes with temperature using software tools like CO2SYS (Lewis and Wallace, 1998) also shows a non-linear relationship, but its exact form and how much it deviates from linearity both vary widely depending on which parameterisation for the equilibrium constants of carbonic acid is used (McGillis and Wanninkhof, 2006). Data compiled from multiple cruises during which fCO₂ was measured at two different temperatures in an equivalent pair of samples, often with accompanying A_T and T_C, showed that a linear (constant υ) adjustment was insufficiently accurate for Δt greater in magnitude than about 1 °C (Wanninkhof et al., 2022a).

Wanninkhof et al. (2022a) recommended calculating υ from A_T and T_C in order to make accurate adjustments with greater Δt, using the carbonic acid dissociation constants of Lueker et al. (2000) and the borate-to-chlorinity ratio of Uppström (1974). However, the requirement for a second marine carbonate system parameter (other than fCO₂) is a major disadvantage because, due to ships of opportunity and the variety of (semi-)autonomous sensors available, seawater fCO₂ is by far the most widely measured marine carbonate system parameter in the ocean and is often the only one available at a particular observation point. A possible workaround would be to use a parameterisation to generate values for a second parameter. For example, across much of the open ocean, A_T can be calculated from temperature and salinity with an uncertainty of less than 10 µmol kg⁻¹ (Lee et al., 2006). Under typical conditions, this component alone would propagate through to a less than 0.1 µatm uncertainty in fCO₂ adjusted by 10 °C (calculated with PyCO2SYS; Humphreys et al., 2022). However, this ignores the much larger contribution of uncertainties in the equilibrium constants used in the calculations. These are not well constrained, but using the set suggested by Orr et al. (2018) makes the uncertainty in the same adjusted fCO₂ value several hundred times greater, at around 15 µatm. An adjustment method based on direct measurements of fCO₂ and independent from a second marine carbonate system parameter, avoiding the uncertainties inherent in calculations involving equilibrium constants, is therefore highly desirable.

Given the evidence that the equations of Takahashi et al. (1993) are not sufficiently accurate for the full range of temperature adjustments of fCO₂ required by ocean scientists, why do they remain so widely used? First, they are straightforward to calculate, relying only on variables that will be available with any fCO₂ measurement (i.e. fCO₂ itself and temperature). The more complex existing parameterisations that also fulfil this criterion (Copin-Montegut, 1988; e.g. Gordon and Jones, 1973; Weiss et al., 1982) are not based on up-to-date parameterisations of the carbonic acid dissociation constants, which can lead to differences in adjusted fCO₂ values of over 30 µatm for larger Δt (McGillis and Wanninkhof, 2006). However, a more fundamental impediment is not knowing what form the relationship between temperature and ln (fCO₂) should be expected to take from a theoretical perspective. Should it be (roughly) linear, with deviations from linearity thus indicating problems with the calculations used to solve the marine carbonate system, or should it be non-linear, thus indicating some deficiency with the measurements of Takahashi et al. (1993)? With respect to the latter, this then poses the follow-up question of which of the non-linear relationships obtained from the different parameterisations of the carbonic acid dissociation constants should be followed.

Here, we aim to provide this missing theoretical basis by developing a new functional form regarding how fCO₂ and, thus, υ vary with temperature, derived from the equations governing the marine carbonate system and the van 't Hoff equation. We compare this new form with υ calculated from a fully resolved marine carbonate system – to which end, we have added the calculation of υ to PyCO2SYS (Humphreys et al., 2022) – and with laboratory and field data. We compute uncertainties in υ and temperature-adjusted fCO₂ values due to experimental uncertainties and variability in υ through space and time. We assess the contributions of the different equilibrium constants of the marine carbonate system to υ, especially those for carbonic acid dissociation, and consider the implications for how these constants are parameterised and how these previously separate aspects can be better integrated in future studies for a more accurate and internally consistent understanding of the marine carbonate system.

2 Methods

2.1 Adjusting fCO₂ to a different temperature

2.1.1 General theory

The temperature sensitivity of seawater fCO₂ (υ, in units of % °C⁻¹) is the first derivative of the natural logarithm of fCO₂ with respect to temperature (t, in °C) under constant total alkalinity (A_T) and dissolved inorganic carbon (T_C) (see Humphreys et al., 2022, and references therein for more detailed definitions of A_T, T_C and fCO₂):

$$\begin{array}{}\text{(1)}& \mathit{\upsilon }={\left.{\displaystyle \frac{\partial \ln \left(f{\mathrm{CO}}_{\mathrm{2}}\right)}{\partial t}}\right|}_{A_{\mathrm{T}},T_{\mathrm{C}}}.\end{array}$$

We note that, here and throughout, fCO₂ must be divided by the unit of fCO₂ (e.g. 1 µatm if fCO₂ is in µatm) before taking the logarithm, and “constant A_T and T_C” implies that all other parameters but temperature are also constant (e.g. salinity and total nutrient contents). We report υ in percent per degree Celsius (% °C⁻¹; i.e. 10⁻² °C⁻¹).

To adjust an fCO₂ value from one temperature (t₀) to another (t₁), υ must be integrated across the relevant temperature range:

$$\begin{array}{}\text{(2)}& \ln \left[f{\mathrm{CO}}_{\mathrm{2}}\right(t_{\mathrm{1}}\left)\right]=\ln \left[f{\mathrm{CO}}_{\mathrm{2}}\right(t_{\mathrm{0}}\left)\right]+\mathrm{{\rm Y}},\end{array}$$

where Υ denotes the integral

$$\begin{array}{}\text{(3)}& \mathrm{{\rm Y}}={\int }_{t_{\mathrm{0}}}^{t_{\mathrm{1}}}\mathit{\upsilon }\mathrm{d}t.\end{array}$$

Equivalent to Eq. (2),

$$\begin{array}{}\text{(4)}& f{\mathrm{CO}}_{\mathrm{2}}\left(t_{\mathrm{1}}\right)=f{\mathrm{CO}}_{\mathrm{2}}\left(t_{\mathrm{0}}\right)\times \exp \left(\mathrm{{\rm Y}}\right).\end{array}$$

Thus, in order to adjust the temperature of an fCO₂ value using Eq. (4), a function for υ that can be integrated with respect to temperature is needed.

2.1.2 Experiments to measure υ

Takahashi et al. (1993) attempted to determine υ by measuring how pCO₂ varies with temperature (eight measurements at seven temperatures from 2.1 to 24.5 °C) in one seawater sample with constant A_T (reported as “2270 to 2295” µmol kg⁻¹) and T_C (2074 µmol kg⁻¹). The practical salinity was 35.380. Phosphate and silicate concentrations were reported as zero while noting the possibility of silicate contamination from the sample bottles. The sample had been part of the laboratory intercomparison exercise reported by Poisson et al. (1990).

We determined a more precise A_T for the sample by finding the value that, when paired with the T_C and other parameters above, returned the least-squares best fit to the reported pCO₂ values. When using the Lueker et al. (2000) carbonic acid constants and the borate-to-chlorinity ratio of Uppström (1974), the A_T obtained with PyCO2SYS (Humphreys et al., 2022) was 2263 ± 3 µmol kg⁻¹, although each set of carbonic acid constants returned its own unique A_T. For consistency with modern measurements and the rest of our analysis, we also converted the Takahashi et al. (1993) pCO₂ values into fCO₂, although this makes no meaningful difference to any of our results.

Takahashi et al. (1993) fitted their pCO₂ dataset to linear (subscript l) and quadratic (subscript q) equations (here written in terms of fCO₂):

$$\begin{array}{}\text{(5)}& {\displaystyle }\ln \left(f{\mathrm{CO}}_{\mathrm{2}}\right)=b_{\mathrm{l}}t+c_{\mathrm{l}},\text{(6)}& {\displaystyle }\ln \left(f{\mathrm{CO}}_{\mathrm{2}}\right)=a_{\mathrm{q}}{t}^{\mathrm{2}}+b_{\mathrm{q}}t+c_{\mathrm{q}},\end{array}$$

where a, b and c are the fitted coefficients. Differentiating with respect to t yields the following for υ:

$$\begin{array}{}\text{(7)}& {\displaystyle }{\mathit{\upsilon }}_{\mathrm{l}}=b_{\mathrm{l}},\text{(8)}& {\displaystyle }{\mathit{\upsilon }}_{\mathrm{q}}=\mathrm{2}{a}_{\mathrm{q}}t+b_{\mathrm{q}}.\end{array}$$

Finally, integrating across the relevant temperature range yields Υ, following, for example, Takahashi et al. (2009):

$$\begin{array}{}\text{(9)}& {\displaystyle }{\mathrm{{\rm Y}}}_{\mathrm{l}}=b_{\mathrm{l}}(t_{\mathrm{1}}-t_{\mathrm{0}}),\text{(10)}& {\displaystyle }{\mathrm{{\rm Y}}}_{\mathrm{q}}=a_{\mathrm{q}}(t_{\mathrm{1}}^{\mathrm{2}}-t_{\mathrm{0}}^{\mathrm{2}})+b_{\mathrm{q}}(t_{\mathrm{1}}-t_{\mathrm{0}}).\end{array}$$

Takahashi et al. (1993) reported values of b_l=4.23 % °C⁻¹, $a_{\mathrm{q}}=-\mathrm{4.35}\times {\mathrm{10}}^{-\mathrm{5}}$ °C⁻² and b_q=4.33 % °C⁻¹ with a root-mean-square deviation (RMSD) of 1.0 µatm in fCO₂ for both fits. Reproducing these fits with their reported data gives an RMSD of 1.2 µatm for the linear fit and 0.9 µatm for the quadratic.

Lee and Millero (1995) conducted a similar set of measurements as part of a wider experimental study of the marine carbonate system. Again using only one single sample, they measured the variation in fCO₂ with temperature at roughly 5 °C intervals from 5.05 to 35 °C (seven data points). The reproducibility of the fCO₂ measurements was reported as ±3 µatm. The properties of the sample were reported as A_T=2320.8 ± 1.6 µmol kg⁻¹, T_C=1966.8 ± 1.2 µmol kg⁻¹ and salinity=34.95. The experimental set-up was a little different; for example, Takahashi et al. (1993) equilibrated the sample with an air headspace that was subsampled for separate measurement, whereas Lee and Millero (1995) continuously circulated the air headspace through the CO₂ analyser. Lee and Millero (1995) fitted a linear equation of the form of Eq. (5) to their data, finding a best-fitting υ_l of 4.26 % °C⁻¹. Refitting their dataset but minimising residuals in fCO₂ rather than ln (fCO₂) returns $b_{\mathrm{l}}={\mathit{\upsilon }}_{\mathrm{l}}=\mathrm{4.23}$ % °C⁻¹.

2.1.3 The van 't Hoff equation approach

Takahashi et al. (1993) did not give a theoretical basis for either of the forms (linear and quadratic; Eqs. 5 and 6) that they fitted to their dataset, nor did they give any reason to choose one over the other. However, as will be discussed later, the form chosen has an important effect on adjustments, especially over greater Δt (Sect. 3.1), and on uncertainty propagation (Sect. 3.2). Here, we provide the theoretical basis for an alternative form based on the equations governing the marine carbonate system and the van 't Hoff equation, which is derived from the well-known thermodynamic equation for the Gibbs free energy of chemical reactions.

Seawater fCO₂ is related to A_T and T_C by (e.g. Humphreys et al., 2018)

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

The asterisks in $K_{\mathrm{1}}^{\ast }$ and $K_{\mathrm{2}}^{\ast }$ show that they are stoichiometric equilibrium constants, based on the substance contents of the equilibrating species, rather than thermodynamic constants based on their activities; the prime symbol in $K_{\mathrm{0}}^{\prime }$ shows that it is a hybrid equilibrium constant, based on the substance content of CO₂(aq) and the fugacity of CO₂(g) (Zeebe and Wolf-Gladrow, 2001). Following Humphreys et al. (2018), A_T and T_C can be approximated with A_x and T_x, which include only the bicarbonate (${{\mathrm{HCO}}_{\mathrm{3}}}^{-}$) and carbonate (${{\mathrm{CO}}_{\mathrm{3}}}^{\mathrm{2}-}$) terms:

$$\begin{array}{}\text{(12)}& {\displaystyle }{A}_{\mathrm{T}}\approx A_x=\left[{{\mathrm{HCO}}_{\mathrm{3}}}^{-}\right]+\mathrm{2}\left[{{\mathrm{CO}}_{\mathrm{3}}}^{\mathrm{2}-}\right],\text{(13)}& {\displaystyle }{T}_{\mathrm{C}}\approx T_x=\left[{{\mathrm{HCO}}_{\mathrm{3}}}^{-}\right]+\left[{{\mathrm{CO}}_{\mathrm{3}}}^{\mathrm{2}-}\right].\end{array}$$

In typical seawater, the approximations in Eqs. (12) and (13) are more than 99 % accurate for T_C and more than 97 % accurate for A_T. These approximations have previously been used to understand how fCO₂ responds to changes in A_T and T_C under constant temperature (Humphreys et al., 2018). Equations (12) and (13) can be rearranged for [${{\mathrm{HCO}}_{\mathrm{3}}}^{-}$] and [${{\mathrm{CO}}_{\mathrm{3}}}^{\mathrm{2}-}$]:

$$\begin{array}{}\text{(14)}& {\displaystyle }\left[{{\mathrm{HCO}}_{\mathrm{3}}}^{-}\right]=\mathrm{2}{T}_x-A_x,\text{(15)}& {\displaystyle }\left[{{\mathrm{CO}}_{\mathrm{3}}}^{\mathrm{2}-}\right]=A_x-T_x.\end{array}$$

Under this approximation, both [${{\mathrm{HCO}}_{\mathrm{3}}}^{-}$] and [${{\mathrm{CO}}_{\mathrm{3}}}^{\mathrm{2}-}$] are constant under constant A_x and T_x, so temperature affects fCO₂ only through the $K_{\mathrm{2}}^{\ast }/K_{\mathrm{0}}^{\prime }{K}_{\mathrm{1}}^{\ast }$ term in Eq. (11). To relate this to υ, we need to quantify this sensitivity to temperature:

$$\begin{array}{}\text{(16)}& {\displaystyle }\begin{array}{rl}{\displaystyle \frac{\mathrm{d}\ln \left(f{\mathrm{CO}}_{\mathrm{2}}\right)}{\mathrm{d}t}}& {\displaystyle }\approx {\displaystyle \frac{\mathrm{d}}{\mathrm{d}t}}\ln \left({\displaystyle \frac{K_{\mathrm{2}}^{\ast }}{K_{\mathrm{0}}^{\prime }{K}_{\mathrm{1}}^{\ast }}}\right)\\ {\displaystyle }& {\displaystyle }={\displaystyle \frac{\mathrm{d}\ln K_{\mathrm{2}}^{\ast }}{\mathrm{d}t}}-{\displaystyle \frac{\mathrm{d}\ln K_{\mathrm{0}}^{\prime }}{\mathrm{d}t}}-{\displaystyle \frac{\mathrm{d}\ln K_{\mathrm{1}}^{\ast }}{\mathrm{d}t}}.\end{array}\end{array}$$

The van 't Hoff equation provides a theoretical basis for the first-order response of equilibrium constants to temperature:

$$\begin{array}{}\text{(17)}& {\displaystyle \frac{\mathrm{d}\ln K}{\mathrm{d}t}}={\displaystyle \frac{{\mathrm{\Delta }}_{\mathrm{r}}{H}^{\ominus }}{R{t}_{\mathrm{K}}^{\mathrm{2}}}},\end{array}$$

where Δ_rH^⊖ is the standard enthalpy of the reaction, R is the ideal gas constant and t_K is the temperature in kelvin. Substituting Eq. (17) into Eq. (16) gives

$$\begin{array}{}\text{(18)}& {\displaystyle \frac{\mathrm{d}\ln \left(f{\mathrm{CO}}_{\mathrm{2}}\right)}{\mathrm{d}t}}\approx {\mathit{\upsilon }}_x={\displaystyle \frac{{\mathrm{\Delta }}_{\mathrm{r}}{H}_{\mathrm{2}}^{\ominus }-{\mathrm{\Delta }}_{\mathrm{r}}{H}_{\mathrm{0}}^{\ominus }-{\mathrm{\Delta }}_{\mathrm{r}}{H}_{\mathrm{1}}^{\ominus }}{R{t}_{\mathrm{K}}^{\mathrm{2}}}}.\end{array}$$

Equation (18) shows that, if the approximations A_x and T_x are sufficiently accurate, then (to first order) ln (fCO₂) should be proportional to $\mathrm{1}/t_{\mathrm{K}}$, rather than to t.

Based on reported standard enthalpies of formation at 298.15 K (Ruscic and Bross, 2023), and assuming that these are constant under the temperature range of interest, the numerator in Eq. (18) is 25 288 ± 98 J mol⁻¹ (Appendix A). Given the approximations involved, Eq. (18) is unlikely to provide correct values for υ in seawater directly by using this value, but it does suggest a revised form that could be fitted to experimental data such as those of Takahashi et al. (1993):

$$\begin{array}{}\text{(19)}& \ln \left(f{\mathrm{CO}}_{\mathrm{2}}\right)=c_h-{\displaystyle \frac{b_h}{R{t}_{\mathrm{K}}}}.\end{array}$$

Following the same steps taken in Sect. 2.1.2 for the linear and quadratic forms of Takahashi et al. (1993), we find the following:

$$\begin{array}{}\text{(20)}& {\displaystyle }{\mathit{\upsilon }}_h={\displaystyle \frac{b_h}{R{t}_{\mathrm{K}}^{\mathrm{2}}}},\text{(21)}& {\displaystyle }{\mathrm{{\rm Y}}}_h={\displaystyle \frac{b_h}{R}}\left({\displaystyle \frac{\mathrm{1}}{t_{\text{K,0}}}}-{\displaystyle \frac{\mathrm{1}}{t_{\text{K,1}}}}\right).\end{array}$$

We fitted Eq. (19) to the Takahashi et al. (1993) dataset, optimising residuals in fCO₂ rather than ln (fCO₂) to avoid implicitly weighting the lower-fCO₂ values more heavily. The least-squares best-fit b_h was 28 995 ± 235 J mol⁻¹, with an RMSD of 2.9 µatm in fCO₂. Fitting Eq. (19) in the same way to the Lee and Millero (1995) dataset gave a least-squares best-fit b_h of 30 794 ± 357 J mol⁻¹, with an RMSD of 2.8 µatm. The different b_h values considered here are summarised in Table 1, including the parameterisation developed later in Sect. 2.4, which gives a range of possible values.

Table 1 Overview of the different values for b_h (Eqs. 19–21) that are discussed in this article.

2.2 Uncertainty in the experiment

2.2.1 Propagation equations

To quantify uncertainties in fCO₂ resulting from temperature adjustments, one needs to know the uncertainty in υ and Υ for whichever form is being used. Here, we derive uncertainty propagation equations following JCGM (2008) for the linear (Eq. 5), quadratic (Eq. 6) and van 't Hoff (Eq. 19) forms.

For the linear and van 't Hoff forms, uncertainty propagation is straightforward:

$$\begin{array}{}\text{(22)}& {\displaystyle }{\mathit{\sigma }}^{\mathrm{2}}\left({\mathit{\upsilon }}_{\mathrm{l}}\right)={\mathit{\sigma }}^{\mathrm{2}}\left(b_{\mathrm{l}}\right),\text{(23)}& {\displaystyle }{\mathit{\sigma }}^{\mathrm{2}}\left({\mathrm{{\rm Y}}}_{\mathrm{l}}\right)={\mathit{\sigma }}^{\mathrm{2}}\left(b_{\mathrm{l}}\right)\times (t_{\mathrm{1}}-t_{\mathrm{0}}{)}^{\mathrm{2}},\text{(24)}& {\displaystyle }{\mathit{\sigma }}^{\mathrm{2}}\left({\mathit{\upsilon }}_h\right)={\displaystyle \frac{{\mathit{\sigma }}^{\mathrm{2}}\left(b_h\right)}{(R{t}_{\mathrm{K}}^{\mathrm{2}}{)}^{\mathrm{2}}}},\text{(25)}& {\displaystyle }{\mathit{\sigma }}^{\mathrm{2}}\left({\mathrm{{\rm Y}}}_h\right)={\displaystyle \frac{{\mathit{\sigma }}^{\mathrm{2}}\left(b_h\right)}{R^{\mathrm{2}}}}{\left({\displaystyle \frac{\mathrm{1}}{t_{\text{K,0}}}}-{\displaystyle \frac{\mathrm{1}}{t_{\text{K,1}}}}\right)}^{\mathrm{2}},\end{array}$$

where σ² is uncertainty expressed as a variance.

For the quadratic form, covariances between the fitted coefficients must be included. In general,

$$\begin{array}{}\text{(26)}& {\mathit{\sigma }}^{\mathrm{2}}\left(p\right)=\mathbf{J}\left(p\right)\cdot {\mathbf{U}}_{\mathrm{q}}\cdot \mathbf{J}(p{)}^{\mathrm{T}},\end{array}$$

where p is the parameter of interest and the uncertainty matrix for the coefficients is

$$\begin{array}{}\text{(27)}& {\displaystyle }{\mathbf{U}}_{\mathrm{q}}=\left[\begin{array}{cc}{\mathit{\sigma }}^{\mathrm{2}}\left(a_{\mathrm{q}}\right)& \mathit{\sigma }(a_{\mathrm{q}},b_{\mathrm{q}})\\ \mathit{\sigma }(a_{\mathrm{q}},b_{\mathrm{q}})& {\mathit{\sigma }}^{\mathrm{2}}\left(b_{\mathrm{q}}\right)\end{array}\right].\end{array}$$

The Jacobian matrices for υ_q and Υ_q are

$$\begin{array}{}\text{(28)}& {\displaystyle }\mathbf{J}\left({\mathit{\upsilon }}_{\mathrm{q}}\right)=\left[\begin{array}{cc}\frac{\partial {\mathit{\upsilon }}_{\mathrm{q}}}{\partial a_{\mathrm{q}}}& \frac{\partial {\mathit{\upsilon }}_{\mathrm{q}}}{\partial b_{\mathrm{q}}}\\ \mathrm{⋮}& \mathrm{⋮}\end{array}\right]=\left[\begin{array}{cc}\mathrm{2}t& \mathrm{1}\\ \mathrm{⋮}& \mathrm{⋮}\end{array}\right],\text{(29)}& {\displaystyle }\mathbf{J}\left({\mathrm{{\rm Y}}}_{\mathrm{q}}\right)=\left[\begin{array}{cc}\frac{\partial {\mathrm{{\rm Y}}}_{\mathrm{q}}}{\partial a_{\mathrm{q}}}& \frac{\partial {\mathrm{{\rm Y}}}_{\mathrm{q}}}{\partial b_{\mathrm{q}}}\\ \mathrm{⋮}& \mathrm{⋮}\end{array}\right]=\left[\begin{array}{cc}(t_{\mathrm{1}}^{\mathrm{2}}-t_{\mathrm{0}}^{\mathrm{2}})& (t_{\mathrm{1}}-t_{\mathrm{0}})\\ \mathrm{⋮}& \mathrm{⋮}\end{array}\right].\end{array}$$

Thus, for a single measurement, Eq. (26) gives

$$\begin{array}{}\text{(30)}& {\displaystyle }{\mathit{\sigma }}^{\mathrm{2}}\left({\mathit{\upsilon }}_{\mathrm{q}}\right)=\mathrm{4}{t}^{\mathrm{2}}{\mathit{\sigma }}^{\mathrm{2}}\left(a_{\mathrm{q}}\right)+{\mathit{\sigma }}^{\mathrm{2}}\left(b_{\mathrm{q}}\right)+\mathrm{4}t\mathit{\sigma }(a_{\mathrm{q}},b_{\mathrm{q}}),\text{(31)}& {\displaystyle }\begin{array}{rl}{\displaystyle }{\mathit{\sigma }}^{\mathrm{2}}\left({\mathrm{{\rm Y}}}_{\mathrm{q}}\right)=& {\displaystyle }(t_{\mathrm{1}}^{\mathrm{2}}-t_{\mathrm{0}}^{\mathrm{2}}{)}^{\mathrm{2}}{\mathit{\sigma }}^{\mathrm{2}}(a_{\mathrm{q}})+(t_{\mathrm{1}}-t_{\mathrm{0}}{)}^{\mathrm{2}}{\mathit{\sigma }}^{\mathrm{2}}\left(b_{\mathrm{q}}\right)\\ {\displaystyle }& {\displaystyle }+\mathrm{2}(t_{\mathrm{1}}^{\mathrm{2}}-t_{\mathrm{0}}^{\mathrm{2}})(t_{\mathrm{1}}-t_{\mathrm{0}})\mathit{\sigma }(a_{\mathrm{q}},b_{\mathrm{q}}).\end{array}\end{array}$$

Finally, uncertainty in Υ can be converted into uncertainty in exp (Υ) for use in Eq. (4) with

$$\begin{array}{}\text{(32)}& {\mathit{\sigma }}^{\mathrm{2}}\left[\exp \right(\mathrm{{\rm Y}}\left)\right]=\left[\exp \right(\mathrm{{\rm Y}}){]}^{\mathrm{2}}\times {\mathit{\sigma }}^{\mathrm{2}}(\mathrm{{\rm Y}}).\end{array}$$

2.2.2 Uncertainties in fitted parameters

To use the equations in Sect. 2.2.1, the uncertainties in the fitted parameters from Eqs. (5), (6) and (19) must be known. We computed these by running Monte Carlo simulations of the Takahashi et al. (1993) experiment, as this dataset is the most widely used in practice.

We used the Takahashi et al. (1993) temperature and fCO₂ measurements as the starting point for the simulations. For each iteration of the simulation (total iterations = 10⁶), we generated a random number for each of temperature and fCO₂, which were added uniformly across the entire dataset to represent systematic errors (trueness), and then a different random number for each temperature and fCO₂ data point, to represent random errors in each measurement (precision). The distributions used for the random numbers are shown in Table 2. The simulated temperature and fCO₂ values (starting point plus systematic and random errors) were then fitted to Eqs. (5), (6) and (19) separately in each iteration. All of the forms of υ and Υ were also computed based on the fit at each iteration.

Table 2 Uncertainties in temperature and seawater fCO₂ assumed for our Monte Carlo simulations of the Takahashi et al. (1993) experiment.

^∗ Spread indicates the total range for uniform distributions or the standard deviation for normal distributions, with both distributions centred at zero.

The variances σ²(b_l) and σ²(b_h) and the variance–covariance matrix U_q (Eq. 27) were computed from the fitted coefficients across the full set of iterations.

2.2.3 Validation of uncertainty propagation

The uncertainty propagation equations in Sect. 2.2.1 have all been validated by two comparisons: (1) with the variances of the υ and Υ values computed across all of the simulations in Sect. 2.2.2 and (2) by computing the derivatives required for propagation using automatic differentiation (Maclaurin, 2016) rather than manually. The values from the equations agreed to within 0.01 % with those based on the simulations and to within 10⁻¹² % with the calculations based on automatic differentiation. The former difference is greater (albeit still negligible) because of non-linearities in the forms of υ that are not accounted for in the first-order derivatives used for uncertainty propagation here. The latter difference is attributable to rounding errors at the limit of computer precision.

2.3 PyCO2SYS extension and calculations

PyCO2SYS is a free, open-source Python package that can be used to solve the marine carbonate system, i.e. to calculate the equilibrium balance of the main acid–base systems in seawater and related properties (Humphreys et al., 2022). PyCO2SYS was originally based on the MATLAB/GNU Octave program CO2SYS.m (v2.0.5; Orr et al., 2018) – itself part of a family of similar software tools beginning with the original CO2SYS for MS-DOS (Lewis and Wallace, 1998).

2.3.1 Calculating υ in PyCO2SYS

In PyCO2SYS, υ can be calculated by finite differences or by automatic differentiation (Maclaurin, 2016). The accuracy of the former depends on the size of the difference used, whereas the latter gives the exact result directly (Humphreys et al., 2022). For consistency with the chemical buffer factors in PyCO2SYS, we used automatic differentiation to calculate υ. We constructed a function that evaluates ln (fCO₂) from A_T, T_C, t and all other relevant parameters (e.g. salinity, pressure and nutrients) and applied automatic differentiation to determine its derivative with respect to temperature. The function runs after A_T and T_C have been determined from whichever pair of known parameters the user provides. This was added to the set of results returned by PyCO2SYS as of v1.8.3 (Humphreys et al., 2024).

Previous versions of PyCO2SYS could already calculate $\partial f{\mathrm{CO}}_{\mathrm{2}}/\partial t$ using forward finite-difference derivatives as part of the uncertainty propagation tool. υ can be calculated from $\partial f{\mathrm{CO}}_{\mathrm{2}}/\partial t$ as follows:

$$\begin{array}{}\text{(33)}& \mathit{\upsilon }={\displaystyle \frac{\partial f{\mathrm{CO}}_{\mathrm{2}}}{\partial t}}(f{\mathrm{CO}}_{\mathrm{2}}{)}^{-\mathrm{1}}.\end{array}$$

We used the finite-difference $\partial f{\mathrm{CO}}_{\mathrm{2}}/\partial t$ and Eq. (33) to validate the values from automatic differentiation. Across 10⁴ random combinations of A_T, T_C, salinity, t, pressure, nutrient contents and carbonic acid dissociation constant parameterisation, we found that the two approaches were consistent, with a mean absolute difference of less than 10⁻⁴ % °C⁻¹ (or 0.01 %) and a maximum absolute difference of less than 0.1 % °C⁻¹. This check has been incorporated into the PyCO2SYS test suite to ensure continued consistency in future versions.

2.3.2 Adjusting fCO₂ to different temperatures

Prior to v1.8.3, PyCO2SYS could already convert fCO₂ (also pCO₂, [CO₂(aq)] or xCO₂) values from one temperature to another when no second carbonate system variable was provided, although using only Eq. (10), the quadratic fit of Takahashi et al. (1993). In v1.8.3, we have added the option to use the other methods discussed here instead, i.e. the linear fit of Takahashi et al. (1993) and the van 't Hoff form with several different options for b_h, including its parameterisation in Sect. 2.4.

As in previous versions of PyCO2SYS, the same adjustment is applied if pCO₂, [CO₂(aq)] or xCO₂ is given as the single known marine carbonate system parameter, as these can all be inter-converted with each other and with fCO₂ without knowledge of a second marine carbonate system parameter (Humphreys et al., 2022).

2.3.3 Components of υ

The effect of temperature on fCO₂ propagates through from all of the different equilibrium constants of the marine carbonate system. The total υ is the sum of the partial derivatives with respect to each constant:

$$\begin{array}{}\text{(34)}& \mathit{\upsilon }=\sum _K\left({\displaystyle \frac{\partial \ln \left(f{\mathrm{CO}}_{\mathrm{2}}\right)}{\partial K}}\cdot {\displaystyle \frac{\partial K}{\partial t}}\right),\end{array}$$

where K is the set of equilibrium constants and other properties that are sensitive to t and involved in calculating fCO₂ from A_T and T_C, specifically, the CO₂ solubility constant $K_{\mathrm{0}}^{\prime }$, the carbonic acid dissociation constants $K_{\mathrm{1}}^{\ast }$ and $K_{\mathrm{2}}^{\ast }$, the borate equilibrium constant $K_{\mathrm{B}}^{\ast }$, and the water dissociation constant $K_{\mathrm{w}}^{\ast }$. The dissociation constants for phosphoric acid, orthosilicic acid, hydrogen sulfide, ammonia, sulfate and hydrogen fluoride also play a role, but their contributions to υ are negligible, being much less than 1 % in the modern surface ocean, so we did not include them further in our analysis. (Their effects are still included in the total υ calculated with automatic differentiation, described in Sect. 2.3.1.)

We used the forward finite-difference derivatives calculated by PyCO2SYS (Humphreys et al., 2022) to evaluate $\partial f{\mathrm{CO}}_{\mathrm{2}}/\partial K$ and $\partial K/\partial t$ and combined these with Eq. (33) to quantify the influence of each component on υ.

2.3.4 Additional carbonic acid parameterisation

We incorporated the parameterisation of Papadimitriou et al. (2018) for the carbonic acid dissociation constants into PyCO2SYS. This parameterisation covers a range of temperature (−6 to 25 °C) and salinity (33 to 100) conditions appropriate for sea-ice brines; previously, the lowest valid temperature for any parameterisation in PyCO2SYS was −1.7 °C and the highest salinity was 50. The calculations were tested against the check values provided by Papadimitriou et al. (2018) in their Table 3, and these tests were incorporated into the PyCO2SYS automatic test suite (Humphreys et al., 2022).

2.4 Variability in υ and b_h

We used the OceanSODA-ETZH data product (Gregor and Gruber, 2021; dataset version 5, downloaded 4 December 2023) to investigate spatiotemporal variability in surface ocean υ. OceanSODA-ETZH provides fields of several marine carbonate system and other hydrographic parameters gridded across the surface ocean (1°×1°) and through time (monthly from 1985 to 2018). The parameter fields were generated using an ensemble cluster–regression approach based on observations in SOCAT (Bakker et al., 2016) and GLODAP (Lauvset et al., 2016). Other similar data products exist, but the main patterns and variability in surface ocean carbonate chemistry are well enough constrained that the choice of a particular data product will not significantly affect the global-scale and time-averaged analyses conducted here.

First, at each grid point and for each month, we computed the mean across all years of temperature, salinity, A_T and T_C. We then computed υ from these monthly mean gridded fields using PyCO2SYS v1.8.3 (Sect. 2.3), with the carbonic acid dissociation constants of Lueker et al. (2000), the bisulfate dissociation constant of Dickson (1990a), the hydrogen fluoride dissociation constant of Dickson and Riley (1979), and the total borate-to-chlorinity ratio of Uppström (1974). OceanSODA-ETZH does not include nutrient fields, so their concentrations were set to zero, but we do not expect this to affect the results because, as noted in Sect. 2.3.2, their effects on υ are negligible.

We also used the OceanSODA-ETZH data product to parameterise b_h. We calculated fCO₂ at 50 evenly spaced temperatures from −1.8 to 35.83 °C (i.e. the full range in the data product) at each monthly mean grid point. This simulates running an experiment like that of Takahashi et al. (1993) at each point. Next, we fit Eq. (19) to the generated t and fCO₂ data to find the best-fitting b_h value at each point. Finally, we collated all of the b_h values and fitted them using a function of the following form:

$$\begin{array}{}\text{(35)}& {\displaystyle }\begin{array}{rl}{\displaystyle }{b}_h=& {\displaystyle }{u}_{\mathrm{0}}+u_{\mathrm{1}}t+u_{\mathrm{2}}s+u_{\mathrm{3}}f+u_{\mathrm{4}}{t}^{\mathrm{2}}+u_{\mathrm{5}}{s}^{\mathrm{2}}\\ {\displaystyle }& {\displaystyle }+u_{\mathrm{6}}{f}^{\mathrm{2}}+u_{\mathrm{7}}ts+u_{\mathrm{8}}tf+u_{\mathrm{9}}sf,\end{array}\end{array}$$

where s is practical salinity, f is fCO₂ and u₀–u₉ are the fitted coefficients. The coefficients and their variance–covariance matrix are provided in the Supplement (Tables S1 and S2). Values of υ computed using Eq. (35) for b_h are denoted υ_p.

Fitting the parameterisation to a data product has an advantage over using regular grids of the predictors in that it guarantees including all combinations of predictors that do occur in the real ocean and avoids combinations that do not. It also causes the fit to be weighted more strongly by more prevalent combinations of conditions, leading to improved performance in global applications. However, some regions are less well represented by this global fit (Fig. S1 in the Supplement), so regionally specialised fits might work better for some studies.

A difficulty with applying our parameterisation of b_h is that in situ fCO₂ is one of the predictors. This may not always be known (e.g. if discrete measurements were conducted at a fixed temperature in the laboratory), and using a different fCO₂ will not yield an accurate b_h. To mitigate this issue, if in situ fCO₂ is not known before the conversion, we use the constant b_h fitted to the Takahashi et al. (1993) dataset (Sect. 2.1.3) to preliminarily adjust fCO₂ to the in situ temperature; we then use this preliminary fCO₂ to compute b_h with Eq. (35) before recalculating in situ fCO₂ with the computed b_h. This approach is incorporated into the PyCO2SYS implementation, with the user required to identify whether the input- or output-condition fCO₂ values represent in situ conditions for this purpose.

2.5 Open-source software

All of the analysis presented here was conducted using the free and open-source Python programming language (Python Software Foundation; https://www.python.org/, last access: 23 October 2024). The main packages used, which are also free and open source, include NumPy (Harris et al., 2020), Autograd (Maclaurin, 2016) and Xarray (Hoyer and Hamman, 2017) for data manipulation and general calculations; SciPy (Virtanen et al., 2020) for curve fitting; Matplotlib (Hunter, 2007) and Cartopy (Met Office, 2010) for data visualisation; and PyCO2SYS (Humphreys et al., 2022) for marine carbonate system calculations. Geographical features on maps were drawn using data that are freely available from Natural Earth (https://www.naturalearthdata.com/, last access: 23 October 2024).

3 Results and discussion

3.1 Functional form of υ

Our first goal is to determine whether the newly proposed form υ_h (Eq. 20) can accurately represent the temperature sensitivity of seawater fCO₂. To this end, we assess the accuracy of the approximations A_x and T_x made in developing Eq. (20) by comparing υ_h with calculations from a fully resolved marine carbonate system (Sect. 3.1.1). Next, we examine how well the different approaches to calculating υ agree with direct measurements; we do this for (1) laboratory experiments with measurements of fCO₂ at multiple temperatures for a single seawater sample (Sect. 3.1.2) and (2) across a global dataset with pairs of measurements at two different temperatures (Sect. 3.1.3).

3.1.1 Accuracy of the approximations

We made several approximations to develop Eq. (20) (Sect. 2.1.3). First, all non-carbonate components of A_T (Eq. 12) and the contribution of [CO₂(aq)] to T_C (Eq. 13) were neglected. Second, the simple integrated van 't Hoff equation (Eq. 17) was assumed to accurately represent how the various equilibrium constants vary with temperature. Third, when calculating b_h from standard enthalpies of formation (Appendix A), these were assumed to be constant across the relevant temperature range. We must therefore check whether Eq. (20) can still capture the sensitivity of fCO₂ to temperature accurately enough for our purposes before applying it to oceanographic data.

The first approximation above was also used to investigate the complementary question of how fCO₂ varies with A_T and T_C under constant temperature (Humphreys et al., 2018). There, it was observed that the gradient of isocaps (i.e. lines of constant fCO₂ through A_T–T_C phase space) was correlated with fCO₂. The approximations A_x and T_x were used as the basis for a simplified representation of the marine carbonate system from which equations were derived to explain the observed correlation. This successful application to a closely related problem provides a first indication that the same approximations may also be suitable here.

We can do a more direct test using software like PyCO2SYS, which does not make any of the approximations above when calculating how fCO₂ varies with temperature under constant A_T and T_C (Humphreys et al., 2022). Instead, the complete equations for A_T and T_C are used, and the equilibrium constants are parameterised with multiple temperature-, salinity- and pressure-sensitive terms, which implicitly include variations in the enthalpies of formation. Such calculations can therefore be used to quantify the effect of the approximations on the accuracy of Eq. (20). Equation (20) cannot be fitted directly to fCO₂ data because υ represents the derivative with respect to temperature; instead, we need to fit Eq. (19) to obtain the unknown b_h which can then be used in Eq. (20) for υ_h. If we use PyCO2SYS to calculate fCO₂ across a range of temperatures with constant A_T and T_C, denoted fCO₂(A_T,T_C), and then fit Eq. (19) to the results (Sect. 2.4), denoted fCO₂(υ_p), the residuals from the fit should reveal any reduction in accuracy due to the approximations. Large residuals would suggest that at least one approximation was unsuitable for modelling υ, whereas small residuals would mean that the approximated system does still capture the main features of the t–fCO₂ relationship that are needed to model υ.

For over 97 % of the monthly mean grid points in OceanSODA-ETZH, the RMSD of the residuals between fCO₂ calculated from Eq. (19) and from A_T and T_C is less than 1 µatm, with maximum deviations of ±2 µatm (Fig. 1). At lower $T_{\mathrm{C}}/A_{\mathrm{T}}$ ($<\sim \mathrm{0.95}$), fCO₂(υ_p) is generally greater than fCO₂(A_T,T_C) at low (<10 °C) and high (>30 °C) temperatures but less than fCO₂(A_T,T_C) at intermediate temperatures (10–30 °C), but this switches to the opposite pattern at higher $T_{\mathrm{C}}/A_{\mathrm{T}}$ (Fig. 1a). We could not find a relationship between the switchover point, which has a very low RMSD, and some property of chemical buffering in seawater, although the $T_{\mathrm{C}}/A_{\mathrm{T}}$ value at which it occurs is correlated with salinity (Fig. 1b).

Figure 1 (a) Residuals in fCO₂ for the best fit of Eq. (19) to the fCO₂ values calculated from a fully solved marine carbonate system at each monthly mean grid point in the OceanSODA-ETZH data product (Sect. 2.4), colour-coded by $T_{\mathrm{C}}/A_{\mathrm{T}}$. Thus, each point in OceanSODA-ETZH appears as a line of points across panel (a) with constant $T_{\mathrm{C}}/A_{\mathrm{T}}$. (b) Variation in the RMSD of the residuals shown in panel (a) with $T_{\mathrm{C}}/A_{\mathrm{T}}$, colour-coded by practical salinity. Thus, each point in OceanSODA-ETZH appears as a single point in panel (b), which shows the RMSD of the corresponding line of points in panel (a). The blue lines in panels (a) and (b) correspond to an RMSD of 1 µatm: points with an RMSD of less than 1 µatm fall between the blue lines in panel (a) and below the horizontal blue in panel (b). Less than 3 % of the points have an RMSD greater than 1 µatm. (c) Variation in the percentage errors in the approximations A_x (Eq. 12) and T_x (Eq. 13) with $T_{\mathrm{C}}/A_{\mathrm{T}}$ at each monthly mean grid point in OceanSODA-ETZH.

As context for interpreting this <1 µatm RMSD, the uncertainty in fCO₂ calculated from A_T and T_C that propagates from 0.1 % uncertainties in each of A_T and T_C (i.e. the GOA-ON climate quality goal; Newton et al., 2015) varies from 2 µatm at −1 °C to 9 µatm at 35 °C. The uncertainties in the equilibrium constants involved in this calculation are poorly known, but using the set proposed by Orr et al. (2018) roughly doubles the total uncertainty in calculated fCO₂ given above. The <1 µatm RMSD between fCO₂(υ_p) and fCO₂(A_T,T_C) found for 97 % of the surface ocean is thus much smaller than the uncertainties in the marine carbonate system calculations themselves. It is also very similar to the RMSD of the linear and quadratic fits to the Takahashi et al. (1993) dataset (Sect. 2.1.2) and better than the GOA-ON climate quality goal for fCO₂ measurements.

In the remaining <3 % of the global surface ocean, the RMSD can increase to over 5 µatm, which is less satisfactory but still within the GOA-ON weather quality goal (Newton et al., 2015). These higher RMSD values are constrained to $T_{\mathrm{C}}/A_{\mathrm{T}}$ values above 0.95 and to certain regions and seasons, primarily the Baltic Sea during winter and the near-coastal Arctic Ocean (Fig. S1). In these regions, the RMSD is correlated with several marine carbonate system properties, but one of the strongest correlations is with $T_{\mathrm{C}}/A_{\mathrm{T}}$ (Fig. 1b). The cause of the increased RMSD is likely inaccuracy in the approximation of T_x (Eq. 13). The fraction of T_C comprised of [CO₂(aq)], which is ignored in T_x, increases with $T_{\mathrm{C}}/A_{\mathrm{T}}$ (Fig. 1c). Consequently, the approximation that $[{{\mathrm{HCO}}_{\mathrm{3}}}^{-}{]}^{\mathrm{2}}/[{{\mathrm{CO}}_{\mathrm{3}}}^{\mathrm{2}-}]$ is constant across different temperatures (Eq. 16), which emerges from the definitions of A_x and T_x (Eqs. 12–15), becomes less accurate with increasing $T_{\mathrm{C}}/A_{\mathrm{T}}$. The fraction of non-carbonate alkalinity, which is ignored in A_x, decreases with increasing $T_{\mathrm{C}}/A_{\mathrm{T}}$, which makes A_x more accurate at higher $T_{\mathrm{C}}/A_{\mathrm{T}}$ (Fig. 1c); thus, A_x should not be the cause of the increased RMSD. The RMSD of the b_h fit is positively correlated with the difference T_C−T_x but negatively correlated with A_T−A_x (Fig. S3), supporting the hypothesis that inaccuracy in T_x rather than A_x drives the increased RMSD observed at higher $T_{\mathrm{C}}/A_{\mathrm{T}}$.

The ∼59 % of data points with $T_{\mathrm{C}}/A_{\mathrm{T}}<\mathrm{0.9}$ (in dark purple in Fig. 1a) all have the same pattern of residuals with an RMSD that is virtually independent of $T_{\mathrm{C}}/A_{\mathrm{T}}$ and less than 1 µatm. Their independence from $T_{\mathrm{C}}/A_{\mathrm{T}}$ suggests that these residuals result from the additional temperature-dependent terms beyond $\mathrm{1}/t_{\mathrm{K}}$ from the van 't Hoff equation that are included in parameterisations of the equilibrium constants (e.g. Weiss, 1974; Dickson, 1990b; Lueker et al., 2000). Including such additional terms in Eq. (19) (and thus also in Eqs. 20–21) might allow it to better capture the pattern of ln (fCO₂) variation with t, reducing the RMSD. Were υ_h to be adapted to include additional fitted parameters, the functional forms included should have a theoretical basis (e.g. Clarke and Glew, 1966) rather than being purely empirical, to ensure accurate uncertainty propagation (Sect. 3.2).

However, to include additional fitted parameters now would be premature. There are major differences between the different carbonic acid parameterisations (e.g. Fig. S2) that are not yet fully understood, and Fig. 1a and b would show different patterns if reproduced using a different parameterisation. Even though it is the parameterisation given by the best-practice guide (Dickson et al., 2007), it is still well-known that Lueker et al. (2000) does not perfectly represent the marine carbonate system to the accuracy with which we can measure it (e.g. Álvarez et al., 2020; Carter et al., 2024; Wang et al., 2023; Woosley, 2021; Woosley and Moon, 2023), and there are outstanding issues to resolve with, for example, the contribution of organic matter to total alkalinity (e.g. Ulfsbo et al., 2015; Hu, 2020; Kerr et al., 2021; Lee et al., 2024) and the borate-to-chlorinity ratio (e.g. Lee et al., 2010; Fong and Dickson, 2019). Furthermore, there are no robust uncertainty estimates for the equilibrium constants that can be propagated through to determine a meaningful uncertainty in υ_Lu00 (Orr et al., 2018). The carbonic acid and other equilibrium constant parameterisations are not optimised for this derivative, emergent property (υ); direct measurements (like those of Takahashi et al., 1993), with which more subtle variations about the first-order form in Eq. (20) can be verified, are scarce (Sect. 3.1.2).

Even with these issues remaining unresolved, we conclude that, for over 97 % of the global surface ocean, Eq. (20) with all its approximations (Sect. 2.1.3) can capture the pattern of fCO₂ sensitivity to temperature according to calculations from A_T and T_C with the carbonic acid dissociation constants of Lueker et al. (2000) to within an uncertainty of 1 µatm.

3.1.2 Direct measurements of υ

Having established how well υ_h from Eq. (20) agrees with calculations from a fully solved marine carbonate system, we next consider how these and other approaches compare with the scarce experimental evidence available.

The linear (υ_l) and quadratic (υ_q) forms of Takahashi et al. (1993) both agree well with their own measurements, so it is not clear which (if either) is more appropriate. The RMSD for υ_q (0.9 µatm) is lower than that for υ_l (1.2 µatm), but this reduction is inevitable because υ_q has an extra adjustable coefficient (Eqs. 5–8). There is a meaningful difference between υ_l and υ_q only above ∼25 °C (Fig. 2a), where υ_q is lower, meaning that fCO₂ would increase less with warming than suggested by υ_l. For example, adjusting from 25 to 35 °C would give an fCO₂ that is about 17 µatm lower with υ_q than with υ_l.

Figure 2 (a) Variation in fCO₂ with temperature, normalised to the linear fit of Takahashi et al. (1993) (υ_l=4.23 % °C⁻¹), according to the measurements of Takahashi et al. (1993) (Ta93; filled circles with vertical 1σ error bars) and the different fitted and theoretical values of υ computed under the conditions of the Takahashi et al. (1993) experiment (Sect. 2.1.2). υ_x is based on Eq. (18), using the theoretical b_h of 25 288 J mol⁻¹ (Appendix A), while ${\mathit{\upsilon }}_h^{\text{Ta93}}$ uses Eq. (19) and the best-fit b_h for the Takahashi et al. (1993) dataset of 28 995 J mol⁻¹ (Sect. 2.1.3). (b) The values of υ calculated using each of the approaches shown in panel (a) as a function of temperature.

Next, we consider the form proposed here in Eq. (19). When b_h is fixed to the value based on standard enthalpies of reaction (Sect. 2.1.3), so only the intercept (c_h) can be optimised, the resulting curve υ_x has a rather different pattern that is inconsistent with the measurements (Fig. 2a). However, also allowing b_h to be fitted to the measurements of Takahashi et al. (1993) (${\mathit{\upsilon }}_h^{\text{Ta93}}$) gives a pattern similar to υ_q, although with stronger curvature. Even though it agrees with the Takahashi et al. (1993) measurements less well than υ_l and υ_q do, with an RMSD of 2.9 µatm as opposed to ∼1 µatm, the υ_h line still falls within their uncertainty window. As for υ_q, the most significant deviation from υ_l occurs above 25 °C, where Takahashi et al. (1993) did not report any measurements. Adjusting from 25 to 35 °C would give an fCO₂ that is about 45 µatm lower with ${\mathit{\upsilon }}_h^{\text{Ta93}}$ than with υ_l.

There are 18 different parameterisations of the carbonic acid dissociation constants in PyCO2SYS v1.8.3, and each one gives a different pattern for υ. However, they are mostly negatively correlated with temperature, like υ_q and ${\mathit{\upsilon }}_h^{\text{Ta93}}$, so they all give lower fCO₂ values at higher temperatures than υ_l does (Fig. S2). Two of the parameterisations are highlighted here.

The parameterisation of Millero (1979) (υ_Mi79) agrees remarkably well with υ_x (Fig. 2). This parameterisation was designed for zero-salinity freshwater and excludes many of the non-carbonate equilibria that are active in seawater, such as that of borate.

The parameterisations that are (at least partly) based on the measurements of Mehrbach et al. (1973) fit the Takahashi et al. (1993) dataset better than those that do not, with the Lueker et al. (2000) parameterisation (υ_Lu00) having one of the best fits (RMSD=2.6 µatm) and falling within the uncertainty of the measurements (Fig. 2a). The pattern of υ_Lu00 is very similar to ${\mathit{\upsilon }}_h^{\text{Ta93}}$ (Fig. 2b). The similarity was expected, given the results discussed in Sect. 3.1.1 and the $T_{\mathrm{C}}/A_{\mathrm{T}}$ of ∼0.92 for the sample used by Takahashi et al. (1993). However, the similarity above 25 °C is more remarkable in this case, as the υ_h curves in Sect. 3.1.1 were fitted to data spanning the full temperature range from −1.8 to 35.83 °C, whereas there were no data above 24.5 °C (nor below 2.1 °C) to guide the shape of the fitted ${\mathit{\upsilon }}_h^{\text{Ta93}}$ curve here. Thus even without having a wide enough temperature range for the data to show convincing curvature in υ, the extrapolation of Eq. (20) still matches the calculations with PyCO2SYS very well, unlike υ_l and υ_q.

There are two main reasons for the differences between υ_x and ${\mathit{\upsilon }}_h^{\text{Ta93}}$ or, equivalently, between υ_Mi79 and υ_Lu00. First, ${\mathit{\upsilon }}_h^{\text{Ta93}}$ and υ_Lu00 include additional equilibria that are not present in υ_x and υ_Mi79. Some of these extra equilibria are explicitly part of the alkalinity equation in PyCO2SYS (e.g. borate); others are implicit, such as the formation of ${{\mathrm{CaCO}}_{\mathrm{3}}}^{\mathrm{0}}$ and ${{\mathrm{MgCO}}_{\mathrm{3}}}^{\mathrm{0}}$ ion pairs (Millero and Pierrot, 1998), which are incorporated within [${{\mathrm{CO}}_{\mathrm{3}}}^{\mathrm{2}-}$] and $K_{\mathrm{2}}^{\ast }$ in PyCO2SYS (Humphreys et al., 2022). The approximation that all non-carbonate species are negligible within A_T (Eq. 12) is thus more accurate for υ_x and υ_Mi79 (hence their mutual agreement) than it is for ${\mathit{\upsilon }}_h^{\text{Ta93}}$ and υ_Lu00. Second, the van 't Hoff equation with standard enthalpies of reaction predicts the temperature sensitivity of thermodynamic equilibrium constants. However, chemical equilibria in seawater are modelled using stoichiometric equilibrium constants, which have an additional temperature dependency related to the major ion composition and ionic strength of the solution (Pitzer, 1991). Both reasons would affect ${\mathit{\upsilon }}_h^{\text{Ta93}}$ and υ_Lu00 (and hence their deviation from υ_x) but not υ_Mi79, which is consistent with υ_x.

The strong agreement between υ_x and υ_Mi79 strengthens our conclusion from Sect. 3.1.1 that the approximations made in deriving Eq. (20) are valid for representing the first-order form of υ. It also means that under certain low-salinity conditions, Eq. (18) can be used directly with Δ_rH^⊖ values from the literature and no further fitting to quantify υ accurately with respect to calculations with an equilibrium model (PyCO2SYS). The strong agreement between ${\mathit{\upsilon }}_h^{\text{Ta93}}$ and υ_Lu00 suggests that the extra equilibria and ionic strength effects can be adequately accounted for by adjusting the value of b_h and without changing the form of the t–ln (fCO₂) relationship.

It is unfortunate that the dataset of Takahashi et al. (1993) has a maximum temperature of 24.5 °C, because extra data points above 25 °C could unambiguously distinguish between the different forms of υ. At lower temperatures, the differences between the forms are similar in size to the uncertainty of the measurements, so many more highly accurate measurements would need to be conducted to robustly distinguish between the forms (Fig. 2a). The Kanwisher (1960) measurements covered an even narrower temperature range from 11 to 21.5 °C, which is also insufficient to convincingly show non-linearity.

One published dataset does cover a wider temperature range, with seven measurements of a single sample from 5.05 to 35 °C (Lee and Millero, 1995). These measurements can be fitted well with both the linear and quadratic forms (Eqs. 5 and 6), with RMSD values of 1.8 µatm (linear) and 1.5 µatm (quadratic). Lee and Millero (1995) reported a constant υ_l of 4.26 % °C⁻¹, stating that the difference between this and the value of Takahashi et al. (1993) was due to a different $T_{\mathrm{C}}/A_{\mathrm{T}}$ for the sample measured (0.85 versus ∼0.92). As with the Takahashi et al. (1993) dataset, ${\mathit{\upsilon }}_h^{\text{LM95}}$ fits less well than either υ_l or υ_q, with an RMSD of 2.8 µatm, although the fitted curve again falls within the uncertainty window of the data points (±3 µatm). However, in this case, ${\mathit{\upsilon }}_h^{\text{LM95}}$ is completely different from the calculation with the Lueker et al. (2000) carbonic acid dissociation constants (Fig. 3). It also predicts significantly less deviation from υ_l or υ_q than in the Takahashi et al. (1993) case.

Figure 3 Variation in fCO₂ with temperature, normalised to the linear fit of Lee and Millero (1995) (υ_l=4.26 % °C⁻¹), according to the measurements of Lee and Millero (1995) (LM95; filled circles with vertical 1σ error bars) and the different fitted and theoretical values of υ computed under the conditions of the Lee and Millero (1995) experiment. ${\mathit{\upsilon }}_h^{\text{LM95}}$ uses Eq. (19) and the best-fit b_h for the Lee and Millero (1995) dataset of 30 794 J mol⁻¹ (Sect. 2.1.3). The line for υ_Lu00 continues to decrease after exiting the figure, reaching −40 µatm at 35 °C.

The fitted b_h value of 30 794 J mol⁻¹ for ${\mathit{\upsilon }}_h^{\text{LM95}}$ is unusual, falling outside the range fitted across the entire monthly mean OceanSODA-ETZH data product (25 139 to 30 731 J mol⁻¹; Sect. 2.4; Fig. S2a). Based on the A_T, T_C and salinity reported by Lee and Millero (1995), and with the carbonic acid dissociation constants of Lueker et al. (2000), we would have expected a b_h of 29 825 J mol⁻¹. Unlike Takahashi et al. (1993), which is not inconsistent with Eq. (20) but lacks data at high enough temperatures to convincingly support or reject it, the Lee and Millero (1995) dataset is clearly inconsistent with Eq. (20).

However, the Lee and Millero (1995) measurements are also inconsistent in a number of other ways not related to Eq. (20). We refitted the Lee and Millero (1995) dataset with the linear form (Eq. 5) and found a lower RMSD (1.6 µatm) using the υ_l value of 4.23 % °C⁻¹ from Takahashi et al. (1993). It appears that Lee and Millero (1995) minimised residuals in ln (fCO₂) in their fit rather than minimising residuals in fCO₂; the latter approach returns almost exactly the Takahashi et al. (1993) value of 4.23 % °C⁻¹. Minimising residuals in ln (fCO₂) implicitly assumes that uncertainty in fCO₂ is a constant fraction of the fCO₂ value rather than an absolute fixed amount, but Lee and Millero (1995) reported the “reproducibility” as ±3 µatm, not as a percentage. While, at first glance, it may seem like an improvement that fitting the data consistently with the uncertainty brings the υ_l value in line with Takahashi et al. (1993), this in fact contradicts the reasoning of Lee and Millero (1995) themselves in that a different value is expected because of the rather different $T_{\mathrm{C}}/A_{\mathrm{T}}$. For $T_{\mathrm{C}}/A_{\mathrm{T}}$ to be the same as that of Takahashi et al. (1993), either T_C or A_T would need to be adjusted by over 160 µmol kg⁻¹, which is far more than permitted by the measurement accuracy.

The Lee and Millero (1995) dataset is also not consistent with PyCO2SYS calculations. As mentioned above, the carbonic acid dissociation constants of Lueker et al. (2000) do not fit well at all, with an RMSD of 20 µatm (Fig. 3). While Lee and Millero (1995) claim that their data are “thermodynamically consistent” with the constants of Roy et al. (1993) and Goyet and Poisson (1989), neither of these fits very well, with both having an RMSD of ∼6 µatm (Fig. 3). However, they did acknowledge that more studies of the temperature sensitivity of the marine carbonate system are needed, pointing to possible problems with the Weiss (1974) formulation for CO₂ solubility, which is still used in PyCO2SYS (Humphreys et al., 2022).

In summary, while comparisons with the experimental datasets do not rule out Eq. (20), they are less convincing than the comparisons with PyCO2SYS calculations in Sect. 3.1.1. For the Takahashi et al. (1993) dataset, it is plausible that extra measurements at higher temperatures could have revealed significant deviations from linearity, as predicted independently by both PyCO2SYS and Eq. (20), although the linear υ_l still fits the data more closely than υ_h does. The Lee and Millero (1995) dataset, which does include higher-temperature measurements, does not support this. This latter dataset does feature a number of other inconsistencies that could cast doubt on its accuracy, but, alternatively, there may be some aspect of the equilibrium chemistry of CO₂ in seawater that is missing from or incorrectly parameterised within existing models such as PyCO2SYS. It is not possible to draw a reliable conclusion either way from only these two small datasets. We need a much bigger body of measurements to be able to robustly distinguish between the different possible forms of υ.

3.1.3 Global data compilation

The Wanninkhof et al. (2022a) compilation provides much more data (∼2000 data points) with which to test the different temperature adjustment approaches, including across variations in salinity and carbonate chemistry (A_T and T_C). The compilation consists of fCO₂ measured in discrete seawater samples at 20 °C (irrespective of the in situ temperature), along with fCO₂ measured at (or close to) the in situ temperature with a continuously flowing underway system at closely corresponding times and locations. Independent measurements of A_T and T_C are also available for many of the data points.

This dataset goes some way towards fixing the data scarcity problem mentioned in Sect. 3.1.2, but it is still not a complete substitute for additional experiments like that of Takahashi et al. (1993). In particular, each sample has measurements of fCO₂ at only two different temperatures, so the form of the t–fCO₂ relationship between those two points cannot be determined for any individual sample. Instead, we must infer the suitability of different forms of υ from how well they work across the entire dataset.

We adjusted the fCO₂ at 20 °C from the discrete samples to the temperature of the in situ measurement using five different approaches: the linear (υ_l) and quadratic (υ_q) equations of Takahashi et al. (1993) (Sect. 2.1.2); the van 't Hoff form with b_h constant as fitted to the Takahashi et al. (1993) dataset (${\mathit{\upsilon }}_h^{\text{Ta93}}$; Sect. 2.1.3) and parameterised as a function of temperature, salinity and fCO₂ (υ_p; Sect. 2.4); and calculation with PyCO2SYS using T_C as the second known parameter and the carbonic acid dissociation constants of Lueker et al. (2000) (υ_Lu00). To compare the accuracy of these adjustments, we computed the rolling mean of the differences between fCO₂ adjusted from 20 °C to the in situ temperature and fCO₂ measured directly at the in situ temperature for each approach, using a ±4 °C window (Fig. 4). Like Wanninkhof et al. (2022a), we excluded outliers falling more than ±2 standard deviations from zero, which was 1.4 % of the total data points.

Figure 4 Difference between fCO₂ adjusted from 20 °C to the in situ temperature and a direct measurement at (or close to) the in situ temperature, shown as a rolling mean (±4 °C window) for each adjustment approach, for the data in the Wanninkhof et al. (2022a) compilation (Wa22). The shaded area shows the 2σ uncertainty in the rolling mean, centred about zero, which is a function of the number of data points within each interval and virtually identical for every approach.

For negative Δt there are relatively few data (about 14 % of the dataset), such that the different rolling mean lines overlap and mostly fall within the large uncertainty window (Fig. 4). It is not possible to meaningfully distinguish between the different approaches. All appear to overestimate fCO₂(t₁) for Δt from about −14 to −5 °C, by a similar amount to each other, with some falling just outside the uncertainty window of the rolling mean. The consistency between the different approaches here is expected, as they are most similar to each other in the corresponding t₁ range from about 10 to 20 °C (Fig. 2a). The minor deviations for negative Δt are therefore probably due to the sparsity of data in this temperature range, which means that errors due to measurement uncertainty and imperfect co-location of samples are not averaged out. We would expect these deviations from zero to diminish were more data added to the compilation.

The Δt range above 0 °C contains the majority of the dataset (86 %) and sufficient data for the approaches to diverge clearly from each other and for some to move well outside the uncertainty window (Fig. 4). Here, the linear form (υ_l) significantly overestimates fCO₂(t₁) as Δt increases. Wanninkhof et al. (2022a) tried fitting υ_l directly to this data compilation to improve the accuracy of the adjustment but were not able to fix this issue. Instead, they concluded that calculations with a fully resolved carbonate system were still required for Δt greater than about 1 °C, recommending the Lueker et al. (2000) carbonic acid dissociation constants with the Uppström (1974) borate-to-chlorinity ratio for the best agreement (υ_Lu00). Like υ_l, the quadratic form (υ_q) overestimates fCO₂(t₁) for positive Δt. The error is about half of that for υ_l but still well outside the uncertainty window. The van 't Hoff form with constant b_h (${\mathit{\upsilon }}_h^{\text{Ta93}}$) has a similar error to υ_q but with the opposite sign; i.e. fCO₂(t₁) is underestimated. However, the van 't Hoff form with a variable, parameterised b_h (υ_p) gives almost perfect agreement, even for the largest Δt, performing even better than the calculation with a fully resolved carbonate system (υ_Lu00), which also falls just within the uncertainty window.

This analysis shows that the van 't Hoff form is at least as capable as the other approaches of adjusting fCO₂ to different temperatures across this data compilation. When variability in b_h is accounted for by the parameterisation proposed in Sect. 2.4, which is completely independent of this data compilation, this approach (υ_p) can be used for accurate Δt adjustments of up to at least 10 °C. More data for negative Δt are still essential to improve our confidence in this conclusion across the full temperature range found in the ocean.

Part of the motivation for the new approach to υ developed here (Sect. 2.1.3) was to avoid needing to know a second marine carbonate system parameter when adjusting fCO₂ to different temperatures. υ_p with fixed b_h does not meet this requirement, because, as for υ_l and υ_q, the adjusted fCO₂ values are not accurate enough. Initially it appears that υ_p does meet the requirement – its only inputs are fCO₂, temperature and salinity (Eq. 35), and it is equally or more accurate than calculations with PyCO2SYS (Fig. 4). However, arguably, the need for a second parameter has been sidestepped rather than truly eliminated. The parameterisation for b_h works because, to first order, b_h, υ_h and fCO₂ are closely correlated with the ratio $T_{\mathrm{C}}/A_{\mathrm{T}}$ but less sensitive to the absolute values of either T_C or A_T (Goyet et al., 1993; Wanninkhof et al., 2022a). In other words, the fCO₂ term in the parameterisation acts as a proxy for $T_{\mathrm{C}}/A_{\mathrm{T}}$, implicitly providing the information about T_C and A_T necessary to calculate υ.

It would be possible to parameterise υ_l and υ_q across the ocean, as we have done for υ_h (i.e. υ_p), which would improve their agreement with the data compilation too. However, we have already established that there is no theoretical reason to expect the t–ln (fCO₂) relationship to be linear, especially at higher t, and Wanninkhof et al. (2022a) tried fitting υ_l to this dataset directly without satisfactory results. Further benefits of υ_h over υ_q are discussed in Sect. 3.2 on uncertainty propagation.

3.2 Experimental uncertainty

We performed Monte Carlo simulations of the Takahashi et al. (1993) experiment (Sect. 2.2.2) to determine the uncertainty in υ computed from that dataset using the linear (υ_l; Eq. 5), quadratic (υ_q; Eq. 6) and van 't Hoff (${\mathit{\upsilon }}_h^{\text{Ta93}}$; Eq. 19) fits.

Were the linear fit a valid representation of the t–ln (fCO₂) relationship, then the 1σ uncertainty in experimentally determined υ_l would be 0.035 % °C⁻¹. The biggest contributor (0.028 % °C⁻¹ on its own) is the precision of the fCO₂ measurement, with 0.020 % °C⁻¹ from systematic bias in fCO₂ alone. Systematic biases can affect υ because it is the gradient of the natural logarithm of fCO₂ (Eq. 1), which is a non-linear transformation. Temperature uncertainties have a much smaller impact, with 0.001 % °C⁻¹ uncertainty in υ arising from t precision alone and no effect from t bias. Propagating through to exp (Υ_l) for a Δt of +1 °C gives an uncertainty of just under 0.04 % (Figs. 5 and S4). Following Eq. (4), the percentage uncertainty in exp (Υ) propagates through to the same percentage uncertainty in the adjusted fCO₂ value at t₁, so this calculation is most relevant for understanding how uncertainties in υ affect temperature-adjusted fCO₂ values.

Figure 5 The 1σ uncertainty in exp (Υ) for $\mathrm{\Delta }t=+\mathrm{1}$ °C, which is equal to the percentage uncertainty in fCO₂ at t₁ caused by uncertainty in υ, due to experimental uncertainties in t and fCO₂ in the Takahashi et al. (1993) dataset only and fitted with the linear (dotted line; Eq. 5) and quadratic (dashed line; Eq. 6) forms as well as the fitted van 't Hoff form (solid line; Eq. 19). The shaded area shows the range of t from the Takahashi et al. (1993) experiment, while the full t-axis range matches OceanSODA-ETZH. The uncertainty in υ itself as a function of temperature has a virtually identical pattern to that shown here for exp (Υ) (e.g. an uncertainty in exp (Υ) of 0.04 %, with $\mathrm{\Delta }t=+\mathrm{1}$ °C, corresponds to an uncertainty in υ of about 0.04 % °C⁻¹; Fig. S4).

With the quadratic fit (Eq. 6), the 1σ uncertainty in the slope (b_q) is 0.127 % °C⁻¹, whereas the corresponding value in the squared term (a_q) is $\mathrm{41.2}\times {\mathrm{10}}^{-\mathrm{6}}$ °C⁻², and there is a covariance between these of $\mathit{\sigma }(a_{\mathrm{q}},b_{\mathrm{q}})=-\mathrm{51}\times {\mathrm{10}}^{-\mathrm{9}}$ °C⁻³. Propagating through to υ_q and exp (Υ_q) gives an uncertainty that is similar to that for the linear fit in the middle of the t range but is up to 4 times greater in colder and warmer waters (Figs. 5 and S4). It is essential to include the covariance term for propagation; failure to do so wrongly inflates the final uncertainty.

The van 't Hoff fit has a 1σ uncertainty in the b_h of 216.4 J mol⁻¹. Propagating through to ${\mathit{\upsilon }}_h^{\text{Ta93}}$ and $\exp \left({\mathrm{{\rm Y}}}_h^{\text{Ta93}}\right)$ gives a similar uncertainty to that for υ_l, decreasing slightly with increasing temperature (Figs. 5 and S4).

The uncertainties calculated here, for $\mathrm{\Delta }t=+\mathrm{1}$ °C, are all lower than the climate quality target uncertainty for seawater fCO₂ measurements of 0.5 % proposed by GOA-ON (Newton et al., 2015). This means that this uncertainty component alone would not cause an fCO₂(t₁) value with Δt up to ±1 °C to breach this target, regardless of which form was correct, although it is probably not small enough to ignore in a complete uncertainty budget. Also, its contribution scales roughly in proportion with Δt; therefore, for greater Δt, this component would become a greater part of the overall uncertainty in fCO₂(t₁).

Why do the different approaches have different uncertainties and how do we choose which approach and corresponding uncertainty to use? The linear and van 't Hoff forms both appear to have similarly low uncertainties, but the temperature adjustments that they predict differ from each other by more than the apparent uncertainty windows, especially above 25 °C (Fig. 2). The quadratic form has a completely different uncertainty profile with apparently much greater uncertainty at the edges of the temperature range (Fig. 5). The assumption implicit within any uncertainty calculation – that the model being fitted to the data is valid and that all of its adjustable parameters are truly unknown and necessary – is key to resolving this.

As we have shown in Sect. 3.1, it is not consistent with existing models of the marine carbonate system for υ to be constant and independent from temperature, so the uncertainty calculated for υ_l is valid only under this premise. In other words, the assumptions that “υ is constant” and “υ varies with temperature” are mutually incompatible. This is how υ_l and ${\mathit{\upsilon }}_h^{\text{Ta93}}$ can deviate from each other by more than their apparent uncertainties.

It makes sense that υ_q has greater uncertainty than υ_l and ${\mathit{\upsilon }}_h^{\text{Ta93}}$, as υ_q has an extra coefficient to be constrained by the same number of data points. Given more data at higher t, υ_q could plausibly be fitted to a shape almost identical to υ_h, yet it would still appear to have greater uncertainty for this reason. When using υ_q, we accept that there is some curvature in the t–ln (fCO₂) relationship, but we implicitly assume that we do not know what form that curvature takes, other than it can be represented by some second-order polynomial. The adjustable parameters that represent the curvature are thus constrained empirically only by the available data, leading to much higher apparent uncertainty towards its edges. However, we now know that, to first order, υ should follow a particular curvature that can be represented with only one adjustable parameter (b_h); it is not an unknown that must be found empirically. The lower uncertainty in ${\mathit{\upsilon }}_h^{\text{Ta93}}$ and $\exp \left({\mathrm{{\rm Y}}}_h^{\text{Ta93}}\right)$ relative to υ_q and exp (Υ_q) is therefore more realistic (Fig. 5), with the caveat that more experimental work is still needed to determine whether the expected relationship from Eq. (20) can be produced in laboratory measurements (Sect. 3.1.2).

While small, this uncertainty component is probably not negligible for the most accurate analyses. More independent reproductions of the Takahashi et al. (1993) experiment, ideally including temperatures up to 35 °C, could reduce this component into a negligible term. However, this component reflects only the experimental uncertainties in fCO₂ and temperature. To use this value as the total uncertainty in υ and/or exp (Υ) would be to assume that the single sample measured by Takahashi et al. (1993) is representative of the global ocean. It does not account for spatiotemporal variability in υ, which we discuss in Sect. 3.3.

3.3 Spatial and temporal variability

3.3.1 Variability in υ_Lu00

The time-averaged theoretical value of surface ocean υ_Lu00 in OceanSODA-ETZH (Gregor and Gruber, 2021) follows a latitudinal distribution, with lower values near the Equator (Fig. 6a). The range in υ_Lu00 is more than 15 % of its mean, from less than 4.0 % °C⁻¹ to over 4.6 % °C⁻¹. The pattern is dominantly driven by temperature. The other variables (A_T, T_C and salinity) together act in the opposite sense, slightly reducing the latitudinal gradient in υ_Lu00. This can be seen from the distribution of ${\mathit{\upsilon }}_h^{\text{Ta93}}$, which uses the constant b_h found by fitting to the Takahashi et al. (1993) dataset (Sect. 2.1.3) and has a similar but stronger latitudinal gradient (Fig. 6b). A counterintuitive consequence of this is that, were the Takahashi et al. (1993) experiment to be repeated systematically across the globe, fitting with a linear equation at each location, this would return the opposite pattern of υ_l to what we see in Fig. 6a: higher values near the Equator and lower near the poles. The variability in υ_l as determined following Takahashi et al. (1993) is independent from the original in situ temperature of the sample, so the experiment would reveal variability in υ only due to the non-thermal parameters.

Figure 6 Spatial and temporal variability in υ in OceanSODA-ETZH using the Lueker et al. (2000) parameterisation of the carbonic acid dissociation constants. (a) Theoretical υ_Lu00 calculated from A_T, T_C and auxiliary variables, averaged through time across the entire OceanSODA-ETZH dataset. (b) ${\mathit{\upsilon }}_h^{\text{Ta93}}$ calculated from OceanSODA-ETZH temperature using the b_h value fitted to the Takahashi et al. (1993) measurements, again averaged through time. Panels (a) and (b) use the same colour scale, and the white contour in both panels shows the Takahashi et al. (1993) linear fit value of υ_l=4.23 % °C⁻¹. (c) Mean seasonal range and (d) multi-decadal trend in υ_Lu00.

The seasonal range in υ is a similar size to, although slightly less than, the spatial range in its time-averaged value. It is less than 0.2 % °C⁻¹ across most of the open ocean, reaching a maximum of just over 0.4 % °C⁻¹, or around 10 % of the mean, in some near-coastal areas with strong seasonal variability in temperature and/or $T_{\mathrm{C}}/A_{\mathrm{T}}$ (Fig. 6c).

Together, the spatial and seasonal variability in υ_Lu00 has a standard deviation of 0.23 % °C⁻¹. This could be considered to represent the 1σ uncertainty in using a constant υ_l globally, keeping in mind that it has strong regional systematic biases rather than being randomly distributed. It is significantly greater than the apparent experimental uncertainty in υ_l or ${\mathit{\upsilon }}_h^{\text{Ta93}}$ of around 0.04 % °C⁻¹ (Sect. 3.2). Experiments like that of Takahashi et al. (1993) should be sufficiently accurate to detect the variability in υ if conducted at different locations around the global ocean. However, this also means that spatiotemporal variability that is unaccounted for would lead to uncertainty in υ and exp (Υ) that is several times greater in practice than that computed in Sect. 3.2.

On multi-decadal timescales, υ has been decreasing globally at an average rate of 0.001 $\mathrm{\%}\mathrm{°}{\mathrm{C}}^{-\mathrm{1}}{\mathrm{yr}}^{-\mathrm{1}}$ (Fig. 6d). The decrease has been primarily driven by increasing $T_{\mathrm{C}}/A_{\mathrm{T}}$ with a smaller contribution from warming of the surface ocean. This trend is consistent with that identified by Wanninkhof et al. (2022a), and it should lead to a minor decrease in the temperature dependency of surface ocean fCO₂ while atmospheric CO₂ levels continue to rise.

Wanninkhof et al. (2022a) pointed out that υ is greater at lower temperatures and lower at lower $T_{\mathrm{C}}/A_{\mathrm{T}}$ and that temperature and $T_{\mathrm{C}}/A_{\mathrm{T}}$ are, to first order, negatively correlated in the global surface ocean (largely due to the influence of temperature on CO₂ solubility), leading to them having opposing effects on υ and thus reducing its variability. We find that this effect does reduce spatial variability in υ, although by 30 % at most relative to a constant-$T_{\mathrm{C}}/A_{\mathrm{T}}$ scenario. For a temperature change from −1.7 to 30 °C, the main body of the OceanSODA-ETZH monthly mean data product has a range in υ_Lu00 from about 3.97 % to 4.60 % °C⁻¹, whereas at a constant $T_{\mathrm{C}}/A_{\mathrm{T}}$ of 0.9, υ_Lu00 varies from 3.85 % to 4.72 % °C⁻¹ (Fig. 7). However, this partial cancellation of the effect of temperature on υ by $T_{\mathrm{C}}/A_{\mathrm{T}}$ is not relevant for making temperature adjustments; it does not mean that the temperature sensitivity of υ can be ignored or is any less important. This is because adjustments are made on single samples under constant $T_{\mathrm{C}}/A_{\mathrm{T}}$. The sensitivity of υ to temperature for any individual sample is not lessened by the negative correlation between temperature and $T_{\mathrm{C}}/A_{\mathrm{T}}$ observed across the global dataset. In other words, coupled temperature–$T_{\mathrm{C}}/A_{\mathrm{T}}$ variations in the surface ocean reduce the spatiotemporal variability in υ under in situ conditions, but they do not affect by how much υ changes with temperature for a particular sample. Consequently, the apparent reduced variability in υ is relevant only for small (e.g. Δt<1 °C) temperature adjustments.

Figure 7 Each coloured line shows how υ_Lu00 varies with temperature at a constant $T_{\mathrm{C}}/A_{\mathrm{T}}$, with T_C always at the median value of the monthly mean OceanSODA-ETZH data product (i.e. 2045 µmol kg⁻¹). $T_{\mathrm{C}}/A_{\mathrm{T}}$ increases from 0.6 (dark purple, top right) to 1.2 (yellow, bottom left). The light-blue points in the background show all monthly mean data from OceanSODA-ETZH.

3.3.2 Components of υ

Variability in υ through space and time arises primarily from the temperature dependence of the equilibrium constants of the marine carbonate system. Across the contemporary global surface ocean, the biggest contributor to υ is $K_{\mathrm{0}}^{\prime }$, the solubility constant for CO₂, which typically represents 58 %–89 % of total υ (Fig. 8). The carbonic acid equilibrium constants $K_{\mathrm{1}}^{\ast }$ and $K_{\mathrm{2}}^{\ast }$ contribute a similar size effect on υ to each other and slightly less than $K_{\mathrm{0}}^{\prime }$. However, the $K_{\mathrm{1}}^{\ast }$ contribution is always negative, whereas the $K_{\mathrm{2}}^{\ast }$ contribution is always positive; therefore, they partly cancel each other out. Smaller contributions are made by the borate (9 %–21 %) and water (1 %–8 %) equilibria ($K_{\mathrm{B}}^{\ast }$ and $K_{\mathrm{w}}^{\ast }$, respectively). All other equilibrium constants have a negligible influence on υ in the contemporary ocean, although they might become important at pH values closer to their pK.

Figure 8 Violin plot of the contributions of each equilibrium system to υ as distributed across the monthly averaged OceanSODA-ETZH data product. From the left, the first three violins show the contributions of the CO₂ solubility constant ($K_{\mathrm{0}}^{\prime }$) and the borate ($K_{\mathrm{B}}^{\ast }$) and water ($K_{\mathrm{w}}^{\ast }$) equilibria. The remaining columns show the contributions of the carbonic acid dissociation constants $K_{\mathrm{1}}^{\ast }$ and $K_{\mathrm{2}}^{\ast }$ as well as the total υ for each parameterisation in PyCO2SYS v1.8.3. The horizontal dashed line through the “Total” group is at the Takahashi et al. (1993) linear fit value (υ_l=4.23 % °C⁻¹). The colours are to help the eye track the groups vertically; they do not have any further meaning. The lower horizontal axis shows the option code within PyCO2SYS for the carbonic acid constants, while the upper axis shows an abbreviation of the corresponding study (presented as “option code/abbreviation” in the following): 1/Ro93 – Roy et al. (1993); 2/GP89 – Goyet and Poisson (1989); 3/DM87-H, 4/DM87-M and 5/DM87-HM – refits of various combinations of Hansson (1973a, b) and Mehrbach et al. (1973) by Dickson and Millero (1987); 6/Me73 and 7/Me73-P – GEOSECS fits of Mehrbach et al. (1973); 8/Mi79 – Millero (1979), for freshwater; 9/CW98 – Cai and Wang (1998); 10/Lu00 – Lueker et al. (2000); 11/MM02 – Mojica Prieto and Millero (2002); 12/Mi02 – Millero et al. (2002); 13/Mi06 – Millero et al. (2006); 14/Mi10 – Millero (2010); 15/Wa13 – Waters and Millero (2013); 16/Su20 – Sulpis et al. (2020); 17/SB21 – Schockman and Byrne (2021); 18/Pa18 – Papadimitriou et al. (2018).

Each parameterisation of the carbonic acid dissociation constants gives a different υ distribution, but there are some similarities (Fig. 8). We could not identify a link between the functional forms of the parameterisations (i.e. which combinations of t and salinity terms were used) and the distributions shown in Fig. 8; rather, the type of seawater, specific dataset(s) and methodology used to determine the constants were the main controls.

In Fig. 8, the two options labelled “GEOSECS” (6 and 7; Me73 and Me73-P) both use the measurements of Mehrbach et al. (1973), as parameterised by Peng et al. (1987) (because Mehrbach et al., 1973, did not provide parameterised functions for their measurements); curiously, total υ with these options has a roughly opposite spatial distribution to all of the other options, with the lowest υ found at high latitudes. Although the GEOSECS options also use different $K_{\mathrm{B}}^{\ast }$ and $K_{\mathrm{w}}^{\ast }$ parameterisations from all of the other options (Humphreys et al., 2022), which modulates the pattern, the inversion is caused by $K_{\mathrm{1}}^{\ast }$ and $K_{\mathrm{2}}^{\ast }$. In any case, these options are present in (Py)CO2SYS to reproduce calculations from the GEOSECS era and are not intended for use with modern measurements. Aside from this pair, all other parameterisations based at least in part on the Mehrbach et al. (1973) dataset (labelled “Mehrbach”) have similar distributions for $K_{\mathrm{1}}^{\ast }$, $K_{\mathrm{2}}^{\ast }$ and total υ. The SB21 parameterisation (Schockman and Byrne, 2021) consists of new, spectrophotometric measurements of the product $K_{\mathrm{1}}^{\ast }{K}_{\mathrm{2}}^{\ast }$ which were paired with the Waters and Millero (2013) $K_{\mathrm{1}}^{\ast }$ to find $K_{\mathrm{2}}^{\ast }$, resulting in virtually zero variability in total υ overall. This low variability in total υ is echoed by the Su20 parameterisation (Sulpis et al., 2020), which is based on field observations during which A_T, T_C, pH and fCO₂ were measured simultaneously, but the low variability is arrived at in a different way, with rather different distributions for the individual $K_{\mathrm{1}}^{\ast }$ and $K_{\mathrm{2}}^{\ast }$ effects. However, this does not mean that these parameterisations support the use of a constant, temperature-insensitive υ, rather just that the effects of temperature and seawater chemistry (e.g. $A_{\mathrm{T}}/T_{\mathrm{C}}$) on in situ υ cancel each other out across the surface ocean more completely than they do for Lueker et al. (2000) (Fig. 7). Indeed, neither SB21 nor Su20 can be fitted to the Takahashi et al. (1993) dataset very well (Fig. S2). Parameterisations Pa18 (Papadimitriou et al., 2018), focused on high-salinity, low-temperature sea-ice brines, and Mi02 (Millero et al., 2002), based on field measurements, both have similar variability in υ to the “Mehrbach” group.

There is a noticeable difference in the distributions, especially for $K_{\mathrm{1}}^{\ast }$, when moving to the group of parameterisations based on synthetic seawater (Dickson and Millero, 1987; Goyet and Poisson, 1989; Roy et al., 1993; DM87-H and DM87-HM, GP89, and Ro93, respectively). This suggests that some components of natural seawater that were not replicated in the synthetic seawater used in their experiments may affect the apparent temperature sensitivity of the equilibrium constants being measured and thus υ. For example, the synthetic seawater recipes used by Goyet and Poisson (1989) and Roy et al. (1993) both did not include any borate, unlike real seawater (Lee et al., 2010; Uppström, 1974). Suspected interactions between borate and carbonate species (Mojica Prieto and Millero, 2002) were not explicitly modelled when evaluating the equilibrium constants, so their effects would be implicitly included within the $K_{\mathrm{1}}^{\ast }$ and $K_{\mathrm{2}}^{\ast }$ values in real seawater but not in the synthetic seawater. Dissolved organic matter could also have an equivalent effect. Ultimately, any dissolved species that reacts in an equilibrium with water (i.e. any weak acid or base) or with any component of T_C (i.e. CO₂(aq), ${{\mathrm{HCO}}_{\mathrm{3}}}^{-}$ or ${{\mathrm{CO}}_{\mathrm{3}}}^{\mathrm{2}-}$) could affect υ.

The parameterisation of Mi79 is also shown for completeness, although its distribution is not meaningful because it is intended for zero-salinity freshwater (Millero, 1979).

The parameterisations of $K_{\mathrm{1}}^{\ast }$ and $K_{\mathrm{2}}^{\ast }$ differ from each other sufficiently that, given the high accuracy with which we can measure fCO₂ and thus experimentally determine υ (Sect. 3.2), more experiments like that of Takahashi et al. (1993) but including a wider range of temperatures (i.e. at least up to 35 °C) would be a powerful tool to help distinguish between the different parameterisations of the equilibrium constants (Fig. S2) and other options (e.g. borate-to-chlorinity ratio; Humphreys et al., 2022) used when solving the marine carbonate system. Taking this a step further, experimental data like those of Takahashi et al. (1993), Woosley (2021), and Woosley and Moon (2023) could be incorporated when generating these parameterisations. These data constrain the derivatives of the equilibrium constants with respect to temperature, which the parameterisations are not normally optimised for. Constraining these derivatives in the fit would mean that the parameterisation can more accurately calculate derivative properties like υ. Such an approach is common in other fields, such as the Pitzer model for chemical activities in aqueous solutions (Pitzer, 1991). There, different measurable properties of the solution constrain the derivatives of fitted coefficients with respect to temperature, pressure and composition, so a single set of coefficients can be found that can accurately calculate all of these properties (e.g. Archer, 1992).

At present, including derivative properties in parameterisations would be limited by the scarcity of measurements like those of Takahashi et al. (1993). One possible alternative would be to use data like those compiled by Wanninkhof et al. (2022a) and datasets in which the marine carbonate system has been overdetermined (i.e. more than two parameters measured) in the field to better parameterise the equilibrium constants (e.g. Sulpis et al., 2020). A practical challenge with this would be that fCO₂ measurements are archived in different compilations (e.g. SOCAT; Bakker et al., 2016) from other carbonate system variables (e.g. GLODAP; Lauvset et al., 2024), and it is not straightforward to match them back together. However, even if this were overcome, this approach would be less useful for parameterising equilibrium constants than using dedicated experiments in which well-characterised samples are manipulated under a range of different conditions. In the ocean, temperature is correlated with most other properties, so it would be difficult to isolate which of several correlated drivers was responsible for observed variability. Using these field data to generate parameterisations would also mean that they could no longer be used for independent validation. A better approach would be to reserve field data compilations for testing parameterisations that are based on more carefully controlled laboratory experiments. Combining new experimental measurements with field data in this way and interpreting the results through the lens of the new theoretical basis for υ developed in this study will lead to more accurate temperature adjustments of seawater fCO₂ measurements and contribute towards a more internally consistent understanding of the marine carbonate system.

3.3.3 Propagation through to temperature adjustments

Expressions for υ, such as υ_l, υ_q and υ_h, are of interest because they can be used to adjust fCO₂ data to different temperatures without knowing a second carbonate system parameter. The variability in υ discussed so far in Sect. 3.3 is therefore most important in terms of how it affects exp (Υ) and consequent uncertainty and biases in the adjusted fCO₂ values if ignored.

After a 1 °C warming adjustment from mean surface ocean conditions in OceanSODA-ETZH, both the υ_l and υ_p forms have a global mean bias (relative to the change calculated from A_T and T_C) of less than 0.1 µatm, which is negligible. That the global mean bias is close to zero for υ_l is a coincidence, due to Takahashi et al. (1993) using seawater with roughly average υ_Lu00 (Fig. 6a). However, there are stronger regional biases in υ_l, with 95 % of the residuals falling in the range from −1.3 to +1.1 µatm, leading to a global RMSD of 0.8 µatm, or equivalently 0.24 % (Fig. 9a). The van 't Hoff form with parameterised b_h (υ_p) has much lower regional biases, with a global RMSD of 0.16 µatm (0.04 %) and with 95 % of the residuals falling from −0.26 to +0.24 µatm (Fig. 9b). The Takahashi et al. (1993) quadratic fit (υ_q) and the van 't Hoff form with constant b_h fitted to the Takahashi et al. (1993) dataset (${\mathit{\upsilon }}_h^{\text{Ta93}}$) (both not shown) had negligible global mean biases and intermediate regional biases, with RMSD values of around 0.4 µatm.

Figure 9 Bias in adjusted fCO₂ for a 1 °C adjustment ($\mathrm{\Delta }t=+\mathrm{1}$ °C) relative to the change in fCO₂ calculated from A_T and T_C with PyCO2SYS and the Lueker et al. (2000) parameterisation of the carbonic acid dissociation constants when the adjustment is applied using (a) the Takahashi et al. (1993) constant υ_l or (b) the van 't Hoff form with parameterised b_h (Sect. 2.4). The colour scale is the same for panels (a) and (b).

The experimental 1σ uncertainty in exp (Υ_l) and $\exp \left({\mathrm{{\rm Y}}}_h^{\text{Ta93}}\right)$ fitted to the Takahashi et al. (1993) dataset was around 0.04 % for a 1 °C warming adjustment in both cases (Sect. 3.2). For υ_l, the additional uncertainty from spatiotemporal variability that was unaccounted for, as well as the fact that the t–ln (fCO₂) relationship is not linear, is therefore 4 times greater than the experimental uncertainty. For υ_h, the experimental and spatiotemporal uncertainties are about equal in size. Combining these components gives a 1σ uncertainty in adjusted fCO₂ of around 0.06 % for υ_p, or 0.24 % for υ_l. A 1 °C adjustment with either approach could therefore still meet the GOA-ON climate quality goal of 0.5 % (Newton et al., 2015), if other uncertainty components (e.g. in the original fCO₂ measurement at t₀) were small enough. However, the spatiotemporal component of the uncertainty in υ_l results from a known bias with strong regional patterns. Such biases should be eliminated by using a more realistic form of υ (e.g. υ_h rather than υ_l or υ_q, and υ_p or similar to account for variability in υ_h), rather than incorporating them into an overall uncertainty value. However, if applying υ_l globally, 0.24 % is a reasonable value for the component of uncertainty due to a +1 °C adjustment. It is important to remember that this component is a systematic bias with strong regional patterns when combining it with other components that may be independent for each measurement, because these two types of uncertainty have very different statistical consequences when propagated through to other calculations.

The global biases and the amplitude of the spatiotemporal variability scale approximately, but not exactly, with Δt (Fig. 10). The υ_p approach has the lowest RMSD throughout the Δt range, whereas υ_l has the greatest. υ_q and υ_h have intermediate RMSD values, with that for υ_h being significantly lower than that for υ_q for negative Δt. The RMSD values are generally greater for positive Δt than for negative Δt, which emphasises the need to better understand the temperature sensitivity of fCO₂ and the wider marine carbonate system, particularly at higher temperatures (e.g. Sect. 3.1.1; Woosley, 2021). For a Δt of 10 °C, as might be used in studies unravelling the drivers of seasonal fCO₂ changes, the adjusted fCO₂ would have an uncertainty due to spatiotemporal variability of around 0.8 % for υ_p, which is outside the GOA-ON climate quality goal (0.5 %) but well within the weather quality goal of 2.5 % (Newton et al., 2015). However, with υ_l, this component of the uncertainty in adjusted fCO₂ would be over 3 %, also exceeding the weather quality goal.

Figure 10 Effect of Δt on the global RMSD of the residuals in the OceanSODA-ETZH data product between fCO₂ adjusted by Δt using the value calculated from A_T and T_C (υ_Lu00) and the respective υ_l, υ_q, υ_h (all fitted to Ta93) and υ_p.

4 Conclusions

Seawater fCO₂ data can be adjusted to different temperatures using either Eq. (2) or Eq. (4) (which are mathematically identical) with an appropriate expression for Υ, which is υ integrated over the temperature range of interest. The new approach proposed here uses Υ_h as defined in Eq. (21), which has one unknown coefficient (b_h). The value of b_h can be found by fitting to experimental data where available (e.g. Table 1), but spatiotemporal variability in b_h should be accounted for – for example, with the parameterisation in Eq. (35). This approach and the earlier linear and quadratic approaches of Takahashi et al. (1993) have been built into the PyCO2SYS software as of v1.8.3.

There is no theoretical reason to expect a linear relationship between temperature and ln (fCO₂) under isochemical conditions (constant A_T and T_C), as was proposed empirically by Takahashi et al. (1993). Rather, calculations of υ with PyCO2SYS and a theoretical analysis based on the equations governing the marine carbonate system and the van 't Hoff equation each independently suggest that ln (fCO₂) should be proportional to $\mathrm{1}/t_{\mathrm{K}}$ to first order (Sect. 3.1.1). Although this means that the t–ln (fCO₂) relationship does seem to be approximately linear within the temperature range of 2.1 to 24.5 °C measured by Takahashi et al. (1993), it should diverge rapidly from linearity, especially at temperatures above 25 °C, although there are some discrepancies with the very few existing laboratory data which require more measurements to resolve (Sect. 3.1.2). Nevertheless, the new form of the relationship between temperature and fCO₂ proposed here (υ_h with a parameterised b_h coefficient, i.e. υ_p; Eq. 19 and Sect. 2.4) agrees better with field-based datasets than any other approach previously used to adjust seawater fCO₂ to different temperatures (Sect. 3.1.3).

The parameterisation for b_h (and thus υ_p) is based on a specific set of carbonate system parameterisations. We know that these are not perfect (e.g. Woosley and Moon, 2023; Carter et al., 2024), so the b_h parameterisation is an intermediate solution. It does agree very well with the available field data (Wanninkhof et al., 2022a) with which it can be tested (and from which it is independent), but it may need to be updated when new data and/or more accurate parameterisations of the marine carbonate system become available. For some applications, accuracy may be improved by developing parameterisations of b_h that are optimised for any particular region(s) of interest.

The uncertainty in υ depends on the method being used to compute it. The Takahashi et al. (1993) experiment gave a 1σ uncertainty of about 0.04 % °C⁻¹ for the linear form (Sect. 3.2), but this does not account for variability in υ with temperature and seawater composition, nor for the fact that the relationship between temperature and ln (fCO₂) is not linear. A 1σ in υ of 0.23 % °C⁻¹ is therefore more appropriate if using this linear approach (Sect. 3.3.1), which corresponds to 0.24 % for fCO₂ adjusted by +1 °C (Sect. 3.3.3). This uncertainty can be reduced by using the new form of the t–ln (fCO₂) relationship developed here and parameterising its fitted constant (b_h) as a function of temperature, salinity and fCO₂ (Sect. 2.4), which gives an uncertainty of 0.06 % °C⁻¹ for fCO₂ after the same +1 °C adjustment (Sect. 3.3.3). While these uncertainties are probably not large enough to have major consequences for global air–sea CO₂ flux budgets, they are not negligible and may be important in regional studies. They should be incorporated into uncertainty budgets and their systematic components should be eliminated by using a more meaningful and accurate form to model the t–ln (fCO₂) relationship, like that proposed here. However, the uncertainty budget is still incomplete: we cannot propagate uncertainties from the equilibrium constants of the marine carbonate system through to υ, as those uncertainties are not yet adequately defined (Orr et al., 2018). However, if υ can be fitted directly to measurements like those of Takahashi et al. (1993), the uncertainties in the equilibrium constants can be ignored when adjusting fCO₂ values to different temperatures.

More measurements like those conducted by Takahashi et al. (1993) are still needed to empirically confirm whether the υ_h form, and indeed calculations from A_T and T_C, are sufficiently accurate, especially above 25 °C where the deviation from υ_l is expected to be greatest. If further measurements appear to confirm the high-temperature linearity observed by Lee and Millero (1995), there must be either (1) a major error in a parameterisation or a missing component of the equilibrium model for seawater chemistry implemented in software tools like PyCO2SYS at higher temperatures or (2) some technical problem with the experimental approach that means it does not measure quite what we think. More measurements could show whether discrepancies between υ_p and υ_Lu00 are due to needing a more sophisticated equation for υ_h with extra fitted parameters or whether they represent inaccuracies in the Lueker et al. (2000) and other carbonate system parameterisations. In principle, we can measure the effect of temperature on fCO₂ accurately enough to distinguish between these different parameterisations. Coupling such new experiments with other measurements such as pH on the same samples may also help to understand and rectify wider discrepancies in equilibrium models of seawater chemistry. Together with field measurements for independent validation, and using a theoretical basis to interpret the results rather than working only with empirical fits, such a holistic approach should help to better understand υ and its application to adjusting fCO₂ data to different temperatures and should also lead to a more accurate and internally consistent model of the marine carbonate system.

Appendix A: Standard enthalpies of reaction

The reactions corresponding to the equilibrium constants of interest are as follows (Humphreys et al., 2022):

$$\begin{array}{}\text{(AR1)}& {\displaystyle }{\mathrm{CO}}_{\mathrm{2}}\left(\mathrm{g}\right)\stackrel{K_{\mathrm{0}}^{\prime }}{\rightleftharpoons }{\mathrm{CO}}_{\mathrm{2}}\left(\mathrm{aq}\right),\text{(AR2)}& {\displaystyle }{\mathrm{CO}}_{\mathrm{2}}\left(\mathrm{aq}\right)+{\mathrm{H}}_{\mathrm{2}}\mathrm{O}\stackrel{K_{\mathrm{1}}^{\ast }}{\rightleftharpoons }{{\mathrm{HCO}}_{\mathrm{3}}}^{-}+{\mathrm{H}}^{+},\text{(AR3)}& {\displaystyle }{{\mathrm{HCO}}_{\mathrm{3}}}^{-}\stackrel{K_{\mathrm{2}}^{\ast }}{\rightleftharpoons }{{\mathrm{CO}}_{\mathrm{3}}}^{\mathrm{2}-}+{\mathrm{H}}^{+}.\end{array}$$

The standard enthalpies of formation (Δ_fH^⊖) for the relevant species are given in Table A1.

Table A1 Standard enthalpies of formation at 298.15 K following Ruscic and Bross (2023).

Using Hess's law, the standard enthalpies of reaction for Reactions (AR1)–(AR3), as used in Eq. (18), are as follows:

$$\begin{array}{}\text{(A1)}& {\displaystyle }\begin{array}{rl}{\displaystyle }{\mathrm{\Delta }}_{\mathrm{r}}{H}_{\mathrm{0}}^{\ominus }& {\displaystyle }={\mathrm{\Delta }}_{\mathrm{f}}{H}^{\ominus }\left[{\mathrm{CO}}_{\mathrm{2}}\left(\mathrm{aq}\right)\right]-{\mathrm{\Delta }}_{\mathrm{f}}{H}^{\ominus }\left[{\mathrm{CO}}_{\mathrm{2}}\left(\mathrm{g}\right)\right]\\ {\displaystyle }& {\displaystyle }=-\mathrm{19}\mathrm{720}\mathrm{J}{\mathrm{mol}}^{-\mathrm{1}},\end{array}\text{(A2)}& {\displaystyle }\begin{array}{rl}{\displaystyle }{\mathrm{\Delta }}_{\mathrm{r}}{H}_{\mathrm{1}}^{\ominus }=& {\displaystyle }{\mathrm{\Delta }}_{\mathrm{f}}{H}^{\ominus }\left({{\mathrm{HCO}}_{\mathrm{3}}}^{-}\right)+{\mathrm{\Delta }}_{\mathrm{f}}{H}^{\ominus }\left({\mathrm{H}}^{+}\right)\\ {\displaystyle }& {\displaystyle }-{\mathrm{\Delta }}_{\mathrm{f}}{H}^{\ominus }\left[{\mathrm{CO}}_{\mathrm{2}}\left(\mathrm{aq}\right)\right]-{\mathrm{\Delta }}_{\mathrm{f}}{H}^{\ominus }\left({\mathrm{H}}_{\mathrm{2}}\mathrm{O}\right)\\ {\displaystyle }=& {\displaystyle }\mathrm{9134}\mathrm{J}{\mathrm{mol}}^{-\mathrm{1}},\end{array}\text{(A3)}& {\displaystyle }\begin{array}{rl}{\displaystyle }{\mathrm{\Delta }}_{\mathrm{r}}{H}_{\mathrm{2}}^{\ominus }=& {\displaystyle }{\mathrm{\Delta }}_{\mathrm{f}}{H}^{\ominus }\left({{\mathrm{CO}}_{\mathrm{3}}}^{\mathrm{2}-}\right)+{\mathrm{\Delta }}_{\mathrm{f}}{H}^{\ominus }\left({\mathrm{H}}^{+}\right)\\ {\displaystyle }& {\displaystyle }-{\mathrm{\Delta }}_{\mathrm{f}}{H}^{\ominus }\left({{\mathrm{HCO}}_{\mathrm{3}}}^{-}\right)\\ {\displaystyle }=& {\displaystyle }\mathrm{14}\mathrm{702}\mathrm{J}{\mathrm{mol}}^{-\mathrm{1}}.\end{array}\end{array}$$

Summing Eqs. (A1)–(A3) as required by Eq. (18) and propagating the uncertainties listed in Table A1 gives a value of ${\mathrm{\Delta }}_{\mathrm{r}}{H}_{\mathrm{2}}^{\ominus }-{\mathrm{\Delta }}_{\mathrm{r}}{H}_{\mathrm{0}}^{\ominus }-{\mathrm{\Delta }}_{\mathrm{r}}{H}_{\mathrm{1}}^{\ominus }=\mathrm{25}\mathrm{288}$ ± 98 kJ mol⁻¹.