Technical Note : The effect of vertical turbulent mixing on gross O 2 production assessments by the triple isotopic composition of dissolved O 2

The 17 O excess ( 17 Δ) of dissolved O 2 has been used, for over a decade, to estimate gross O 2 production (G 17 OP) rates in the mixed layer (ML) in many regions of the ocean. This estimate relies on a steady-state balance of O 2 fluxes, which include air–sea gas exchange, photosynthesis and respiration but notably, not turbulent mixing with O 2 from the thermocline. In light of recent publications, which showed that neglecting the turbulent flux of O 2 from the thermocline may lead to inaccurate G 17 OP estimations, we present a simple correction for the effect of this flux on ML G 17 OP. The correction is based on a turbulent-flux term between the thermocline and the ML, and use the difference between the ML 17 Δ and that of a single data-point below the ML base. Using a numerical model and measured data we compared turbulence-corrected G 17 OP rates to those calculated without it, and tested the sensitivity of the GOP correction for turbulent flux of O 2 from the thermocline to several parameters. The main source of uncertainty on the correction is the eddy-diffusivity coefficient, which induces an uncertainty of ∼50%. The corrected G 17 OP rates were 10–90% lower than the previously published uncorrected rates, which implies that a large fraction of the photosynthetic O 2 in the ML is actually produced in the thermocline.


Introduction
Gross O 2 production (GOP) in the ocean is a fundamental process in the global cycling of O 2 .As such, accurate estimates of GOP rates are essential in order to understand and model the global cycles of oxygen and carbon.During the past decade, the application of the triple isotope composition ( 16 O, 17 O, and 18 O) of dissolved O 2 as a tracer for GOP (G 17 OP), which was first presented by Luz and Barkan (2000), has become widespread and has been used to estimate GOP rates in many regions of the ocean (Juranek and Quay, 2013, and refs. within).

GOP from 17
Estimating GOP rates from the isotopic composition of dissolved O 2 is based on the 17 O-excess ( 17) that photosynthetically produced O 2 has in comparison to atmospheric O 2 (Luz et al., 1999).The 17 has been defined in several ways (Kaiser, 2011).A common definition (Miller, 2002;Luz and Barkan, 2005), which we use here is λ is the slope of a reference line on a ln(δ 17 O+1) versus ln(δ 18 O+1) plot, which represents the expected slope of the relevant processes.Following Luz and Barkan (2005), most studies use λ = 0.518 for the calculation of 17 and atmospheric O 2 as a standard for the isotopic measurements.To derive a steady-state expression for GOP in the mixed layer (ML), Luz and Barkan (2000) used an O 2 and 17 1-box model.Their derivation yielded the following equation: where K is the air-sea gas-exchange coefficient, (O 2 ) eq the equilibrium concentration of O 2 with the atmosphere, 17 dis E. Wurgaft et al.: The effect of vertical turbulent mixing on gross O 2 production assessments the 17 value of dissolved O 2 , 17 eq the equilibrium 17  with respect to atmospheric O 2 , and 17 p at steady state between photosynthesis and respiration.Recently, Luz and Barkan (2000) method for ML GOP estimation (hereafter G 17 OP LB ) was revised by Prokopenko et al. (2011) and Kaiser (2011), who derived equations for G 17 OP that use measured δ 17 O and δ 18 O, the isotopic composition of dissolved O 2 in air-seawater equilibrium (δ 17 O eq and δ 18 O eq ), and the isotopic composition of photosynthetic O 2 (δ 17 O p and δ 18 O p ). Unlike Eq. ( 2), their equations avoided a number of numerical approximations.They did, however, rely on δ 17 O p and δ 17 O eq which are still subject to disagreement among researchers in the field (e.g.Kaiser and Abe, 2012).However, Luz and Barkan (2011a, b) and Nicholson (2011) showed that if proper δ 17 O p and δ 18 O p are assigned, the differences between G 17 OP LB and the G 17 OP estimated by the revised versions are small.

