Biogeosciences Technical note : Consistent calculation of aquatic gross production from oxygen triple isotope measurements

Oxygen triple isotope measurements can be used to calculate aquatic gross oxygen production rates. Past studies have emphasised the appropriate definition of the 17O excess and often used an approximation to derive production rates from the17O excess. Here, I show that the calculation can be phrased more consistently and without any approximations using the relative 17O/16O and18O/16O isotope ratio differences (delta values) directly. I call this the “dual delta method”. The17O excess is merely a mathematical construct and the derived production rate is independent of its definition, provided all calculations are performed with a consistent definition. I focus on the mixed layer, but also show how time series of triple isotope measurements below the mixed layer can be used to derive gross production. In the calculation of mixed layer productivity, I explicitly include isotopic fractionation during gas invasion and evasion, which requires the oxygen supersaturation s to be measured as well. I also suggest how bubble injection could be considered in the same mathematical framework. I distinguish between concentration steady state and isotopic steady state and show that only the latter needs to be assumed in the calculation. It is even possible to derive an estimate of the net production rate in the mixed layer that is independent of the assumption of concentration steady state. I review measurements of the parameters required for the calculation of gross production rates and show how their systematic uncertainties as well as the use of different published calculation methods can cause large variations in the production rates for the same underlying isotope ratios. In particular, the17O excess of dissolved O 2 in equilibrium with atmospheric O2 and the17O excess of photosynthetic O 2 need to be re-measured. Because of these uncertainties, all calculaCorrespondence to: J. Kaiser (j.kaiser@uea.ac.uk) tion parameters should always be fully documented and the measured relative isotope ratio differences as well as the oxygen supersaturation should be permanently archived, so that improved measurements of the calculation parameters can be used to retrospectively improve production rates.

1 Introduction Luz et al. (1999) first suggested that the triple-isotope composition of atmospheric oxygen (O 2 ) could be used as a tracer of biological productivity.They showed that photosynthetic O 2 has a small, but measurable excess of the oxygen isotope 17 O with respect to atmospheric O 2 , after normalisation for 18 O/ 16 O isotope ratio differences.The magnitude of the 17 O excess ( 17 ) depends on the chosen normalisation, which is meant to account for so-called mass-dependent isotope fractionation (see Sect. 2.2).
It is clear that stratospheric isotope exchange reactions between ozone (O 3 ) and carbon dioxide (CO 2 ) are responsible for an enhanced 17 O isotope transfer over and above a massdependent relationship with 18 O (Yung et al., 1991(Yung et al., , 1997)).Specifically, the relative 17 O/ 16 O isotope ratio difference of stratospheric CO 2 to tropospheric CO 2 is about 1.7 times that of the 18 O/ 16 O isotope ratio difference (Lämmerzahl et al., 2002), significantly higher than the factor of 0.516 ± 0.015 expected for many mass-dependent fractionation processes (Kaiser, 2008).Since O 2 is the source of the oxygen atoms in the short-lived O 3 molecule, this leads to corresponding 17 O depletion in atmospheric O 2 .
In the aquatic realm, the relative isotope ratio difference between atmospheric and photosynthetic oxygen can be used to calculate gross oxygen production in the mixed layer, using the gas exchange rate as a "timekeeper".Based on the simplified budget by Luz and Barkan (2000), the following equation should hold for the ratio (g) between gross O 2 production (P ) and gross oxygen influx from the atmosphere (kc sat ), where k is the gas exchange coefficient and c sat the saturation concentration of O 2 : g ≡ P kc sat = 17 − 17 sat 17 P − 17  (1) 17 , 17 sat and 17 P refer to the 17 O excess of dissolved O 2 , dissolved O 2 in equilibrium with the atmosphere and photosynthetic oxygen, respectively.The prime symbol distinguishes area-based rates (P ) from volume-based rates (P ).The mathematical treatment in the following is framed in terms of volume-based rates.They can be converted to areabased rates by multiplication by the depth interval for which they are calculated, e.g. the mixed-layer depth for mixedlayer production rates.Equation (1) has been used in numerous studies to calculate gross O 2 production rates (e.g.Sarma et al., 2005Sarma et al., , 2006aSarma et al., , 2008;;Juranek and Quay, 2005;Stanley et al., 2010;Luz andBarkan, 2005, 2009).Usually, statistical uncertainties in the measurement of 17 and in the calculation of k from wind speeds are cited as main contributors to the overall uncertainty in P , which may be between 15 % (Stanley et al., 2010) and 40 % (Quay et al., 2010;Reuer et al., 2007).To eliminate the uncertainty introduced by k, I focus here on the calculation of the dimensionless gross production variable g, which is independent of k.
The advantage of the oxygen triple isotope technique over 18 O/ 16 O isotope ratio measurements in determining production is that the calculated rates are independent of the respiratory isotope effect that is not well known and that would otherwise lead to significant uncertainties (Hendricks et al., 2004;Quay et al., 1993;Venkiteswaran et al., 2008).If implemented correctly, the additional information from the 17 O/ 16 O isotope ratio measurement allows elimination of the respiratory isotope effect from the production calculation.
However, different studies have used different definitions of 17 , without appropriately adjusting 17 sat and 17 P .The use of different definitions of 17 means that the same measurements will give different results for g, causing systematic uncertainty.
Moreover, in an effort for a more rigorous derivation of Eq. ( 1), Hendricks et al. (2004) demonstrated that calculations of g based on this equation were in error.The authors solved the exact equations iteratively, making assumptions for certain parameters and using the biological oxygen supersaturation (O 2 /Ar) (Kaiser et al., 2005) as additional constraint.Specifically, the ratio of (O 2 /Ar) and g was assumed to be equivalent to the ratio of net (N) to gross oxygen production, i.e. f = N/P = (P -R)/P = (O 2 /Ar)/g, where R stands for respiration.The same iterative approach was adopted in subsequent studies (Hendricks et al., 2005;Reuer et al., 2007;Juranek and Quay, 2010;Quay et al., 2010).Quay et al. (2010) and Juranek and Quay (2010) stated that the iterative approach gave on average 10 % higher values for g than Eq. ( 1), without exploring the underlying reasons.None of these iterative calculations considered the effect of inconsistent 17 definitions and the uncertainty in the input parameters used in the calculation of g.
The main goal of the present study is to explore the systematic uncertainty in the calculation of g from triple isotope measurements in dissolved O 2 .This will take into account methodological differences between past studies and the uncertainty in parameters required for the calculation.The scope of the study is limited to these aspects and neither extends to the mass-spectrometric measurement uncertainty in 17 , nor to the uncertainty in the gas exchange coefficient k, nor other systematic errors of the triple oxygen method such as the neglect of horizontal and vertical transport.
I first compare different definitions of the oxygen isotope excess 17 (Sect.2), followed by a derivation of solutions to mass balance equations for dissolved O 2 and its isotopologues in systems without (Sects.3.1 to 3.4) and with (Sect.3.5) gas exchange.These equations are derived without approximations to contrast them with previously published versions.I also show that g can be derived from triple isotope measurements without recourse to iterative solutions or assumptions with respect to f (Sect.4).Then I assess the systematic uncertainty due to the input parameters (Sect.6.1).This is followed by a comparison of g values calculated by different published methods from a range of synthetic data (Sect.6.2).I will make suggestions with respect to which input parameters need to be constrained better, to reduce systematic uncertainties in the calculation of g.I also show how environmental data can help constrain these parameters (Sect.6.3).All uncertainties stated here represent one standard deviation of the mean.Typesetting conventions would require all physical quantities to be in italics.However, for technical reasons this was not possible for the capital Greek delta designating the triple isotope excess 17 .

