Oxygen corrections applied to delayed-mode data

This notebook performs analysis and correction of oxygen data from the following C-PROOF glider deployment:

************
** Calvert Island Line:  glider dfo-bb046**
************
* Deployment: dfo-bb046-20220707

* Optode: Rinko
* Using oxy sensor, with serial number: 0054
* Calibration date: 2021-11-01 00:00:00
* Using Milo filter

* Processing steps will be saved in: Oxygen_dfo-bb046-20220707.html
* Processed by AH, IOS Data Products team, Ocean Sciences Division, Fisheries and Oceans Canada
* Processing date: 20250109

1.0 Preamble

This document describes dissolved oxygen data processing steps applied to delayed mode data collected using a Rinko ARO-FT-CE dissolved oxygen optode mounted on an Alseamar SeaExplorer autonomous ocean glider. This sensor is mounted in the nose cone of the glider, adjacent to the CTD, and measures oxygen concentration. It has a response time of less than 1 s. Corrections applied use data collected by the CTD on the glider. Note that when converting between oxygen units, potential density is used.

1.1 Set up the processing

The processing steps below are applied to delayed mode data stored in a single netcdf timeseries file created using the Pyglider data processing package (https://github.com/c-proof/pyglider). This data has already undergone CTD sensor corrections.

The metadata and a detailed report describing the corrections being applied will be made available via the deployment page on the C-PROOF website: https://cproof.uvic.ca/gliderdata/deployments/dfo-bb046/dfo-bb046-20220707/

The autoreload extension is already loaded. To reload it, use:
  %reload_ext autoreload

Loading coefficients for Rinko sensor, SN 0054

2.0 Corrections applied to delayed mode data for this mission

This document covers the following steps:

1. Compensate for pressure and salinity effects
2. Convert from concentration in $\mu$mol/L to phase
3. Determine and apply the sensor response time lag correction
4. Convert from phase to concentration in $\mu$mol/L
5. Determine and apply the optode drift correction

2.1 Compensate for pressure and salinity effects

The Rinko ARO-FT User's manual states that corrections must be applied for pressure and salinity effects to the dissolved oxygen concentration in 𝜇mol/L (oxygen_concentration).

The pressure compensation factor is estimated as: $p_factor = (|P| / 100) ⋅ Cp + 1$

where P is the pressure from the CTD in dbars, and Cp is the calibration coefficient for pressure correction, listed on the calibration certificate for each sensor.

The salinity compensation factor is estimated as: $s_{}factor = e^{[S ⋅ (B0 + B1 ⋅ Ts + B2 ⋅ Ts^{2} + B3 ⋅ Ts^{3}) + C0 ⋅ S2]}$

where S is the fully corrected salinity from the CTD in PSU, B0 to B3 and C0 are coefficients specified in the Rinko manual, Ts is the scaled temperature: $Ts = ln[(298.15 - T) / (273.15 + T)]$; and T is the temperature from the CTD in degrees Celsius.

Unlike the Aanderaa sensors, the Rinko optodes apparently do not have an internal reference salinity. The manual states “Note that the physical value output of the instrument does not consider salinity and/or pressure.” This is equivalent to sal_ref = 0 for the purposes of the salinity compensation correction above.

Then compensated oxygen concentration is calculated: $oxygen_{}comp = oxygen_{}uncomp ⋅ s_{}factor ⋅ p_{}factor$

The units of oxygen_comp are 𝜇mol/L.

Top left: Some of the oxygen profiles, without corrections. Top right: Difference between compensated and uncompensated oxygen for the first part of the mission. Bottom left: uncompensated and compensated (red) oxygen by depth (pressure). Bottom right: Salinity vs oxygen, showing compensated and uncompensated values for the first part of the glider mission.

Mean percent difference: -18.39 +/- 0.96
Left: Histogram of oxygen values before and after compensation. Right: Histogram of the change in oxygen before and after compensation, for the first part of the mission.

2.2 Convert oxygen concentration in $\mu$mol L$^{-1}$ to phase

The uncalibrated phase (Phase) is the difference between the phase obtained when the foil is excited by Blue light and by Red light. These hexadecimal values are referred to as ‘AD values’, as opposed to calculated temperature and dissolved oxygen concentration, which are referred to as ‘Physical values’.

The time response correction (5.1) is applied to the AD value of dissolved oxygen (AD value of DO), which is functionally equivalent to the calibrated phase. Then the dissolved oxygen concentration in 𝜇mol/L is calculated using the Stern-Volmer relationship proposed by Uchida et al., 2010:

$$ DO [{\mu}mol/L] = {((1 + d0 ⋅ T)/(d1 + d2 ⋅ N + d3 ⋅ t + d4 ⋅ t ⋅ N))^{e0} - 1} ⋅ {1/(c0 + c1⋅T + c2⋅T^2)} $$

where calibration coefficients $c0$, $d0$, $e0$ etc. are printed on the calibration certificate for that sensor, $T$ is temperature, and $N$ is the AD value of $DO/10000$. Typically $e0 = 1$ and $d3 = d4 = 0$, so that we do not require values of t, the LED accumulation time [s].

Our sensors are normally set to output dissolved oxygen concentration in 𝜇mol/L (oxygen_concentration), which is then converted to calibrated phase (AD value of DO) using the temperature measurement from the oxygen sensor, the time response correction is applied, and then the resulting field is converted back to DO[𝜇mol/L] using the fully corrected temperature measurement from the CTD.

Oxygen in phase (left), after conversion from concentration (right):

2.3 Sensor reponse time lag correction

Determine the response time lag:

The time response correction is applied to the calibrated optode phase (AD value of DO for Rinko sensors) to account for the slow sensor response time relative to the diving rate of the glider.

Following Nicholson et al., 2008 and Howatt et al., 2018, a response time constant $\tau$ is determined by minimizing the error between all measured upcast values vs. all measured downcast values. A first-order lag correction is determined by calculating $\tau$ as in Fofonoff et al., 1974 based on the equation:

$$Phase_{true} = Phase_{measured} + \tau ⋅ d/dt ⋅ Phase_{measured}$$

where $Phase$ is the calibrated optode phase, ‘true’ is the true value and ‘measured’ is the measured value, and $d/dt$ indicates a first derivative w.r.t. time during a given upcast or downcast. Derivatives are calculated using a centered finite difference scheme after successive upcasts and downcasts are interpolated onto the same depth grid. To avoid large derivatives caused by points at the start and end of individual profiles, any $dt$ value greater than one standard deviation from the mean is excluded from the calculation.

Then a least squares linear regression is performed on the following equation: $$Phase_{measured}^{upcast} -Phase_{measured}^{downcast} = - \tau ⋅ (d/dt ⋅ Phase_{measured}^{upcast} - d/dt ⋅ Phase_{measured}^{downcast})$$

which assumes that successive upcasts and downcasts are identical.

Expected ranges for the response time constant for each type of sensor are provided in Bittig et al., 2014. For our Rinko sensors, the calculated value of $\tau$ is expected to fall in the range 1 to 3 seconds. If the calculated response time constant falls outside of these ranges, manual inspection of the mission data is required to determine the cause.

Given the rapid response time of Rinko sensors, it is expected that $\tau$ may be too small to necessitate applying a time lag correction. We verify this below.

Grouping upcasts and downcasts, stand by.
Grouping upcasts and downcasts complete.

Interpolating fields.

100%|████████████████████████████████████████████████████████████████████████████████| 445/445 [00:07<00:00, 61.31it/s]

Now we use the upcasts and downcasts that have been interpolated to a regular depth grid to estimate $\tau$ using least-squares regression applied to the entire time series.

****
mean phase diff between up and downcasts: 0.014 +/- 0.073
****
Sensor response lag coefficient [in seconds] tau = 0.447
99th percentile significance = 0.176
R-value for the linear regression = -0.056

Phase difference between up and downcasts, compared with the up-down difference of its first derivative. Phase difference vs d/dt (left), shown with histogram of d/dt values (middle) and difference values (right):

Applying the sensor response time lag correction:

The value of $\tau$ is small as expected, given the 1s response time of the sensor. We note that the least squares fit had a relatively low R-value, suggesting that the upcasts and downcasts were not offset by a consistent amount throughout the time series. Given this, and the small offset, we elect not to apply a sensor reponse time lag correction. This choice additionally avoids applying the smoothing and filter required for this correction, better preserving the signals in the data.

Oxygen concentration for up and downcasts, shown against presssure (left) and adjusted salinity (right):

2.4 Convert from phase to oxygen concentration in $\mu$mol/L

Given that we are not applying a time response correction, we use 'oxygen_comp' (oxygen concentration in $\mu$mol/L) for the next step in the processing.

2.5 Optode drift correction

The sensor drift correction is applied to the time-response corrected, compensated dissolved oxygen concentration data 'DO' to account for the tendency of oxygen optodes to drift over time as the sensor ages. Although drift increases as a function of the number of samples collected, sensors tend to drift more quickly during storage than while in the field by about a factor of 10.

The sensor drift correction attempts to correct for the drift that does occur by comparing the oxygen optode data to a trustworth profile of dissolved oxygen concentration measurements collected nearby in space and time. Ideally, this comparison would be with Winkler titrated bottle samples collected at the same time and location as one of the glider’s dives. In practice, a hierarchy of ‘trustworthy’ data sources is used, within some chosen tolerance for distance and time from the glider’s sampling locations.

The drift correction itself is a linear adjustment applied to DO[$\mu$mol/kg] with slope p1 and intercept p2:

oxygen_concentration_adjusted [$\mu$mol/kg] = p1 DO[$\mu$mol/kg] + p2

Following Bittig et al., 2018, the drift is characterized by a reduced sensor response plus a small offset in the opposite direction. The slope should be positive, and within about 10% of unity. The intercept should be within a few [$\mu$mol/kg] of zero and negative. In practice, the correction is determined with the intercept set to zero a priori. We have found that an intercept of zero typically produces better agreement in the near-surface mixed layer and in deep water, the most appropriate segments of the water column for comparison.

If the slope is too large, the comparison may be rejected. It is expected that the drift for an individual sensor will increase with time (though not linearly). If this is not the case, or if a sudden large increase in drift is detected, the comparison should be re-evaluated and/or the sensor should be examined.

Determining optode drift:

To determine the drift correction, a least-squares fit is performed in density space between all of the suitable comparison profiles and the associated glider optode measurements. The glider profiles to be used are bin-averaged to the densities associated with the oxygen measurements from the comparison profile, with a bin size of +/-0.1 kg/m$^3$ (which can be adjusted higher if needed). The linear least-squares fit is applied to the binned oxygen optode measurements from the glider vs. the measurements from the comparison profile. The associated correlation coefficient for the fit is expected to be above r$^2$ = 0.95, and in practice is often above r$^2$ = 0.99.

Note that potential density is used, not in-situ density, since a given volume of water may move vertically by a significant amount due to isopycnal heaving, and we wish to remove this effect from the comparison. The portion of the profiles within the pycnocline is excluded from the comparison by visually identifying a depth at the base of the pycnocline and retaining only data below it, as well as data in the very near surface that will be close to the saturation value.

The goodness of fit is evaluated by examining the agreement between the deep and shallow drift-corrected glider oxygen measurements with those from the comparison profile. For the shallowest measurements, the saturation value of oxygen at the comparison profile location is calculated and used to provide context for the expected near-surface oxygen concentration values.

Finally, the linear sensor drift correction is applied to the DO measurements (not to the data binned by density).

For this mission, two comparison profiles are available:

1. A profile of Winkler-titrated bottle samples collected by the Hakai Institute at the time of glider recovery in Hakai Pass, and,
2. A profile using an oxygen sensor on a lowered CTD, also collected by the Hakai Institute at the time of glider recovery.

We examine both these profiles compared to the available glider data.

Loading in situ data
The bottom depth for the profile is 340.55 metres
CTD cast and bottle sample data shown with depth for salinity (left), temperature (middle), and oxygen (right):

Temperature vs salinity, and oxygen vs potential density for the bottle and CTD cast samples:

We see good agreement between the comparison profile from the oxygen sensor used during the CTD cast and the Winkler bottle samples, with the exception of the densest, deepest water. For the deepest Winkler sample, it appears that the bottle salinity is quite different from that measured by the CTD, which affects the oxygen agreement in potential density space as well. We now compare to the available glider data.

Using minimum depth of: 250

****
Glider must be within 0.5 km
Profile must be within 1 day(s) of sample
****
Number of glider profiles available for comparison:  1
Profile indices of profile(s) available for comparison:  [893.]

[-128.05360833 -128.05360666 -128.05360498 ... -128.0444552  -128.04445629
 -128.04445737]

Top: latitude and longitude of glider profiles, with the comparison profiles and chosen profile indicated. Bottom left: Oxygen from glider profile (black) shown with CTD (red points) and bottle (red stars) samples. Bottom right: Potential density vs oxgen for the glider profile and sample data from CTD and bottles.

The comparison profile was collected at the time of glider recovery, during the final glider profile, and includes high-quality Winkler-titrated bottle samples as well as the profile collected by the oxygen sensor attached to the CTD. Given the agreement between the Winkler samples and the oxygen from the CTD profile, and given the small number of Winkler samples available for comparison, we determine the correction to the glider using the oxygen sensor data.

We compare to the last full-depth glider profile collected, and determine the comparison using all available data, not restricting the densities. There is a slight offset from the glider oxygen profile visible in density space, which is what we expect to correct.

Using densities between 1020 and 1027

Binned and unbinned oxygen collected by the CTD, plotted against potential density (left). Potential density vs oxygen also shown for the binned and unbinned glider profile (right).

Potential density vs oxygen concentration for the binned glider and CTD profiles (left). Oxygen from the CTD vs oxygen from the glider profile (right).

We now apply a least-squares fit between the two profiles in isopycnal space. We force the intercept to pass through zero, as this has been demonstrated to result in more accurate corrections. As oxygen optodes drift, they tend to read lower values of oxygen [Bittig et al., 2015], so we expect the correction to produce a slight increase in oxygen concentration.

slope m = [1.097] and R2 = 0.97
p: 0.0, stnerr: 0.0096, intercept: 3.518

2.5 Optode drift correction

Sanity checks:

We now apply our drift correction and confirm that the drift correction has had the intended effect. We check the following:

1. Do we see better agreement with the full comparison profile, particularly in the deep water and at the surface?
2. Are values at the surface consistent with saturated (or over-satured) water?
3. Are values at 1000 m consistent with other deep oxygen samples or at least reasonable?

Applied slope of 1.097
Oxygen from the CTD profile vs glider oxygen, shown before (black) and after (red) applying a linear correction:

Glider oxygen before and after the linear correction:

Potential density vs oxygen shown for the glider and sampled oxygen data before correction (left) and after correction applied to the glider profile (right):

Agreement with comparison profile: The agreement with the comparison profile has improved significantly, except near the surface where there is high variability.

Depth vs oxygen concentration for the glider profile (black), CTD profile (red), and oxygen saturation (blue) both before (left) and after (right) the drift correction:

Surface values: Both the glider profile and comparison profile appear over-saturated right at the surface. Before the drift correction, the glider is in decent agreement with the comparison profile below 5 m but less saturated above 5 m. After the correction, agreement below 5 m is not quite as good, but the glider profile is closer to 100% saturation above 5 m and matches to the comparison profile more closely. This type of surface comparison provides information, but is not conclusive, particularly in a dynamic environment such as Hakai Pass.

The minimum oxygen value before any oxygen corrections is 5.0 umol/kg
The minimum oxygen value before the drift correction is 4.22 umol/kg
The minimum oxygen value after the drift correction is 4.51 umol/kg
Glider oxygen shown at three stages: with no corrections applied (grey), directly before the drift correction with the compensations and other corrections applied (light blue), and directly following the linear correction using the sample data (dark blue):

Comparison of deep oxygen: There are no good DFO profiles from beyond the shelf break to compare with during the glider mission. However, we see that the drift correction has increased the deep oxygen values slightly. This is a good sanity check, since oxygen optode sensor drift tends to cause lower values. Furthermore, the change caused by the drift correction is smallest for the deep oxygen data, as expected. However, the oxygen has overall slightly decreased relative to the delayed-mode data, when considering all corrections performed here.

Change in oxygen concentration before and after corrections applied shown over time (left), and by concentration (right).

3.0 Summary of corrections applied

Given the high-quality comparison profile, and the supporting evidence suggesting the estimated drift correction produces reasonable results, we choose to apply the drift correction. However, note that the optode on this glider was calibrated only 8 months before the mission, and the drift correction estimated is large.

The change in oxygen concentration after applying the drift correction is between 0 and 10$\mu$mol/kg below 400m, and up to nearly 40$\mu$mol/kg at the surface.

Corrected data saved to file: e:/Glider/data//dfo-bb046/dfo-bb046-20220707/dfo-bb046-20220707_corrected_v2.nc

Gridfile created: e:/Glider/data//dfo-bb046/dfo-bb046-20220707/dfo-bb046-20220707_grid_corrected_v2.nc

The corrected salinity and temperature, shown with filtered conductivity and corrected oxygen:

The delayed-mode oxygen, corrected oxygen, and difference between the two. The red star indicates the profile used for comparison with sample data:

Overall, the oxygen declined relative to the delayed-mode data when all conversions and adjustments were applied, though we can see a few profiles that increased. These may require further steps to examine further if the data is valid.

References

1. AANDERAA TD 269, ‘OPERATING MANUAL OXYGEN OPTODE 4330, 4831, 4835’. June 2017. Aanderaa Data Instruments AS

2. RINKO A-0133-05, ‘Rinko Fast optical DO sensor for integration RINKO FT (Model name: ARO-FT, AROD-FT) User’s Manual’. JFE Advantech Co., Ltd.

3. Uchida H., Kawano, T., Kaneko, I., and Fukasawa, M., ‘In-situ calibration of optode-based oxygen sensors’. 2008. Journal of Atmospheric and Oceanic Technology, 25: 2271-2281.

4. Uchida, H. Johnson, G.,C. and McTaggart, G,. C., ‘CTD Oxygen Sensor Calibration Procedures’. In, The GO-SHIP Repeat Hydrography Manual: A Collection of Expert Reports and Guidelines. 2010. Version 1, (eds Hood, E.M., C.L. Sabine, and B.M. Sloyan), 17pp. (IOCCP Report Number 14; ICPO Publication Series Number 134). DOI: https://doi.org/10.25607/OBP-1344.

5. Bittig, H. C., Fiedler, B., Scholz, R., Krahmann, G., Körtzinger, A., ‘Time response of oxygen optodes on profiling platforms and its dependence on flow speed and temperature’. 2014. Limnology and Oceanography: Methods, 12. DOI: 10.4319/lom.2014.12.617.

6. Howatt, T., Palter, J. B., Matthews, R., deYoung, B., Bachmayer, R. and Claus, B., ‘Ekman and Eddy Exchange of Freshwater and Oxygen across the Labrador Shelf Break’. 2018. J. Phys. Oceanogr., 48, 1015–1031. DOI: 10.1175/JPO-D-17-0148.1.

7. Nicholson, D., Emerson, S. and Eriksen, C. C., ‘Net community production in the deep euphotic zone of the subtropical North Pacific gyre from glider surveys’. 2008. Limnol. Oceanogr., 53, 2226–2236. DOI: 10.4319/lo.2008.53.5_part_2.2226.

8. Fofonoff, N. P., Hayes, S. P. and Millard Jr., R. C., ‘W.H.O.I./Brown CTD microprofiler: Methods of calibration and data handling’. 1974. Woods Hole Oceanographic Institution Tech. Rep. WHOI-74-89, 64 pp. DOI: 10.1575/1912/647.

9. Miloshevich, L. M., Paukkunen, A., Vömel, H., and Oltmans, S. J., ‘Development and Validation of a Time-Lag Correction for Vaisala Radiosonde Humidity Measurements’. 2004. J. Atmos. Oceanic Tech., 21, 1305–1327. DOI: 10.1175/1520-0426(2004)021<1305:DAVOAT>2.0.CO;2.

10. Bittig, H. C., Körtzinger, A., Neill, C., Van Ooijen, E., Plant, J. N., Hahn, J., Johnson, K. S., Yang, B. and Emerson, S. R., ‘Oxygen optode sensors: principle, characterization, calibration, and application in the ocean’. 2019. Frontiers in Marine Science, 4, p.429.