The effect of turbulent mixing on G 17 OP estimation
As in Luz and Barkan (2000), the ML G 17 OP equations that were presented by Prokopenko et al. (2011) and Kaiser (2011) were derived with a 1-box representation of the ML in which a steady-state balance exists between the O 2 fluxes of GOP, respiration and air-sea gas exchange, but without accounting for the flux of turbulent mixing with O 2 from the thermocline.This was, in spite of the fact that vertical 17 profiles often show a pronounced increase below the ML base (Fig. 1; Luz and Barkan, 2000;Juranek and Quay, 2005;Quay et al., 2010).Juranek and Quay (2005) estimated that vertical turbulence had a negligible affect over G 17 OP.However, Nicholson et al. (2012) showed that mixing of ML O 2 with high 17 O 2 from the thermocline into the ML (by either entrainment due to ML deepening, or by turbulent flux) may result in an overestimation of up to 80 % in ML G 17 OP.Nicholson et al. (2012) further suggested that this overestimation was the likely explanation for the higher ratio of G 17 OP to 14 C-based net primary productivity (G 17 OP : NPP( 14 C); Marra, 2002;Juranek and Quay, 2005;Quay et al., 2010), compared to the ratio of GOP estimated from 18 O incubations to the same net productivity estimate (G 18 OP : NPP( 14 C)).In addition, Jonsson et al. (2013) found that the turbulent mixing had a considerable effect on estimation of net O 2 production, using O 2 : Ar measurements.
As noted above, high 17 O 2 from the thermocline can mix into the ML either by entrainment of water from the thermocline into the ML, which takes place as the ML base deepens, or by vertical turbulence of O 2 from the thermocline.The latter process dominates when the ML depth is constant (Fig. 1).Nicholson et al. (2012) did not consider these two process separately; however, their findings showed that G 17 OP overestimated GOP even when the ML depth was relatively constant, which indicates that vertical turbulence also affects G 17 OP.