Notation
Isotope ratio differences (δ values) of a sample relative to a reference are defined as follows with the solidus ( / ) separating the species of interest in numerator and denominator on the right hand side of the equation.
Using the oxygen isotopes 17 O and 16 O as an example, the isotope-amount ratio r (or, shorter, the isotope ratio) is defined as where n stands for the amount of substance.For clarity and simplicity, I therefore use a notation where only the minor isotope is listed as a left superscript index (just like in nuclide notation), with the species of interest given as a right subscript index, e.g. 17 δ sample/reference .
It is common practice to use atmospheric oxygen ("Air-O 2 ") as reference material for dissolved O 2 in aquatic systems.The use of atmospheric oxygen as international measurement standard has been endorsed by the Commission on Isotopic Abundances and Atomic Weights (Wieser and Berglund, 2009), a commission under the Inorganic Division of the International Union of Pure and Applied Chemistry (IUPAC).When Air-O 2 is the reference in the following, I will also omit it from the quantity symbol, e.g. 17 δ sample .
Finally, when the sample of concern is dissolved O 2 in water, the corresponding index is also omitted from the quantity symbol, e.g. 17 δ.

Quantification of deviations from mass-dependent isotope ratio relationships
As discussed by Kaiser et al. (2004), deviations from massdependent isotope ratio relationships have been defined using four functional relationships between 17 δ and 18 δ.≡ ln(1+ 17 δ)−λ ln(1+ 18 δ) (Angert et al., 2003) Superscript indices such as " †" have been added to distinguish between different 17 definitions.I do not make a distinction between the symbols ( 17 O), 17 O and 17 by way of definition.
The coefficients κ and λ are meant to reflect the "expected" mass-dependent isotope fractionation, but strictly speaking their choice is entirely arbitrary, as these are merely definitions.Generally, κ and λ may be derived from empirical relationships, e.g.κ = 0.515 in the case of a study on N 2 O (Cliff and Thiemens, 1997) or λ = 0.5279 for meteoric waters (Barkan and Luz, 2007), or may be based on theoretical predictions (Young et al., 2002;Kaiser, 2008).
The designation of 17 values as "isotope anomalies" may be misleading, especially when the 17 values are small, because non-zero 17 values might just be due to the way they were defined.More neutral terms such as " 17 O excess" (Kaiser et al. 2003;Angert et al., 2004;Barkan and Luz, 2007) and " 17 O balance" (Kaiser, 2008) have been suggested." 17 O excess" has been adopted most widely and I will use this term here for 17 values.
From a theoretical point of view, Eq. ( 6) is the most satisfactory because it obeys the basic isotope delta "addition theorems", e.g.δ A/C = δ A/B + δ B/C + δ A/C δ B/C .However, Eq. ( 4) also has merits because of its mathematical simplicity and ease of use with mass-balance and mixing calculations.
Per se, none of the definitions or coefficients is better or worse than others -all of them are merely mathematical constructs.However, different definitions give different 17 values for the same underlying 17 δ-18 δ pairs.Any subsequent calculations or manipulations have to bear this in mind and follow a consistent mathematical treatment.Moreover, any 17 value should not be cited in isolation, i.e. not without its definition and, crucially, not without the corresponding 18 δ and/or 17 δ values.These caveats have not always been followed in the past, which, as we will see below, is partly responsible for different g values obtained for the same 17 values, depending on the calculation method.
For the mathematical treatment of the budget equations in the present paper I choose Eq. ( 4) because it simplifies the discussion.Unless required for clarity, I drop the index † from 17 † .I adopt κ = 0.5179, based on the weighted average ratio between the 17 O/ 16 O and 18 O/ 16 O isotope fractionations during respiration (Sect.5.1) so that 17 If a coefficient other than 0.5179 is used, I will indicate this explicitly, e.g. 17 (κ = 0.521), corresponds to Eq. ( 8), but with the coefficient 0.521.

Units
The 17 O isotope excess ( 17 ) of O 2 in the atmosphere and aquatic environment is always less than 10 −3 .The measurement precision is usually between 1 × 10 −6 and 8 × 10 −6 .17 values of O 2 are therefore conveniently expressed in multiples of 10 −6 , for which Luz et al. (1999) used the symbol "per meg", following a similar practice adopted for gas delta values related to O 2 /N 2 ratios (Keeling et al., 1998).The same symbol was also chosen by subsequent studies on oxygen triple isotopes.
However, "per meg" appears to be an awkward replacement for the symbol "ppm" (short for "parts per million").The symbol "ppm" has traditionally been used to represent the value 10 The symbol "per meg" has no such recognition.Therefore, I will use the symbol ppm with 17 values.For other δ values and isotope fractionations (designated with the symbol ε), I will use multiples of 10 −3 , abbreviated ‰ (short for "per mill" or "parts per thousand").

Budget calculations
In the following, I give a consistent mathematical treatment of isotope budgets in increasing order of complexity, with a view to derive gross oxygen production (P ) and to gain an understanding of how uncertainties in the calculations parameters can affect P .Firstly, I discuss respiration and production on their own, then the combination of both processes, in general and under isotopic steady-state conditions, and finally explain how both processes can be combined with diffusive and bubble-mediated gas exchange.
Extensive properties (such as rates and concentrations) without an index refer to the isotope 16 O only.They can be related to the total value (sum over all isotopes) using the reference isotope ratios (r r ) and δ values.For example, the 16 O production ( 16 P = P ) is related to the total production (P total ) via the following relationship: Since 17 r r ≈ 0.000387, 18 r r = 0.002053 (Kaiser, 2008), 17 δ P 1 and 18 δ P 1, P ≈ P total .The correction from P to P total is less than 0.25 % and therefore negligibly small compared with other uncertainties that enter into the calculation of P from triple isotope measurements.

Respiration only
One of the simplest budgets comprises respiration only.It has the mass balance equation The concentration of the major isotope 16 O is represented by the symbol c.The 16 O respiration (R) is assumed to be of zeroth order, i.e. independent of the oxygen concentration (Bender, 1990).The corresponding equation for 17 O is (11) The ratio between 17 R and R = 16 R is assumed to follow the isotope distribution in dissolved O 2 , but modified by a respiratory isotope effect 17 ε R (Bender, 1990): Substituting this into Eq.( 11) and using δ notation (omitting the index "17" from δ) gives The derivative is expanded and the constant r r cancelled on both sides: Equation ( 10) is substituted into the preceding equation to give and which can be integrated to the well-known Rayleigh fractionation equation: As pointed out by Angert et al. (2003), the resulting 17 # value is With λ set equal to the ratio of the 17 O/ 16 O and 18 O/ 16 O isotope fractionations (γ R ), i.e.
This led Angert et al. (2003) and Luz and Barkan (2005) to suggest that a 17 definition following Eq.( 7) with λ = γ R would be more appropriate than others because it removes the influence of the respiratory isotope effect on the measured 17 O excess.However, this assertion fails when production is included in the oxygen budget, as shown in Sects.3.2 and 3.3.

Production only
An oxygen budget that includes production, but not respiration, is given by with the corresponding isotopic relationship Combining the two equations gives which can be integrated to In this case, there is no simple relationship between ln (1 +17 δ) and ln (1 + 18 δ) that would single out Eq. ( 7) over other possible definitions of the 17 O excess, no matter what λ value is chosen.However, Eq. ( 4) allows writing down a straightforward relationship between the corresponding 17 O excesses: This shows that as the concentration c increases due to production, the influence of the initial composition ( 17 0 ) decreases and 17 approaches 17 P asymptotically.

Production and respiration
The following budget equation combines production and respiration: This equation can be readily integrated to give the concentration c: with the ratio of net to gross production f = N /P = (P -R)/P .The corresponding equation for the minor isotope is and combination of the equations for major and minor isotopes gives This can be integrated to for P = R (and therefore f = 0) and to for P = R (and therefore f = 0).There is no simple relationship between ln (1 + 17 δ) and ln (1 + 18 δ) or 17 δ and 18 δ that would single out a certain definition of 17 over other possible definitions.

Isotopic steady state between production and respiration
The case dδ/dt = 0 corresponds to "isotopic steady state".Then, for any combination of P and R, it follows from Eq. ( 28) Isotopic steady state can be attained even when the concentrations vary (i.e. if P = R and therefore f = 0).In particular, for t → ∞, the δ value attained according to Eqs. ( 29) and ( 30) is equal to δ S .However, the concentration may still increase or decrease according to Eq. ( 26).
The steady-state 17 O isotope excess could be defined with Thus, with a choice of λ appropriate for a certain value of f , it may be possible to argue that the corresponding 17 # definition is preferable because the respiration isotope effect does not appear in 17 # S .Obviously, this is a misleading conclusion, since the value of f is usually not known.Note that the value of 17 # P would also depend on the value adopted for f if this "tuned" value of λ was adopted in the definition of the 17 O excess.
Luz and Barkan (2005) argued that the global biosphere could be considered in steady state (P = R, and therefore f = 0) and that a 17 # definition with λ = ln (1 + 17 ε R )/ln (1 + 18 ε R ) would be the most suitable in this case because 17 # S would always be equal to 17 # P , independent of the respiration fractionation.This is in contrast to their suggestion that λ = 17 ε R / 18 ε R , should be chosen in other cases (see Sect. 3.1) and confirms the notion that there is no definition of the 17 O excess that is inherently "better" than others and that it is essentially possible to adopt any definition.
The fact that the 17 O excess in isotopic steady state is dependent on the f ratio has important consequences for the triple isotope technique because it means that the 17 O excess is not only influenced by production and gas exchange, but also by respiration.This was recognised by Hendricks et al. (2004) who adopted an iterative approach to derive g.A simpler approach to derive g is shown in Sect.4.2.1.
To illustrate the effect of f on the steadystate 17 O excess, I use a numerical example with  7), i.e. 17 # , and Eq. ( 4) i.e. 17 , with triple oxygen isotopes cannot rely on the 17 O excess alone, but has to consider the underlying 17 δ and 18 δ values (see Sect. 5.2).