. (a)
A schematic illustration of a typical mid-ocean 17 vertical profile.The dashed lines, which extend below the mixed-layer base define the 17 gradient, which in turn, is correlated with the 17 "flux" into the mixed layer.Depending on the profile shape, the choice of a "thermocline" reference point at some depth below the mixed-layer depth, results in a 17 gradient which is different than the real 17 gradient.(b) A conceptual model of the O 2 isotopologues fluxes in and out of the ML.
While Sarma et al. (2006) and Quay et al. (2010) discussed the effect of entrainment on G 17 OP and suggested nonsteady-state corrections, and Castro-Morales et al. (2012) suggested a correction for the vertical flux of O 2 into the ML, the effect of turbulence on the triple isotopic composition and the resulting effect on ML G 17 OP estimations has not been explicitly examined.In light of Nicholson et al. (2012) and Jonsson et al. (2013) results, it is clear that to accurately estimate G 17 OP rates, the magnitude of the turbulent mixing effect should be evaluated, and if large, corrected for.The aims of this work were to evaluate the magnitude of the effect of turbulent mixing on ML G 17 OP estimation, and to derive an analytical correction for this effect.
2 Derivation of a GOP equation with a turbulent mixing term Prokopenko et al. (2011) presented a new equation for G 17 OP (G 17 OP PRO ), which was obtained by rigorous derivation of time variations of O 2 isotopologues.We added a turbulent-flux term to Eq. ( 4) and Eq. ( 5) in Prokopenko et al. (2011).The turbulent flux was calculated between the base of the ML and a single point along the 17 gradient below the ML, which was assigned as "thermocline" (Fig. 1).Consequently, the 17 gradient below the ML was assumed to be linear with depth (we will revisit this assumption in sensitivity tests section).The resulting equation for the rate of change in 17 in the ML was (full derivation of the equation can be found in Appendix A): (3) where h is the ML depth, (O 2 ) is the dissolved O 2 concentration in the ML, κ is the eddy-diffusivity coefficient, and Z is the vertical distance between the base of the ML and the depth assigned as "thermocline".For convenience, we use D = κ/Z hereafter.X * represents the ratio * O/ 16 O.The subscripts "p" and "eq" denote "photosynthetic" and "equilibrium", respectively.Note that as was shown by Luz and Barkan (2009) and Prokopenko et al. (2011), the rate of change of 17 is independent of respiration.When steadystate conditions in the ML are assumed, the resulting term for turbulent-flux corrected G 17 OP is where G 17 OP C is the G 17 OP corrected for turbulent flux of O 2 from the thermocline.For convenience we will abbreviate the GOP correction for turbulent flux of O 2 from the thermocline (the second term on the right-hand side in Eq. 4) to "TFC" hereafter.The numerator of the TFC represents the contribution of the turbulent flux to the GOP estimated from 17 in the ML.mixing (model equations, parameters and MATLAB files are given in the Supplement).An additional layer at the top of the water column represented the ocean surface.In this layer, (O 2 ) and its isotopic composition were kept in airsea equilibrium.The eddy-diffusivity coefficient in the ML (2.5 × 10 −3 m 2 s −1 ) was assigned so as to let (O 2 ) in the ML be fully mixed.The concentration of each isotopologue in the ML was affected by turbulent mixing with the uppermost layer, photosynthesis, respiration and turbulent mixing with the seasonal thermocline residing below.For the seasonal thermocline, we used κ = 10 −4 m 2 s −1 (see sensitivity tests section).
We ran two simulations to examine the effect of turbulent mixing on ML G 17 OP.In the "fixed mixing depth" simulation (Table 1) we ran the model with a constant ML depth of 40 m.The layer directly below the ML base, at 50 m, was assigned as the "thermocline" data point for TFC calculation.ML δ 17 O, δ 18 O and 17 values were calculated every 30 model time steps (model month).In the first model month of the simulation, G 17 OP LB and G 17 OP PRO slightly overestimated (by 7 %) the GOP assigned in the model (GOP M ).
In the following model months, the overestimation of both G 17 OP LB and G 17 OP PRO increased, reaching ∼ 40 % after 5 model months.As shown in Fig. 2, the increase in GOP overestimation was closely related to the increase in 17 in the seasonal thermocline.However, G 17 OP C , which corrects for the turbulent mixing flux of O 2 from the thermocline, remained constant with a slight underestimation of ∼ 7 % throughout the entire simulation period.When we effectively shut down turbulent mixing in the model by reducing κ within the thermocline to 10 −6 m 2 s −1 , G 17 OP LB and G 17 OP PRO were in good agreement with GOP M .This indicates that turbulent flux was indeed the cause of the overestimation.
In the "varying mixing depth" simulation, ML depths were allowed to change, roughly according to the typical monthly variations in BATS station ( All GOP rates are in mmol m −2 day −1 .In the columns marked with t, κ = 10 −4 m 2 s −1 , whereas in the columns marked with nt, κ = 10 −6 m 2 s −1 , which effectively shuts down the turbulent flux of O 2 between the thermocline and the mixed layer.

Sensitivity tests
In addition to simulating the effect of turbulent mixing on GOP estimations, we used the 1-D model to test the sensitivity of the TFC to the depth of the "thermocline" reference point, to the analytical error associated with 17 measurements, and to the ML depth.The choice of the depth which represents the "thermocline" point can affect the resulting GOP (Fig. 1).While the TFC assumes a linear 17 gradient between the ML and the "thermocline", the actual 17 gradient is not necessarily so.Therefore, the closer the "thermocline" point to the ML base, the better it represents the actual fluxes of O 2 isotopolouges between the ML and the thermocline (Fig. 1).On the other hand, the difference in the isotopic composition between the ML and the "thermocline" has to be considerably larger than the analytical error associated with 17 measurements (∼ 7 per meg; e.g.Reuer et al., 2007).In order to assess the sensitivity of the TFC to the depth of the "thermocline" reference point and to the analytical error on 17 , we calculated the magnitude of the turbulence correction with different "thermocline" reference points.As illustrated in Fig. 3a, the magnitude of the TFC decreased as the distance between the ML base and the "thermocline" reference point increased.Assuming that the TFC calculated using the "thermocline" immediately (10 m) below the ML is the most accurate, we consider the difference between this value and the TFC calculated for other "thermocline" depths as the error induced on the TFC by the selection of the "thermocline" depth.Using the same 1-D model simulations, we calculated the effect of the analytical error associated with 17 measurements on the magnitude of the TFC.For this end, we calculated the TFC three times per each "thermocline" depth.The first TFC was calculated without error on 17 values, whereas for the other two, a 7 per meg error was either added or subtracted from the 17 values of the ML and the "thermocline".The combinations which yielded the maximal deviations from the noerror TFC are illustrated by the vertical error bars in Fig. 3a.The magnitude of this error decreased with increasing depth of the "thermocline".Finally, we estimated the combined effect of the choice of the "thermocline" depth and the analytical error on the TFC, by propagating these two errors.The resulting error is illustrated in Fig. 3a.For the scenario used in our tests, the minimal error (∼ 20 %) was observed at 20 m below the ML base.We note that this error is likely to be an overestimation of the actual error which results from these two factors, since we used the maximal analytical error on both ML 17 and 17 thr values simultaneously, and in reverse directions, which yielded the maximal effect on the resulting 17 gradient.Moreover, as the ML 17 value is usually an average of several measurements, it is likely to be subject to a much smaller error than 7 per meg.
The magnitude of the TFC is linearly dependent on the value of κ used in Eq. ( 4).However, our simulations also showed that using κ values smaller than 1 × 10 −4 m 2 s −1 , yielded 17 values in the thermocline that were higher than 200 per meg, and when we used κ = 0.5 × 10 −4 m 2 s −1 17 Table 2. Results of model simulation of GOP with varying mixed-layer depth ML -mixed layer, GOP-M -input GOP rate, G 17 OP LB - Luz and Barkan (2000), G 17 OP PRO -Prokopenko et al. (2011), G 17 OP C -this work includes a correction for turbulent flux of O 2 from the thermocline.All GOP rates are in mmol m −2 day −1 .In the columns marked with t, κ = 10 −4 m 2 s −1 , whereas in the columns marked with nt, κ = 10 −6 m 2 s −1 , which effectively shuts down the turbulent flux of O 2 between the thermocline and the mixed layer.2012) used similar values (8-9 × 10 −5 m 2 s −1 ) to reproduce the physical conditions (ML depth, heat content and sea surface temperature) in the upper 1000 m in the BATS and in HOT.Therefore, we estimate that our choice of κ was rather accurate for the processes of turbulent mixing between the ML and the seasonal thermocline.Apparently, in spite of the fact that 1 × 10 −4 m 2 s −1 is almost an order of magnitude higher than the value estimated from SF 6 release experiments in the permanent thermocline (∼ 300 m, Ledwell et al., 1993), κ values near the interface between the ML and the thermocline are higher than those which characterize the thermocline at greater depths where SF 6 release experiments were conducted.Given the highly unrealistic profiles we obtained for lower values and the agreement between our model and the model used by Nicholson et al. (2012), we assume an uncertainty of ∼ 50 % on this value, and consequently, a 50 % uncertainty on the TFC term.
Finally, we tested the sensitivity of the TFC to the ML depth.The results (Fig. 3b) show that as the ML depth increases, the contribution of turbulence to the uncorrected GOP decreases.This implies that correcting G 17 OP rates for turbulent fluxes of O 2 isotopologues from the thermocline is especially important in ocean areas in which the summer ML is relatively shallow, such as in BATS.