Production, respiration and gas exchange
I now consider production and respiration together with diffusive gas exchange; the model generally adopted for gross oxygen production calculations from triple oxygen isotopes.
A similar mass balance equation was used by Hendricks et al. (2004), but the latter study did not consider the influence of isotopic fractionation during oxygen invasion from and evasion to the atmosphere.Following Luz et al. (2002), I include these fractionations here explicitly.I also show how the model could be extended to include bubble injection, which was mentioned but apparently disregarded in the gross production calculations of Stanley et al. (2010).
The mass balance equation for the major isotope (without bubble injection) is where the gas exchange frequency of the mixed layer is the ratio of gas exchange coefficient and mixed layer depth, i.e.
The corresponding equation for the minor isotope is where ε E and ε I are the isotopic fractionations during evasion and invasion, respectively.The combination of both equations gives Substituting f = (P -R)/P , g = P /(ν mix c sat ) and the super- This can be re-arranged to isolate the δ value: The isotopic fractionations during evasion and invasion are related via the δ value of dissolved oxygen at saturation, δ sat , such that 1 + With ε E = 0 we recover from Eq. ( 36) the corresponding equation of Hendricks et al. (2004): In Sect.4.2, I discuss how Eq. ( 36) can be used to derive g.Bubble-mediated transfer can contribute to air-sea exchange.Usually, two bubble transfer mechanisms are distinguished: bubble injection due to complete dissolution of small bubbles and bubble exchange due to partial dissolution of larger bubbles, with bubble exchange contributing only 0 to 10 % to the total bubble flux (Stanley et al., 2009).It is relatively straightforward to include bubble injection in the O 2 mass balance.The smaller contribution from bubble exchange is neglected here, but could be treated in a mathematically similar way to diffusive gas exchange.The amended mass balance equation is: where F inj is the air injection flux and χ is the mixing ratio of atmospheric O 2 .Since the δ values are expressed relative to atmospheric O 2 , the corresponding mass balance equation for the minor isotope (Eq.36) has the same additional term F inj χ.This gives

Calculation of gross production rates
In this section, I explain how the mass balance equations derived in Sects.3.3 and 3.5 are used to compute gross oxygen production below and within the mixed layer.

Production below the mixed layer
Disregarding vertical and horizontal transport, oxygen below the mixed layer is only influenced by production and respiration (Eq.28).By measuring temporal changes in the isotope composition of dissolved O 2 , it is possible to determine gross production (Luz and Barkan, 2009).The corresponding two budget equations for 17 O/ 16 O and 18 O/ 16 O isotope ratios are combined to eliminate R and to compute P : Using 17 = 17 δκ 18 δ, this equation can also be re-written as Compare this with the approximation given by Luz and Barkan (2009) in their Eq.( 5) (re-written for a single depth and corrected for two errors in the numerator -the index of their second term should be "in" -here replaced by index "0" -and the third term should be subtracted rather than added): which can be written in non-discretised form, In other words, the approximate solution only agrees with the exact solution in Eq. ( 43) if κ = γ R (1 + 17 δ) / (1 + 18 δ).Since 17 δ 1 and 18 δ 1, the approximate solution will often be sufficiently precise if we choose κ = γ R .However, there does not appear to be any advantage in using the approximate solution because 17 δ and 18 δ are available anyway and used to compute 17 .I therefore suggest choosing Eq. ( 42) to compute P based on oxygen triple isotope measurements below the mixed layer.

Production within the mixed layer
In the mixed layer, we have to consider production, respiration and gas exchange in the oxygen budget (again, neglecting horizontal and vertical transport).In principle, 18 δ or 17 δ alone would be sufficient to derive g from Eq. ( 36): This equation requires the temporal trend in the isotopic composition, dδ/dt, to be known.Often, the temporal trend is neglected because it is small.This allows calculating g based on the analysis of the oxygen triple isotope composition of a single sample.Then, we have Note that concentration steady state, i.e. dc /dt = 0, is not required.
Most parameters (ε E , ε I , δ P ) and variables (s, δ) can be measured with sufficient precision to calculate g.However, ε R is often not sufficiently well-constrained to give precise results for g (Quay et al., 1993).Moreover, f has to be estimated from s and g via f = s/g (Hendricks et al., 2004), which can introduce additional uncertainties, because the relationship f = s/g is derived from assuming dc/dt = 0 and therefore P -R = sν mix c sat , cf.Eq. ( 33).

Direct calculation of g from 17 δ, 18 δ and O 2 supersaturation s
Just as for production below the mixed layer, oxygen triple isotopes allow elimination of ε R .Based on the 17 O and 18 O equivalents of Eq. ( 47) we obtain where γ R = 17 ε R / 18 ε R .Eliminating ε R has also removed the (1-f ) term because it always appears as a product with ε R .

Iterative calculation of g from 17 and O 2 supersaturation s
With the additional assumption of concentration steady state, i.e. dc/dt = 0 (and therefore s = gf ), it is possible to determine g based on only two variables, for example, 17 and s.This corresponds to the iterative approach used by Hendricks et al. (2004).In this case, 18 δ is calculated using the following equation derived from Eq. ( 37) and an initial guess for g: Then, 17 δ is derived from 17 and 18 δ and an improved value of g is calculated via Eq.( 48).All steps are repeated until g converges.If the assumption dc/dt = 0 is valid, this approach gives the same result as Eq. ( 48).However, since the use of Eq. ( 48) has fewer caveats and since 17 δ and 18 δ values are available anyway, the use of Eq. ( 48) is preferable.
Often, it is argued that instead of the oxygen supersaturation s, the biological supersaturation ] -1 should be used to calculate g because it corrects s for physical processes such as bubble-mediated gas transfer.While the validity of this argument was rigorously demonstrated for net production (Kaiser et al., 2005), it is not clear that it also applies for gross production.In particular, bubble injection influences the δ values in a different way than diffusive gas exchange (see Sect. 3.5), which makes the use of s bio inappropriate.For the present study, I will disregard the influence of bubble processes, not least because g values calculated according to Eq. ( 48) are not very sensitive to s (for | s | 1, as commonly found in the surface ocean).In contrast, the iterative approach is affected to a larger degree by the choice between s and s bio (see also Sect.6.3).

Non-steady state conditions -calculation of gas exchange coefficients
For completeness, I would like to mention how the "dual delta" calculation method based on 17 δ and 18 δ can be used for non-steady state conditions.This requires the disequilibrium terms, dδ/dt, to be measured.The corresponding equation to calculate g is derived from Eq. ( 46): The dual delta method is also suitable for the calculation of the gas exchange coefficient k based on the diurnal cycle of oxygen triple isotopes (Sarma et al., 2010).This approach assumes that g is zero during the night.Measurements of 17 δ and 18 δ throughout the night can then be used to derive k according to

Input parameters
In this section I review the parameters that are required to calculate g and f according to Eqs. ( 48) and (50).Table 2 gives an overview of the parameters used and their ranges.They will be used to estimate the systematic uncertainty in g.I focus on parameters appropriate for the marine mixed layer because this is where the oxygen triple isotope technique has found the widest application.However, with appropriate adjustments of the parameters, the same technique can also be used for freshwater environments (Luz and Barkan, 2009).

Respiration
The ratio γ R = 17 ε R / 18 ε R has been measured for individual oxygen consumption pathways, single species and community cultures (Luz and Barkan, 2005;Helman et al., 2005;Angert et al., 2003).Luz and Barkan (2005) derived a weighted average value of γ R = 0.5179 ± 0.0006 from these studies and I adopt this value here.For consistency, I also set κ ≡ 0.5179 in the 17 definition (Eq.8).Some oxygen consumption reactions were found to have γ R values outside this range.For example, the Mehler reaction was shown to have γ R = 0.526 ± 0.002 in isolated pea thylakoids and γ R = 0.497 ± 0.004 in a Synechocystis species, indicating different reaction mechanisms in plants and cyanobacteria (Helman et al., 2005).However, in most natural phytoplankton communities the Mehler reaction does not play a role and I do not consider it here.For the respiration fractionation 18 ε R , a large range from −6 ‰ for fish and human respiration to −24 ‰ for the alternative respiratory pathway (AOX) has been reported (Luz and Barkan, 2005;Helman et al., 2005).For the purposes of the present study, I use 18 ε R = (−20 ± 4) ‰, which is more representative of values reported for marine communities: (−22 ± 3) ‰ in the Southern Ocean (Hendricks et al., 2004), (−21 ± 2) ‰ in the Equatorial Pacific (Hendricks et al., 2005), (−22 ± 6) ‰ in the subarctic Pacific (Quay et al., 1993) and (−20 ± 3) ‰ for different marine organisms (Kiddon et al., 1993).
The value of 17 P = (150 ± 13) ppm derived for marine photosynthetic O 2 disagrees with the 17 O excess of (249 ± 15) ppm reported by Luz and Barkan (2000) for photosynthetic O 2 , which has been adopted by all subsequent studies of gross oxygen production using the triple isotope technique.To understand the disagreement, I will now look into how the value of 249 ppm was obtained, bearing in mind that in the year 2000, none of the triple isotope studies referred to in the previous paragraph had been published yet.
Luz and Barkan (2000) measured the 17 value of O 2 produced by cultures of the planktonic alga Nannochloropsis and the coral Acropora, which gave on average 17 P (κ = 0.521) = (249 ± 15) ppm.All subsequent studies using the triple isotope technique have adopted this value, even though they used different definitions of the 17 O excess.This leads to inconsistent results as shown in Sect.6.2.Unfortunately, Luz and Barkan (2000) did not give the 17 δ or 18 δ values corresponding to their reported 17 O excess, so that it is impossible to recalculate the 17 O excess using other definitions, without making some assumptions.
The important aspect to notice here is the large difference between the steady-state 17 S value for f = 0 and 17 P (for f = 1), cf.Fig. 1. 17 S value would only agree with 17 P for an appropriately "tuned" definition of the 17 O excess (see Sect. 3.4).Coincidentally, in the case of Acropora, 17 # S (λ = 0.518) = 236 ppm for f = 0 and the inferred 17 # P (λ = 0.518) value is 246 ppm, i.e. they are quite close to each other.The "tuned" λ value to make them match exactly would be ln (1-0.519× 0.0138)/ln (1-0.0138)= 0.517.In another coincidence, 17 # P (λ = 0.518) = 246 ppm is very close to 17 (κ = 0.521) = 252 ppm reported by Luz and Barkan (2000), which has been the basis of the average value of (249 ± 15) ppm for the photosynthetic 17 O excess used in all subsequent studies of triple isotope-based gross production.
What remains unexplained is the difference between the 17 # P (λ = 0.518) = (221 ± 14) ppm inferred at the beginning of this section and the value of (246 ± 5) ppm inferred based on the Acropora culture (a similar calculation may be possible for Nannochloropsis, but would have similar uncertainties because of the lacking 17 δ and 18 δ values).For the purposes of the present study and for consistency with other studies using the triple isotope technique, I assume The kinetic isotope fractionation during O 2 invasion was measured as 18 ε I = (−2.8± 0.2) ‰ (Knox et al., 1992).I am not aware of any triple oxygen isotope studies of the corresponding 17 O/ 16 O isotope fractionation and therefore adopt 17 ε I = (1 + 18 ε I ) θ -1, with θ = 0.516 ± 0.015, covering the theoretically predicted range for mass-dependent isotope effects (Kaiser, 2008).
The kinetic isotope fractionation during evasion, ε E , is calculated from ε I and the δ value of dissolved O 2 in equilibrium with the atmosphere, δ sat , according to ε E = (ε I δ sat )/(1 + δ sat ). 18δ sat has been reported by Benson et al. (1979) and I use fit (2) in their Table VIII, i.e. 18 δ sat = e −0.00072951+0.42696K/T −1, where T is the thermodynamic temperature.The fit is reported to have an uncertainty of 0.017 ‰ for the temperature range from 0 to 60 • C.
Previous studies have neglected isotopic fractionation during evasion (i.e.ε E = 0) and assumed that ε I = δ sat .I will consider the influence of this assumption on the calculated g value in Sect.6.1. 17δ sat is usually not reported directly, but as 17 O excess with respect to 18 δ sat .There have been a number of measurements of the 17 O excess that do not all agree.The published data are summarised in Table 1.After adjusting the values reported in the literature to a single 17 definition following Eq.( 4) with κ = 0.5179, the 17 sat values at room temperature cluster around two values: 17 to 18 ppm (Luz andBarkan, 2000, 2009;Sarma et al., 2006b;Juranek and Quay, 2005) and 8 to 9 ppm (Stanley et al., 2010;Reuer et al., 2007).The measurement of (13 ± 5) ppm by Sarma et al. (2003) could be reconciled with both values, but has a high error.As already pointed out by Stanley et al. (2010) there is no consistent pattern that could explain this, such as the preparation method (bubbling versus stirring) or the water type used.The 18 δ sat values were not reported in all studies, but mostly agree with those calculated using the parameterisation of Benson et al. (1979).One of the studies in the 8 to 9 ppm-cluster found an 18 δ sat value that is about 0.06 ‰ lower than the parameterisation (Reuer et al., 2007); the other study in the same cluster did not report 18 δ sat (Stanley et al., 2010).Interestingly, the 17 sat value derived from Reuer et al. (2007) for 11 • C (row 4a) is close to that of Luz and Barkan (2009) for 12 • C (row 5b), even though the same two studies disagree near room temperature (rows 4b and 5c).
In theory, bubbling should lead to slightly higher δ values.This is perhaps counter-intuitive because bubble injection adds O 2 with δ = 0, i.e. lower than for dissolved O 2 at saturation, and might therefore be expected to decrease the δ value of the O 2 in the bubbled solution.Bubble exchange should have a lesser effect and is neglected here.
The reason for the enhanced δ value due to bubble injection is that bubbling leads to an enhanced O 2 concentration and thus an additional evasion flux.At equilibrium, the bubble flux F inj has to match the diffusive gas exchange ν mix (c c sat ).Using the terminology from Sect.3.5, we require F inj = ν mix c sat s.The corresponding δ value at equilibrium is given by Whether the corresponding 17 value is higher or lower than 17 sat depends on θ, but for s <2 %, the difference should be less than 1 ppm.If s were greater than 2 %, this would lead to an enhancement of the 18 δ sat value by more than 0.07 ‰, which cannot be reconciled with the published data.Luz and Barkan (2009) have also reported a temperature dependence of the 17 O excess at saturation, i.e. 17 # sat (λ = 0.518)/ppm = 0.6 ϑ/ • C + 1.8.At 24 • C, this gives 17 sat = 16 ppm.Since a value of 16 ppm that has been used for 17 sat in most previous studies, I will also adopt it here for a study of the systematic uncertainty due to the calculation method of g.However, I will also test a scenario with 17 sat = 8 ppm.