The effect of turbulent mixing on measured GOP rates
To estimate the effect of turbulent mixing on G 17 OP estimations in the ocean, we compared previously published GOP rates from BATS (Luz and Barkan, 2009) with equivalent G 17 OP C rates.In addition, we used published data (Nicholson et al., 2012) (Nicholson et al., 2012).Following our conclusions from the 1-D model simulations, we compared months in which there were no major changes in ML depth.The κ values for BATS and HOT were both taken from Nicholson et al. (2012).ML depth was determined as the depth in which a difference of 0.5 % in O 2 concentration relative to the sea surface was observed (Castro-Morales and Kaiser, 2012).The "thermocline" reference point was assigned as the depth nearest to the ML base in which 17 was at least 14 per meg larger than 17 in the ML.In BATS, K values were obtained from Luz and Barkan (2009).In HOT, wind speed data from ∼ 10 days before the cruise were obtained from QuickScat database, and K was estimated according to Sweeney et al. (2007).The results (Table 3) showed that G 17 OP C rates were 65-100 % lower than G 17 OP LB in BATS, and 10-40 % lower in HOT.These results are in agreement with Nicholson et al. (2012), who estimated the effect of mixing as 60-90 % of G 17 OP.We note that Luz and Barkan (2009) used Wanninkhof (1992) parameterization for K in BATS.Applying more recent parameterization (Ho et al., 2006;Sweeny et al., 2007), which gives lower estimates of gas-exchange rates, would result in an even greater contribution of the O 2 turbulent flux from the thermocline.