Systematic uncertainty of production within the mixed layer
In this section I evaluate the systematic uncertainties in the calculation of mixed-layer g values due to the input parameters (Sect.6.1) and the calculation method (Sect.6.2).The base case calculation and the evaluation of the uncertainty due to the input parameters follow Eq. ( 48).Since the uncertainty in g depends on the values of the parameters as well as their uncertainty, I phrase this section not in terms of formal error propagation, but rather show the relative deviations from the base case, for different input parameters and calculation methods.Finally, I use a published data set of triple oxygen isotope measurements in the Southern Ocean to show how different calculation methods can affect the calculation of g in practice (Sect.6.3).I also demonstrate how concomitant isotope measurements and O 2 supersaturation data can be used to check the consistency of the calculation method and to potentially improve the input parameters.

Input parameters
The values of the input parameters used for the base case are shown in Table 2.The 17 values are shown for reference only, but are not used in the calculations."Synthetic" 17 δ and 18 δ values are computed for a range of g and f values according to Eq. ( 37), assuming dδ/dt = 0 (isotopic steady  Benson et al. (1979).The last column corresponds to the recalculated 17 O excess using the definition given by Eq. ( 8).Uncertainties are ±1 standard error.n = number of measurements.A value in brackets has not been reported and was assumed.A dash (-) means that the corresponding parameter was not reported.Table 2. Input parameters used as base case in the calculation of g (Sect.6.1) and their uncertainties (Sect.5).All δ values are relative to Air-O 2 .The 17 values are defined as 17 = 17 δ -0.5179 18 δ (cf.Eq. 8) and expressed relative to Air-O 2 .However, they are not needed for the calculation according to Eq. ( 48) and are listed for reference only.All values have been adjusted to the same decimal for clarity, irrespective of their actual uncertainty.

Quantity
state).The oxygen supersaturation s is assumed to correspond to concentration steady state and is calculated from Eq. ( 33) with dc/dt = 0, so that s = gf .
These synthetic 17 δ and 18 δ values are then used to derive g according to Eq. ( 48), with one input parameter at a time increased or decreased by the corresponding uncertainty stated in Table 2.Not all parameters listed in Table 2 need to be considered: 18 ε R and 17 ε R have been eliminated from Eq. ( 48) and are therefore irrelevant.The parameters 17 δ P , 17 δ sat , 17 ε I , 18 ε E and 17 ε E are calculated from other values in Table 2 and are therefore also disregarded.This leaves the following seven parameters to test how much their associated uncertainties contribute to systematic errors in g: γ R , 18 δ P , 17 P , 18 ε I , θ,  2) for different parameters in Eq. ( 48).Panel (a) corresponds to g = 0.4 and a range of f from −1.0 to +1.0 (negative values correspond to net heterotrophy, positive value to net autotrophy).Panel (b) corresponds to f = 0.1 and range of g from 0.01 to 10 (logarithmic axis).
To cover a large range of environmental conditions, I consider two cohorts of synthetic 17 δ and 18 δ values corresponding to (a) "g = 0.4, varying f ": a fixed g value of 0.4 with f varying from −1.0 to +1.0 (giving s values from −40 % to +40 % and 17 values from 46 to 71 ppm) and (b) "f = 0.1, varying g": a fixed f value of 0.1 with g varying from 0.01 to 10 (giving s values from 0.1 % to 100 % and 17 values from 18 to 185 ppm).The range used here is larger than typically encountered for oceanic mixed layer conditions, for which f is more likely to be in the range from −0.1 to +0.4 and g in the range from 0.01 to 1.However, under certain conditions, this range may be exceeded and I have therefore chosen to cover a wider range of f and g values.For example, recent work by Prokopenko et al. (2011) on gross oxygen production during a bloom in the coastal Bering Sea has found f values up to 1 and g values greater than 1.I express the g values calculated according to different uncertainty scenarios in terms of their relative deviation from the base case g values.Since the absolute deviations from the base case scale approximately with g, this means that the relative deviations for the "g = 0.4, varying f " are also representative for other values of g.The results are shown in Fig. 2.
The systematic uncertainties due to 18 δ P (corresponding to the 18 δ value of the source water and the photosynthetic isotope effect), 18 ε I and 18 δ sat are negligible: the relative deviations from the base case are always less than 5 % for 18 δ P (for an uncertainty of ±0.5 ‰), less than 0.5 % for 18 ε I (for ±0.2 ‰), and less than 1.5 % (for ±0.017 ‰).Only when isotopic fractionation during gas evasion is completely neglected (i.e.ε I = δ sat , ε E = 0; the default for previous studies) do the g values deviate noticeably from the base case.However, deviations greater than 10 % are only reached for | f | > 0.5 or g >9.
Since 17 P and 17 sat enter directly into the approximated calculation of g according to Eq. (1), it is not surprising that their uncertainties lead to the largest relative deviations from the base case and therefore the largest systematic uncertainty in g.For the "f = 0.1, varying g" cohort (Fig. 2b), the relative deviations from the base case exceed 15 % for g < 0.1 due to the 2 ppm uncertainty in 17 sat and exceed Table 3.Comparison between different calculation methods for g.A dash (-) or values in brackets mean that the corresponding parameters are not used in the calculation.The "used" 17 O excess values are used by the different calculation methods.The "implied" 17 O excess values are calculated using the definitions adopted by the different calculation methods, based on the listed 17 δ P , 18 δ P , 17 δ sat and 18 δ sat values.Where the calculation method does not require these δ values, the values for the "best case" in Table 2  15 % for g > 1 due to the 15 ppm uncertainty in 17 P .For the "g = 0.4, varying f " (Fig. 2a) cohort, the relative deviation in the oceanographically most relevant range of f from −0.1 to +0.4 is less than 15 %.However, for both cohorts, the relative deviations in g reach 20 % or more for the special cases (1) ( 17 P = 150 ppm) and (3) ( 17 sat = 8 ppm).
It is perhaps surprising to see the noticeable effect of the uncertainty in γ R and θ on g, in particular since the value of γ R = 0.5179 is thought to have an uncertainty of only ±0.0006.The relative deviations due to γ R are not symmetric about the x-axis because γ R enters Eq. ( 48) in both numerator and denominator.For | f | >0.4 or g > 1, the relative deviations of g from the base case can exceed 15 %; however, they stay below 15 % for the oceanographically most relevant ranges of f and g.The uncertainty in θ has a particularly noticeable effect on g for the "g = 0.4, varying f " cohort (Fig. 2a).This can be explained by looking at Eq. ( 49), an equivalent formulation of Eq. ( 48).The bracketed term in the numerator of Eq. ( 49) can be approximated by (θ − γ R ) 18 ε I − 17 sat .Therefore, the further away from zero the O 2 supersaturation s is, i.e. the further away from zero f is, and the larger the difference between θ and γ R , the more pronounced the effect on g.
In summary, the evaluation of the uncertainty in g caused by the input parameters shows that, in principle, these parameters can be measured to within a range that does not cause systematic uncertainty >15 % for typical oceanic mixed layer conditions (−0.1< f < +0.4,0.01 < g <1).This means that the systematic uncertainty in g would not contribute more to the overall uncertainty of gross production P than the lower end of uncertainty estimates for the gas exchange coefficient k (Stanley et al., 2010).However, independent measurements of 17 P = 150 ppm and 17 sat = 8 ppm are not covered by the stated uncertainties in the input parameters and indicate a need for these parameters to be re-measured.Also, for more "extreme" values of f and g, the triple oxygen isotope method is significantly impaired by the uncertainty in the input parameters of the calculation.