Discussion
Our results showed that vertical turbulence of O 2 from the thermocline, affects ML G 17 OP estimations, such that corrected G 17 OP rates are considerably lower than the uncorrected rates.This indicates that a large fraction of the ML photosynthetic O 2 is not produced in the ML itself, but rather in the thermocline below it.However, our results also show that the turbulent contribution to the ML 17 can be corrected in a rather simple manner.While the correction itself was derived using several approximations, such as constant ML depth, and constant 17 gradient between the thermocline and the ML, most of the resulting uncertainties also exist in the uncorrected equations for G 17 OP (Luz andBarkan 2000, Prokopenko et al., 2010).These uncertainties are intrinsic to any extrapolation in time and space of GOP rates estimations from a "snap-shot" profile.Moreover, our sensitivity tests showed that the uncertainty on the TFC is smaller than 50 %, and therefore, we concluded that in spite of the approximations and uncertainties involved, using the TFC will improve the accuracy of G 17 OP estimations.Below, we discuss the technical aspects and the implications of these findings.

Applying GOP correction for turbulent flux of O 2 from the thermocline in oceanographic measurements
Our results showed that when the ML depth does not change considerably, the effect of turbulent flux on ML GOP estimation can be corrected.However, to apply the correction, at least one point of data, which includes [O 2 ] and its isotopic composition below the ML is necessary.Such data points are easy to obtain in time-series study sites, such as BATS and HOT, but can complicate basin-wide G 17 OP estimations, which usually rely on underway seawater systems installed on ships of opportunity for sampling (Juranek and Quay, 2010;Juranek et al., 2012).In the future, such studies would need to either collect several representative "thermocline" samples from the thermocline along the cruise route, or use existing data if available to apply the TFC.Such corrections could also be applied to existing data.
The TFC is sensitive to the exact depth of the "thermocline" data point (Fig. 3a), and to the analytical error on 17 measurements.Technical improvements which would reduce the error on 17 measurements would also increase the accuracy of G 17 OP estimations in general, and the accuracy of the TFC in particular.As shown in Fig. 3a, the error induced on the TFC by the analytical error, decreases with the depth selected for the "thermocline", while the error induced by deviation of the 17 profile from a linear one increases with depth.In the depth range immediately below the ML base, these effects cancel each other to some extent, and an optimal depth can be estimated.This depth is dependent on the vertical distribution of 17 in each study-site, and on the analytical error associated with the 17 measurement process applied.In practical terms, we suggest using the minimum depth in which the 17 value is significantly (within the analytical uncertainty) different than the ML value.
Previous knowledge of 17 dynamics for the study site would help with choosing the optimal depth to collect the "thermocline" sample.For example, Juranek and Quay (2005) and Quay et al. (2010) have shown that during summer months in HOT, differences greater than 60 per meg (an order of magnitude higher than the analytical error) between the ML and thermocline 17 could be found within 20-40 m below the ML base.In BATS, 17 gradients tend to be smaller, and a difference of 60 per meg can usually be found 40-60 m below the ML depth (Nicholson et al., 2012).Therefore, the "thermocline" optimal depth is likely to be shallower in HOT than in BATS, and consequently, associated with a smaller error.

Parameterization of gas exchange and eddy diffusivity
The TFC, and consequently G 17 OP C are sensitive to the accuracy of both K and κ.While the sensitivity to K characterizes G 17 OP in general, a combination of errors on K and  -Luz and Barkan (2000), "thermocline" -depth of the data point representing the thermocline.G 17 OP C -G 17 OP corrected for turbulent flux of O 2 from the thermocline.All the parameters used for the calculations made in BATS were taken from Luz and Barkan (2009).κ values and the raw data used for the calculation of G 17 OP LB and G 17 OP C , and the parameters used for the calculations in HOT were taken from Nicholson et al. (2012).ML depth was determined as the depth in which a difference of 0.5 % in O 2 concentration relative to the sea surface was observed.Gas-exchange rates in HOT were calculated according to Sweeney et al. (2007)  κ may yield inaccurate G 17 OP C rates.We assume that the negative G 17 OP C rate that we obtained in BATS in September 2000 (Table 3) was the result of inaccurate choice of K and κ.We also acknowledge the fact that the uncertainty on κ is the largest source of error on the magnitude of the TFC.Since the main aims of this work were to show the importance of turbulence effects and to suggest a correction, rather than to perform accurate GOP estimations, we used crude estimations of K and κ.However, the fact that 17 in the ML depends upon GOP, gas-exchange and turbulent flux from the thermocline, means that in future studies any one of these three parameters could be estimated by performing simultaneous measurements of the other two parameters.For example, if GOP is estimated by 18 O incubations and gas exchange is estimated from wind speed measurements, 17 profiles could be used to estimate κ in the seasonal thermocline.Moreover, since 18 O incubations are not affected by turbulence, it is likely that provided that K and κ are accurately parameterized, G 17 OP C : N 14 CP would agree with G 18 OP : N 14 CP (Marra, 2002).

Conclusions
1. Turbulent fluxes of O 2 isotopologues from the thermocline have a pronounced effect over 17 values in the ML, and consequently, over the accuracy of G 17 OP estimations.
2. An accurate G 17 OP estimate can be obtained by using a simple correction for the effect of the turbulent fluxes.
3. The main source of uncertainty on the GOP correction for turbulent flux of O 2 from the thermocline is the eddy-diffusivity coefficient, which causes ∼ 50 % uncertainty.
4. The GOP correction for turbulent flux of O 2 from the thermocline is applicable when the mixed-layer depth does not change sharply, and requires measurements of O 2 and its isotopic composition in a single point below the mixed layer, in addition to the standard measurements of these values in the mixed layer.where h is the mixed-layer depth, t is time, GOP is the gross O 2 production, R is the respiration rate and K is the piston velocity.α * is the fractionation factor associated with respiration for each isotopologue.The subscripts "p" and "eq" denote "photosynthetic" and "equilibrium", respectively.The remainder of the derivation is carried out by straightforward applications of the steps outlined in Prokopenko et al. (2011), and repeated here for the sake of completion.The left-hand side of Eq. (A4) can be written explicitly as Upon substituting the left-hand side, and the first term on the right-hand side of Eq. (A3) with Eq. (A4) and Eq.(A3), respectively, rearranging and dividing by X * one gets
photosynthesis respiration turbulent mixing with the thermocline

Fig. 3 .
Fig. 3. Sensitivity of the turbulent-flux correction.(a) The sensitivity of the GOP correction for turbulent flux of O 2 from the thermocline (TFC) to the depth of the "thermocline" reference point below the mixed layer (ML).The error bars represent the maximal error on the TFC, induced by the analytical error associated with 17 measurements (7 per meg).The error resulting from the combination of these two factors (circles) on the TFC shows a minimum at 20 m below ML.(b) The relative contribution of turbulent mixing flux to estimated G17OP as a function of the ML depth.Model conditions for each sensitivity test are described in Sect. 4 in the text.

Table 1 .
-Prokopenko et al. (2011)on with a constant mixed-layer depth of 40 m GOP-M -input GOP rate, G 17 OP LB -Luz and Barkan (2000), G 17 OP PRO-Prokopenko et al. (2011), G 17 OP C -this work includes a correction for turbulent flux of O 2 from the thermocline.
17OP C in Hawaii Ocean time series (HOT).To use the TFC, we chose months in which vertical profiles of 17 and O 2 were measured and published to calculate and compare between G 17 OP LB www.biogeosciences.net/10/8363/2013/Biogeosciences,10, 8363-8371, 2013and G

Table 3 .
Comparison of G 17 OP rates from BATS and HOT ML -mixed layer, G 17 OP LB from wind speed data obtained from QuickSCAT.