Calculation method
In this section, I evaluate the impact of different calculation methods on the derived value of g for the same set of synthetic "measurements" of 17 δ, 18 δ and oxygen supersaturation s used in Sect.6.1.None of the previously published works has used the exact input parameters I have adopted  3).Panel (a) corresponds to g = 0.4 and a range of f from −1.0 to +1.0 (negative values correspond to net heterotrophy, positive value to net autotrophy).Panel (b) corresponds to f = 0.1 and range of g from 0.01 to 10 (logarithmic axis).Black curves correspond to calculation methods based on Eq. ( 1).Red curves correspond to iterative methods.
here (Table 2) although two of them are close (Juranek and Quay, 2010;Quay et al., 2010).I have compiled examples of different calculation methods from the literature in Table 3.Even though the same "measurements" were used for all calculation methods, the 17 O excess values differed, depending on the definition of the 17 O excess adopted in the particular study.
Most studies have used Eq. ( 1) with 17 P (or 17 # P ) = 249 ppm and 17 sat (or 17 # sat ) = 16 ppm, but varying coefficients in the definition of 17 .In none of the studies were the input parameters adjusted to the definition adopted for the 17 O excess.The effect of this is shown in the bottom rows of columns 1, 2, 3 and 8 in Table 3: for the same underlying 17 δ P , 18 δ P , 17 δ sat and 18 δ sat values, the resulting "implied" 17 P and 17 sat values differ by up to 116 ppm and 3 ppm, respectively.
A few studies have adopted an iterative approach (Sect.4.2.2) based on the 17 O excess and the biological oxygen supersaturation s bio (columns 4 to 6 in Table 3).Different λ values in the definitions of the 17 O excess and different 18 δ P , 18 ε R and 18 δ sat values were used (Hendricks et al., 2004;Reuer et al., 2007;Juranek and Quay, 2010).Because of this, even studies using essentially the same 17 # P and 17 # sat values give different results for g.The distinction between s and s bio is irrelevant for the present set of synthetic measurements.Iterative method and a calculation based on Eq. ( 48) therefore give the same value for g.However, the distinction between s and s bio may be relevant for environmental samples, in particular if f values are derived from the ratio of g and s in case of the iterative method.
Columns 7 and 8 of Table 3 correspond to the base case calculation according to Eq. ( 48) and an approximated calculation following Eq.( 1), with 17 P and 17 sat and the 17 definition made consistent with the base case.
Figure 3 illustrates the relative difference between g calculated by the different methods and the base case.For a large range of underlying f and g values, the g values calculated according to Eq. ( 1) are >25 % below the base case (Luz and Barkan, 2000;Sarma et al., 2005;Juranek and Quay, 2005).This is not due to an approximation error in the derivation of Eq. (1) as illustrated by the reasonable agreement between base case calculation and approximation if 17 P and 17 sat parameters consistent with the base case calculation are used ("this paper, approx.").In the case of Luz and Barkan (2000) the difference arises due to the choice of κ = 0.521, which does not agree with γ R = 0.5179 ± 0.0006.In the case of Sarma et al. (2005) and Juranek and Quay (2005)  oxygen production rates.Clearly, in addition to allowing a more consistent calculation of g (Sect.4.2.1), the pairs of 17 δ and 18 δ values provide a wealth of additional information that cannot be subsumed by a single 17 O excess value that, after all, is simply a mathematical construct.

Conclusions
I have reformulated the calculation of mixed-layer gross oxygen production in terms of relative isotope ratio differences ( 17 δ and 18 δ).This so-called dual delta method is based on Eq. ( 48).It avoids mathematical approximations and iterative solutions and gives an explicit result in terms of a dimensionless gross oxygen production rate g.In addi-tion to the parameters identified previously (i.e. the ratio of 17 O/ 16 O respiration fractionation to 18 O/ 16 O respiration fractionation, γ R = 17 ε R / 18 ε R ; the 17 O excess for photosynthetic O 2 , 17 P ; and the 17 O excess of O 2 in water saturated with atmospheric O 2 , 17 sat ), g also depends (to a lesser extent) on the 18 O/ 16 O isotope delta of photosynthetic O 2 ( 18 δ P ), the 18 O/ 16 O isotope fractionation during O 2 invasion ( 18 ε I ), the 18 O/ 16 O isotope delta of O 2 in water saturated with atmospheric O 2 ( 18 δ sat ) and the triple isotope fractionation coefficient for O 2 invasion (θ).A similar dual-delta approach can also be used to calculate gross production below the mixed layer, see Eq. ( 42).Prokopenko et al. (2011) also found the approximated calculation of g according to Eq. ( 1) to be deficient.If 17 values are used in the calculation, adoption of a consistent 17 definition is crucial for the accuracy of the calculation.In the past, 17 O excess values based on different 17 definitions have often been mixed indiscriminately, which leads to systematic errors.It is not important which 17 definition is chosen -any single one will do, provided all calculations are performed consistently.This means essentially keeping track of at least one isotope delta (e.g. 18δ) in addition to 17 and a clear definition of 17 whenever it is used.When an approximated calculation of g based on 17 values and Eq. ( 1) is to be performed, a linear definition of 17 according to Eq. ( 4) with κ = γ R should be chosen to minimise errors.However, since the mass-spectrometric methods used to determine the oxygen triple isotope composition of dissolved O 2 , actually yield 17 δ and 18 δ rather than the derived quantity 17 , the dual-delta calculation method described here seems to be the preferable approach.
A marginal aspect to the 17 definition is the question about the symbol for the frequently used value of 10 −6 .The name "per meg" has been used in the past, but is inconsistent with international metrological practice.I suggest using the internationally recognised symbol "ppm" (short for "parts per million") instead.
Table 2 gives a summary of the input parameters and their uncertainties used in the calculation of g.These uncertainties are based on the best single study to have measured the corresponding parameter.Using a set of synthetic 17 δ and 18 δ measurements, I have evaluated the systematic uncertainty introduced into g by the different input parameters.The good news is that for the oceanographically most relevant range, the achievable measurement quality is sufficient to keep the error in g at or below the minimum error estimate of 15 % for the gas exchange coefficient k, which also enters into the calculation of the gross oxygen production according to P = k/z mix c sat g.However, in some cases, other measurements of input parameters exist that are of similar quality to the best single study, but that are not compatible with the corresponding uncertainty bands (Fig. 2).For example, the discussion about the correct 17 sat value is on-going (see Table 1).In case of 17 P , little experimental detail was provided with the only reported measurement (Luz and Barkan, 2000) and other measurements from the same group give a significantly lower value (Sect.5.2).Both 17 sat and 17 P should be re-measured independently.
Considerable differences can arise from using different input parameters and 17 O excess definitions, as shown by the evaluation of g values based on different calculation methods (Table 3, Fig. 3).In the absence of an accepted recommendation on which input parameters and calculation method to use, it will be best to archive the isotope delta values themselves ( 17 δ and 18 δ), together with the oxygen supersaturation s and/or the biological oxygen supersaturation s bio , so that future methodological improvements can be applied retrospectively to existing measurements.One such improvement may be the inclusion of bubble transfer processes in the calculation of g, as shown for bubble injection in Eq. ( 41).
Another advantage of the calculation method suggested here is that it is independent of the assumption of concentration steady state (Eq.48).Moreover, an estimate of the net to gross oxygen production ratio f may be derived from isotope measurements and the oxygen supersaturation s alone (Eq.50).This may be used to check the method for internal consistency and potentially to derive improved estimates of the input parameters based on concomitant measurements of 17 δ, 18 δ and s.

Fig. 2 .
Fig. 2.Relative deviation of g from the base case (see Table2) for different parameters in Eq. (48).Panel (a) corresponds to g = 0.4 and a range of f from −1.0 to +1.0 (negative values correspond to net heterotrophy, positive value to net autotrophy).Panel (b) corresponds to f = 0.1 and range of g from 0.01 to 10 (logarithmic axis).

Fig. 3 .
Fig. 3.Relative deviation of g from the base case for different calculation methods (see Table3).Panel (a) corresponds to g = 0.4 and a range of f from −1.0 to +1.0 (negative values correspond to net heterotrophy, positive value to net autotrophy).Panel (b) corresponds to f = 0.1 and range of g from 0.01 to 10 (logarithmic axis).Black curves correspond to calculation methods based on Eq. (1).Red curves correspond to iterative methods.
O/ 16 O) to fully characterise the corresponding δ value, which is impracticable for more lengthy mathematical expressions.

Table 1 .
Triple oxygen isotope composition of dissolved O 2 at saturation with atmospheric air reported in the literature.Most samples were analysed by headspace equilibration, except row 3b (analysed by membrane extraction). 18δ(B) is 18 δ at saturation according to have been used for the "implied" 17 O excess.
, more suitable λ values of 0.516 and 0.518 were chosen, which are more suitable for