A turbulence data reduction scheme for autonomous and expendable profiling floats

Autonomous and expendable profiling-float arrays such as those deployed in the Argo Program require the transmission of reliable data from remote sites. However, existing satellite data transfer rates preclude complete transmission of rapidly sampled turbulence measurements. It is therefore necessary to reduce turbulence data on board. Here we propose a scheme for onboard data reduction and test it with existing turbulence data obtained with a modified SOLO-II profiling float. First, voltage spectra are derived from shear probe and fast-thermistor signals. Then, we focus on a fixed-frequency band that we know to be unaffected by vibrations and that approximately corresponds to a wavenumber band of 5–25 cpm. Over the fixed-frequency band, we make simple power law fits that – after calibration and correction in post-processing – yield values for the turbulent kinetic energy dissipation rate ϵ and thermal-variance dissipation rate χ. With roughly 1 m vertical segments, this scheme reduces the necessary data transfer volume 300-fold to approximately 2.5 kB for every 100 m of a profile (when profiling at 0.2 m s⁻¹). As a test, we apply our scheme to a dataset comprising 650 profiles and compare its output to that from our standard turbulence-processing algorithm. For ϵ, values from the two approaches agree within a factor of 2 87 % of the time; for χ, they agree 78 % of the time. These levels of agreement are greater than or comparable to that between the ϵ and χ values derived from two shear probes and two fast thermistors, respectively, on the same profiler.

1 Introduction

Measurements of oceanic turbulence have been made since the 1950s using platforms and sensors of various shapes and sizes (Lueck et al., 2002). Complete resolution of the turbulence requires measuring temperature and velocity gradients at millimeter-to-centimeter scale. Hence, sampling turbulence is data intensive. Whereas conventional profiling measurements of temperature, conductivity, and pressure are typically sampled at 1 Hz (e.g., Argo floats; Roemmich et al., 2019a), a turbulence profile involves sampling multiple sensors at 100–1000 Hz. A relatively minimal requirement of five separate signals sampled at 100 Hz and recorded at 16-bit resolution equates to 1 kB s⁻¹ or 500 kB per 100 m of profiling range at 0.2 m s⁻¹ profiling speed. For floats, this is not a trivial volume of data. For example, transmitting only 3 kB of data from a Deep SOLO float takes 100–200 s (Roemmich et al., 2019b). Extended surfacings are also presented with the danger of surface vessels and vandals. Ultimately, raw turbulence profiles are 2–3 orders of magnitude too large to transmit in a reasonable amount of time.

One approach to reducing turbulence data is given by Rainville et al. (2017), who use it for multi-month glider missions. On board the glider, they calculate spectra of raw voltage signals reported by the shear probes and fast thermistors and then band average each of these spectra into 12 bins. After transmission, these binned values are calibrated and fit to model spectra. In other words, they (i) postpone calibration and (ii) minimize the data manipulation and processing that happens on board. These two strategies are shared by our scheme (and also shared by the reduction scheme developed for χpods; Becherer and Moum, 2017). Otherwise, however, our scheme differs from that of Rainville et al. (2017), as it does not suit our scientific goals of measuring turbulence over the upper ∼ 120 m at high vertical resolution (e.g., ∼ 1 m) and as frequently as possible (see Sect. 2). Two shear probes and two fast thermistors would produce $(\mathrm{2}+\mathrm{2})\times \mathrm{12}=\mathrm{48}$ spectral values per segment. Even without considering the other profiling quantities, the spectral values could add up to >20 kB per dive for our scenario. We would be spending as much time transmitting the data as actually measuring the ocean.

Our scheme is for profiling instruments that contain shear probes and, optionally, fast thermistors (Sect. 2). First, we document the necessary calibration details (Sect. 3). Next, we compress raw shear voltages by way of simple power law fits and show how ϵ is derived from these fits in post-processing (Sect. 4). A test of the scheme employing 650 profiles demonstrates that little accuracy is sacrificed in return for a large reduction in data volume (Sect. 5). Knowing ϵ from the shear probe measurements makes possible a similar method for deriving χ from fast-thermistor measurements (Sects. 6 and 7). Adapting the scheme to a different profiler requires minimal modification (Sect. 8). For our setup, the scheme reduces the dataset size by a factor of ∼ 300: only 2.5 kB for each 100 m of a profile (Sect. 9).

2 The Flippin' χSOLO (FCS)

We intend our data reduction scheme to be sufficiently general so as to be portable to all vertical turbulence profilers that contain shear probes. It can also be used with gliders if a measure of flow speed past the sensors is available (e.g., Greenan et al., 2001; Merckelbach and Carpenter, 2021). In a general sense, some of the values specified herein ought to be considered variables (Sect. 8). However, we do have a particular platform for which we are developing the scheme, the Flippin' χSOLO (FCS), and the values used here are chosen for the objective of detailed upper-ocean profiling.

FCS is a descendant of the SOLO-II profiling float (Roemmich et al., 2004) with the addition of a turbulence package and extra functionality. The turbulence package includes two shear probes (Osborn, 1974) to measure small-scale velocity gradients from which ϵ is computed, two fast thermistors to measure small-scale temperature fluctuations from which χ is computed, and a pressure sensor from which profiling speed is derived. FCS also includes a three-axis accelerometer that is used to measure the surface wave field (although with a method not described in this paper). Accelerometer data from when FCS is profiling are not used in our reduction scheme. FCS and its measurements are described more completely in the companion paper (Moum et al., 2023).

To reverse profiling direction, FCS adjusts buoyancy and flips via internal shifting of the battery pack. This causes the turbulence sensors to always point into undisturbed fluid. Flipping therefore permits profiling on both descent and ascent, including sampling of the upper 5 m on the ascent.

As a prototype, a standard (non-flipping) SOLO float with a modified χpod (Moum and Nash, 2009) attached was deployed in the Bay of Bengal to measure the suppression of turbulence by salinity stratification (Shroyer et al., 2016). This instrument – named χSOLO – did not have shear probes and therefore could not have provided estimates of ϵ within mixed layers. (Values of ϵ can be approximated from χ, but only if there is stratification.) Nevertheless, χSOLO's success motivated the development of the FCS units with flipping capabilities and fully integrated turbulence packages. These new instruments retained the χ or C in their name despite their ability to also directly measure ϵ from shear probes.

Two FCS units were vetted over 4 days in May 2019 off the Oregon coast. During this period, each unit profiled from the surface to ∼ 120 m and back at a typical speed of 0.2 m s⁻¹. Adding time at the surface, each dive cycle took ∼ 30 min, and we obtained 650 profiles in total. In this 2019 experiment, one of the shear probes on one of the two units malfunctioned. Hence, the dataset for this paper contains approximately 25 % less shear data than fast-thermistor data.

3 Conversion of measured voltages to physical units

The core of our data reduction scheme uses power law fits of voltage spectra that are calculated on board and subsequently converted to meaningful turbulence quantities in post-processing. Additional voltage quantities are also recorded to determine temperature, pressure, and profiling speed.

3.1 Nomenclature and conventions

• All quantities measured by FCS that are discussed in this paper are sampled at 100 Hz.

• All voltage spectra are frequency spectra and are denoted by Ψ_x(f) (where x is a label, such as s for shear) with units of V² Hz⁻¹.

• Physical spectra of shear and temperature gradient are wavenumber spectra and are denoted by Φ_x(k) with units of s⁻² cpm⁻¹ and K² m⁻² cpm⁻¹, respectively. Figure 4 is an exception in which shear spectra are frequency spectra: Φ_s(f).

• Wavenumber k has the unit cycles per meter (cpm). Expressions quoted from other papers may differ by factors of 2π for wavenumbers in radians per meter.

• The Kraichnan model spectrum Φ_Kr primarily depends on the dissipation rates of turbulent kinetic energy and temperature variance (ϵ and χ), but it also depends on the molecular viscosity ν and molecular thermal diffusivity D_T. For brevity, we write ${\mathrm{\Phi }}_{\mathrm{Kr}}(k,\mathit{ϵ},\mathit{\chi })$ rather than the more complete ${\mathrm{\Phi }}_{\mathrm{Kr}}(k,\mathit{ϵ},\mathit{\chi },\mathit{\nu },D_{\mathrm{T}})$. Similarly, the Nasmyth spectrum is written as Φ_Na(k,ϵ) rather than ${\mathrm{\Phi }}_{\mathrm{Na}}(k,\mathit{ϵ},\mathit{\nu })$. In cases where the arguments are unambiguous or unimportant, we simply write Φ_Na and Φ_Kr.

• A pair of angle brackets, 〈⋅〉, denotes the mean value over a segment of length: N_seg=512 data points. This equates to ∼ 1 m at our nominal profiling speed of 0.2 m s⁻¹.

• To calculate spectra for a given 512-element voltage segment, we first remove the linear trend, then we use three half-overlapping, Hamming-windowed, 256-element subsegments (i.e., N_fft=256; N_overlap=128).

In general, the values of N_seg and N_fft are variables. Our choices are based on the 100 Hz sampling (∼ 500 cpm) and the goals of FCS, which include obtaining high-vertical-resolution turbulence data, especially near the surface. For different turbulence profilers or different scientific goals, longer segments and/or more overlapping subsegments may be more appropriate (see Sect. 8).

We do not pursue the possibility of using accelerometers to decontaminate spectra (e.g., Levine and Lueck, 1999); three subsegments is too few for this to work well. Rather, we focus on a frequency band that we know to be unaffected by vibration.

3.2 Shear calibration

The voltage reported by the shear probe V_s is linearly proportional to shear:

$$\begin{array}{}\text{(1)}& {\displaystyle }{u}_z={\displaystyle \frac{\mathit{\alpha }}{W^{\mathrm{2}}}}{V}_{\mathrm{s}},\text{(2)}& {\displaystyle }\mathit{\alpha }=\mathrm{1}/\left(\mathrm{2}\sqrt{\mathrm{2}}\mathit{\rho }{G}_{\mathrm{s}}{T}_{\mathrm{s}}{S}_{\mathrm{s}}\right),\end{array}$$

where W is the flow speed past the sensor. The overall engineering calibration α includes the seawater density ρ, the analog circuit gain G_s (equal to 1 for FCS circuitry), the probe sensitivity S_s (∼ $\mathrm{0.25}\times {\mathrm{10}}^{-\mathrm{3}}$ V m² N⁻¹), and the differentiator time constant T_s (∼ 1 s). The linearity in V_s admits a simple link between the physical and voltage spectra:

$$\begin{array}{}\text{(3)}& {\mathrm{\Phi }}_{u_z}\left(k\right)={\displaystyle \frac{\mathrm{1}}{H_{\mathrm{s}}^{\mathrm{2}}\left(k\right)}}{\displaystyle \frac{{\mathit{\alpha }}^{\mathrm{2}}}{W^{\mathrm{3}}}}{\mathrm{\Psi }}_{\mathrm{s}}\left(f\right),\end{array}$$

where $H_{\mathrm{s}}^{\mathrm{2}}\left(k\right)$ is the transfer function that accounts for (i) spatial averaging by the shear probe of high-wavenumber motions and (ii) analog and digital filtering of the raw voltage signal (see Appendix A). Note also the use above of the following relation:

$$\begin{array}{}\text{(4)}& \mathrm{\Psi }\left(f\right)=\mathrm{\Psi }\left(k\right){\displaystyle \frac{\mathrm{d}k}{\mathrm{d}f}}={\displaystyle \frac{\mathrm{\Psi }\left(k\right)}{W}}.\end{array}$$

3.3 Temperature and temperature gradient calibration

Two voltage signals are recorded for each fast thermistor. V_T is the voltage output directly related to T, and V_Tt is the differentiated output, which improves resolution at high frequencies (≳ 10 Hz). Temperature is related to V_T through a quadratic calibration:

$$\begin{array}{}\text{(5)}& {\displaystyle }T=C_{\mathrm{1}T}+C_{\mathrm{2}T}{V}_T+C_{\mathrm{3}T}{V}_T^{\mathrm{2}},\text{(6)}& {\displaystyle }\langle T\rangle =C_{\mathrm{1}T}+C_{\mathrm{2}T}\langle V_T\rangle +C_{\mathrm{3}T}\langle V_T{\rangle }^{\mathrm{2}},\end{array}$$

where C_1T, C_2T, and C_3T are coefficients determined from laboratory calibrations. Equation (6) is technically an approximation because it contains 〈V_T〉² and not $\langle V_T^{\mathrm{2}}\rangle$, but over 5 s time scales this changes 〈T〉 by ≲ 0.001 ^∘C (estimated from the 2019 dataset).

The gradient of the temperature calibration is

$$\begin{array}{}\text{(7)}& {\displaystyle \frac{\partial T}{\partial V_T}}=C_{\mathrm{2}T}+\mathrm{2}{C}_{\mathrm{3}T}{V}_T=C_{\mathrm{2}T}+\mathrm{2}{C}_{\mathrm{3}T}\langle V_T\rangle .\end{array}$$

Consequently, the small-scale vertical temperature gradient T_z is linearly proportional to the differentiated voltage V_Tt. To demonstrate, we first rewrite T_z in terms of more directly measured quantities:

$$\begin{array}{}\text{(8)}& T_z={\displaystyle \frac{\partial T}{\partial z}}={\displaystyle \frac{\partial T}{\partial V_T}}{\displaystyle \frac{\partial V_T}{\partial t}}{\displaystyle \frac{\partial t}{\partial z}}.\end{array}$$

The first quantity on the right-hand side is Eq. (7), the last is $\mathrm{1}/W$, and the second is

$$\begin{array}{}\text{(9)}& {\displaystyle \frac{\partial V_T}{\partial t}}={\displaystyle \frac{V_{Tt}}{C_{Tt}}},\end{array}$$

where C_Tt is the gain of the analog differentiator.

Rewriting Eq. (8), the aforementioned linear relationship between T_z and V_Tt becomes

$$\begin{array}{}\text{(10)}& T_z=\left({\displaystyle \frac{C_{\mathrm{2}T}+\mathrm{2}{C}_{\mathrm{3}T}\langle V_T\rangle }{C_{Tt}W}}\right)V_{Tt}.\end{array}$$

The relationship between physical and voltage spectra is therefore

$$\begin{array}{}\text{(11)}& {\mathrm{\Phi }}_{T_z}\left(k\right)={\displaystyle \frac{\mathrm{1}}{H_{Tt}^{\mathrm{2}}\left(k\right)}}{\left({\displaystyle \frac{C_{\mathrm{2}T}+\mathrm{2}{C}_{\mathrm{3}T}\langle V_T\rangle }{C_{Tt}}}\right)}^{\mathrm{2}}{\displaystyle \frac{\mathrm{1}}{W}}{\mathrm{\Psi }}_{Tt}\left(f\right).\end{array}$$

Again, we have invoked Eq. (4), and the transfer function $H_{Tt}^{\mathrm{2}}\left(k\right)$ is defined in Appendix A.

3.4 Pressure and profiling-velocity calibration

Pressure has a linear calibration:

$$\begin{array}{}\text{(12)}& P\left[\mathrm{dbar}\right]={\displaystyle \frac{C_{\mathrm{1}P}+C_{\mathrm{2}P}{V}_P}{\mathrm{1.45}\mathrm{psi}{\mathrm{dbar}}^{-\mathrm{1}}}}-p_{\mathrm{atm}}.\end{array}$$

In our usage, the coefficients C_1P and C_2P are recorded in units of psi and psi V⁻¹, respectively, and calibrated under total pressure. Subtracting atmospheric pressure makes P=0 at the sea surface. The constant C_1P must account for the vertical position of the pressure sensor on the instrument relative to the shear probes and thermistors. Hence, C_1P differs between upcasts and downcasts. For the reduced dataset, we record the last pressure voltage in each segment. For example, with N_seg=512, we save the 512th, 1024th, and 1536th values of V_P for the first three segments. The average pressures in the second and third segments are 0.5(V_P(1024)+V_P(512)) and 0.5(V_P(1536)+V_P(1024)), respectively. The average pressure in the first segment is found by extrapolation.

The flow speed past the sensors, denoted W, is derived from the rate of change of the pressure voltages just described:

$$\begin{array}{}\text{(13)}& \langle W\rangle ={\displaystyle \frac{C_{\mathrm{2}P}}{\mathrm{1.45}}}{\displaystyle \frac{\left|\mathrm{\Delta }{V}_P\right|}{N_{\mathrm{seg}}\mathrm{\Delta }t}},\end{array}$$

where Δt is the sampling period (here 0.01 s), and ΔP is V_P(1024)−V_P(512) and V_P(1536)−V_P(1024) for the second and third segments, respectively. Extrapolation is again used for the first segment.

No smoothing is necessary before calculating ΔV_P because its magnitude is so much larger than the quantization of the signal (this being the limiting factor for precision of pressure recorded by FCS). In physical units, P is precise to 0.003 dbar, which is 𝒪(300) times smaller than ΔP.

Wave orbitals can introduce variability when W is small (≲ 0.15 m s⁻¹). As a diagnostic we calculate and record the minimum value of W for each segment. This also helps to identify the beginning and end of profiles, as shown in Appendix B. In standard processing, we would derive W(t) from the pressure signal low passed at 2 Hz. To avoid the need to low-pass filter the signal on board, we instead make 10 estimates of W(t) per segment and take the minimum of these:

$$\begin{array}{}\text{(14)}& W_{\min }={\displaystyle \frac{C_{\mathrm{2}P}}{\mathrm{1.45}}}\underset{\text{Calculate on board}}{\underbrace{min\left({\displaystyle \frac{\left|\mathrm{\Delta }{V}_P\left(t_i\right)\right|}{\mathrm{50}\mathrm{\Delta }t}}\right)}},\end{array}$$

where t_i=1, 51, 101, …, 501. Even with this sampling of every 50th element, which follows from subsampling a 100 Hz signal at 2 Hz, ΔV_P(t_i) is large enough that smoothing is unnecessary.

In this paper, we immediately discard all segments in which W_min<0.05 m s⁻¹. This threshold is reached only at the top and bottom of profiles, if at all. Note, however, that this does not imply that a segment with W_min>0.05 m s⁻¹ is trustworthy. Even segments with W_min closer to 0.15 m s⁻¹ should be treated with particular caution. Signs that a segment is questionable are that 〈W〉 and W_min differ by ∼ 20 % or more and that spectral fit scores are low (see Sects. 4.3 and 6.3). These two issues often co-occur because of the nonlinear relationship between shear and profiling speed (Eq. 1). Low-frequency variations in W ultimately lead to spectra that are redder than expected and hence have low fit scores. Such segments should be discarded. There is not a simple way to correct the spectra given the nonlinearity.

4 Reduction of shear data

In this section, we are ultimately going to fit measured spectra to an inertial subrange model that does not necessarily apply at the relevant frequencies or wavenumbers. We will elaborate as we go, but we want to emphasize in advance that measured spectra do not need to conform to an inertial subrange model for us to obtain accurate values of ϵ. The inertial subrange is merely a convenient starting point.

4.1 Summarizing Nasmyth spectra with $f^{\mathrm{1}/\mathrm{3}}$ fits

Shear measurements ideally capture both the inertial and viscous subranges and hence use a wide band of the measured spectrum to derive values for ϵ. In practice, noise and sensor resolution limit how well the true environmental spectrum is resolved. Conventional work-arounds exploit the Nasmyth model spectrum Φ_Na(k,ϵ) (Nasmyth, 1970; Oakey, 1982). One approach is to iterate toward a solution in which the integral of Φ_Na over a specific wavenumber band matches that of the measured spectrum ${\mathrm{\Phi }}_{u_z}$ (e.g., Moum et al., 1995). Another is to find the best fit of ${\mathrm{\Phi }}_{u_z}$ to Φ_Na by using maximum-likelihood estimation together with a model of the expected statistical distribution of the spectral coefficients being fitted (e.g., Bluteau et al., 2016).

Here we develop a new and simpler two-stage approach to fitting shear spectra to Φ_Na. In the first stage, we use an $f^{\mathrm{1}/\mathrm{3}}$ power law fit over a fixed-frequency range of f_l to f_h = 1–5 Hz, where $f^{\mathrm{1}/\mathrm{3}}$ follows from the assumption that we are fitting over the inertial subrange. In the second stage, we correct for this often-invalid assumption.

In the inertial subrange, shear spectra are proportional to $k^{\mathrm{1}/\mathrm{3}}$ and hence also to $f^{\mathrm{1}/\mathrm{3}}$, since f=Wk. With N_fft=256 and 100 Hz sampling (Sect. 3.1), spectral coefficients are separated by frequency increments of 100 Hz $/$ 256 = 0.39 Hz, so there are 10 coefficients between 1 and 5 Hz. (Our processing code will actually use bounding frequencies of 0.98 and 4.88 Hz, as these are half-integer multiples of 0.39 Hz, but for brevity we will write these as 1 and 5 Hz throughout.)

Our choice of f_l=1 Hz is dictated by a requirement that we avoid low-frequency contamination induced by (i) advection by wave orbital motion and (ii) pitch and roll motions of the profiler. Together, these dominate below 0.3 Hz. Setting f_l=0.5 Hz would add only one more spectral coefficient. Our choice of f_h=5 Hz is a trade-off between maximizing the bandwidth of the fit and minimizing how much the measured spectra are subject to either noise or viscous roll off. Other profilers may benefit from different frequency bounds (see Sect. 8).

Our inertial subrange assumption is often false. Indeed, “assumption” is perhaps a misnomer, as we do not expect it to be true; we know that viscous roll off will often occur at frequencies lower than 5 Hz (25 cpm for a nominal value of W=0.2 m s⁻¹). However, because there exists an analytical expression for the viscous roll off, we are able to derive an exact expression that quantifies how much ϵ is underestimated. This is the second stage of our approach. We derive an expression for the correction function F_Na in such a way that it can be calculated in post-processing. The benefits of this approach are that (i) we can fit uncalibrated (i.e., voltage) spectra, and (ii) it simplifies the actual onboard fitting routine (Sect. 4.2).

The full Nasmyth spectrum and its inertial range approximation are as follows (Lueck, 2013):

$$\begin{array}{}\text{(15)}& {\displaystyle }{\mathrm{\Phi }}_{\mathrm{Na}}(k,\mathit{ϵ})={\displaystyle \frac{{\mathit{ϵ}}^{\mathrm{3}/\mathrm{4}}}{{\mathit{\nu }}^{\mathrm{1}/\mathrm{4}}}}{\displaystyle \frac{\mathrm{8.05}(k\mathit{\eta }{)}^{\mathrm{1}/\mathrm{3}}}{\mathrm{1}+(\mathrm{20.6}k\mathit{\eta }{)}^{\mathrm{3.715}}}},\text{(16)}& {\displaystyle }{\mathrm{\Phi }}_{\mathrm{Na}}(k\lesssim \mathrm{0.02}/\mathit{\eta },\mathit{ϵ})=\mathrm{8.05}{k}^{\mathrm{1}/\mathrm{3}}{\mathit{ϵ}}^{\mathrm{2}/\mathrm{3}},\end{array}$$

where $\mathit{\eta }=({\mathit{\nu }}^{\mathrm{3}}/\mathit{ϵ}{)}^{\mathrm{1}/\mathrm{4}}$ is the Kolmogorov length scale.

Let ϵ_init denote the initial value of ϵ that comes from fitting a measured spectrum to the approximate form in Eq. (16) using the simple power law fitting method in Appendix C rather than fitting to the full form in Eq. (15). As noted earlier, the fit will be over the f_l–f_h = 1–5 Hz range which, given our nominal value of W=0.2 m s⁻¹, equates to k_l–k_h = 5–25 cpm.

Consider two contrasting examples of low and high turbulence with $\mathit{ϵ}=\mathrm{1}\times {\mathrm{10}}^{-\mathrm{9}}$ and $\mathrm{1}\times {\mathrm{10}}^{-\mathrm{6}}$ W kg⁻¹, respectively (Fig. 1a). For now, assume the measured spectrum to be fit is itself a Nasmyth spectrum. For $\mathit{ϵ}={\mathrm{10}}^{-\mathrm{6}}$ W kg⁻¹, the $f^{\mathrm{1}/\mathrm{3}}$ fit lies on top of Φ_Na. Conversely, the $f^{\mathrm{1}/\mathrm{3}}$ fit for the smaller ϵ value is seemingly meaningless: the $f^{\mathrm{1}/\mathrm{3}}$ fit (dashed line) does not even match the sign of the slope of Φ_Na. Worse yet, naively inverting this initial fit produces the underestimate ${\mathit{ϵ}}_{\mathrm{init}}=\mathrm{1.2}\times {\mathrm{10}}^{-\mathrm{10}}$ W kg⁻¹, which is 8 times smaller than the true value of ϵ. However, by adjusting by a factor of $\mathrm{1}/F_{\mathrm{Na}}$, defined in the following paragraph, the fit (dotted line) looks like a hypothetical extrapolation of the inertial subrange. Equivalently, ϵ_init is corrected to the true value of ϵ as

$$\begin{array}{}\text{(17)}& \mathit{ϵ}={\mathit{ϵ}}_{\mathrm{init}}/F_{\mathrm{Na}}^{\mathrm{3}/\mathrm{2}}.\end{array}$$

In our example, $\mathrm{1}\times {\mathrm{10}}^{-\mathrm{9}}$ W kg⁻¹ = $\mathrm{1.2}\times {\mathrm{10}}^{-\mathrm{10}}$ W kg⁻¹ $/$ ${\mathrm{0.238}}^{\mathrm{3}/\mathrm{2}}$. The value of 0.238 is the solution to an implicit equation, derived below, that depends on ϵ_init and W. For clarity, our demonstration starts by assuming that we know ϵ rather than ϵ_init.

Figure 1 Calculation of the correction function F_Na for two values of ϵ. For $\mathit{ϵ}=\mathrm{1}\times {\mathrm{10}}^{-\mathrm{6}}$ W kg⁻¹, a $k^{+\mathrm{1}/\mathrm{3}}$ power law is a good approximation of the Nasmyth spectrum over the frequency range f_l–f_h (1–5 Hz) for a profiling speed of W=0.2 m s⁻¹. Although the same is not true for $\mathit{ϵ}=\mathrm{1}\times {\mathrm{10}}^{-\mathrm{9}}$ W kg⁻¹, we can account for the roll off with a factor of F_Na. F_Na can be defined in terms of either ϵ (Eq. 18) or ϵ_init (Eq. 19). Panel (b) takes the former approach. In practice, we must take the latter approach, since we do not know ϵ until after it is derived from ϵ_init and F_Na.

Nasmyth spectra can be flattened to unity over the inertial subrange with the normalization $\mathrm{8.05}{k}^{\mathrm{1}/\mathrm{3}}{\mathit{ϵ}}^{\mathrm{2}/\mathrm{3}}$ (Fig. 1b). Values of F_Na are based on the mean of these flattened spectra over the wavenumber range k_l–k_h ($=f_{\mathrm{l}}/W$–$f_{\mathrm{h}}/W$):

$$\begin{array}{}\text{(18)}& F_{\mathrm{Na}}={\displaystyle \frac{\mathrm{1}}{k_{\mathrm{h}}-k_{\mathrm{l}}}}\underset{k_{\mathrm{l}}}{\overset{k_{\mathrm{h}}}{\int }}{\displaystyle \frac{{\mathrm{\Phi }}_{\mathrm{Na}}(k,\mathit{ϵ})}{\mathrm{8.05}{k}^{\mathrm{1}/\mathrm{3}}{\mathit{ϵ}}^{\mathrm{2}/\mathrm{3}}}}\mathrm{d}k.\end{array}$$

To remove the dependence of the true value of ϵ, we substitute Eq. (17) to produce an implicit function for F_Na, which can be solved numerically:

$$\begin{array}{}\text{(19)}& {\displaystyle \frac{\mathrm{1}}{k_{\mathrm{h}}-k_{\mathrm{l}}}}\underset{k_{\mathrm{l}}}{\overset{k_{\mathrm{h}}}{\int }}{\displaystyle \frac{{\mathrm{\Phi }}_{\mathrm{Na}}\left(k,{\mathit{ϵ}}_{\mathrm{init}}/F_{\mathrm{Na}}^{\mathrm{3}/\mathrm{2}}\right)}{\mathrm{8.05}{k}^{\mathrm{1}/\mathrm{3}}{\mathit{ϵ}}_{\mathrm{init}}^{\mathrm{2}/\mathrm{3}}/F_{\mathrm{Na}}}}\mathrm{d}k-F_{\mathrm{Na}}=\mathrm{0}.\end{array}$$

Note how the two forms of F_Na (Eqs. 18 and 19) are defined with different arguments. For our example, F_Na ($\mathit{ϵ}=\mathrm{1}\times {\mathrm{10}}^{-\mathrm{9}}$ W kg⁻¹) = $F_{\mathrm{Na}}({\mathit{ϵ}}_{\mathrm{init}}=\mathrm{1.2}\times {\mathrm{10}}^{-\mathrm{10}}$ W kg⁻¹) = 0.238. Hereafter, we use the latter: F_Na(ϵ_init).

With f_l and f_h fixed, F_Na is a function of three variables: ϵ_init, W, and ν. F_Na is closer to 1 (less of a correction) for larger values of ϵ_init (Fig. 1). It is also closer to 1 for higher values of W (Fig. 2), since k_l and k_h decrease with increasing W (i.e., k_l–k_h move closer to the inertial subrange).

Figure 2 The correction function F_Na for $\mathit{\nu }=\mathrm{1}\times {\mathrm{10}}^{-\mathrm{6}}$ m² s⁻¹ and f_l–f_h = 1–5 Hz.

To simplify calculations in the upcoming section, we make one final change to Eq. (19) using the following substitution:

$$\begin{array}{}\text{(20)}& {\mathrm{\Phi }}_{\mathrm{Na}}(k,\mathit{ϵ})\to H_{\mathrm{s}}^{\mathrm{2}}\left(k\right){\mathrm{\Phi }}_{\mathrm{Na}}(k,\mathit{ϵ}).\end{array}$$

Therefore,

$$\begin{array}{}\text{(21)}& {\displaystyle \frac{\mathrm{1}}{k_{\mathrm{h}}-k_{\mathrm{l}}}}\underset{k_{\mathrm{l}}}{\overset{k_{\mathrm{h}}}{\int }}{\displaystyle \frac{H_{\mathrm{s}}^{\mathrm{2}}\left(k\right){\mathrm{\Phi }}_{\mathrm{Na}}\left(k,{\mathit{ϵ}}_{\mathrm{init}}/F_{\mathrm{Na}}^{\mathrm{3}/\mathrm{2}}\right)}{\mathrm{8.05}{k}^{\mathrm{1}/\mathrm{3}}{\mathit{ϵ}}_{\mathrm{init}}^{\mathrm{2}/\mathrm{3}}/F_{\mathrm{Na}}}}\mathrm{d}k-F_{\mathrm{Na}}=\mathrm{0}.\end{array}$$

Think of this substitution in Eq. (20) as inverting the conventional way that $H_{\mathrm{s}}^{\mathrm{2}}\left(k\right)$ is invoked. Usually, a measured shear spectrum is amplified at high wavenumbers by $\mathrm{1}/H_{\mathrm{s}}^{\mathrm{2}}\left(k\right)$ and is then fit to the model spectrum Φ_Na. Here, instead of amplifying the measured spectrum, we reduce the model spectrum. With this latter approach, $H_{\mathrm{s}}^{\mathrm{2}}\left(k\right)$ is calculated and applied only during the post-processing stage. (It changes F_Na by only ∼ 5 %, since we fit over relatively low wavenumbers.)

4.2 Obtaining ϵ_init from shear voltage spectra

Since ϵ can be reconstructed from ϵ_init, we require an expression linking ϵ_init to the shear voltage spectrum Ψ_s. Equating Eqs. (3) and (16) gives

$$\begin{array}{}\text{(22)}& {\displaystyle \frac{{\mathit{\alpha }}^{\mathrm{2}}}{W^{\mathrm{3}}}}{\mathrm{\Psi }}_{\mathrm{s}}\left(f\right)=\mathrm{8.05}{k}^{\mathrm{1}/\mathrm{3}}{\mathit{ϵ}}_{\mathrm{init}}^{\mathrm{2}/\mathrm{3}},\end{array}$$

where we have left out $H_{\mathrm{s}}^{\mathrm{2}}\left(k\right)$, since it has been incorporated into F_Na. Rearranging and substituting $k=f/W$ gives

$$\begin{array}{}\text{(23)}& {\mathit{ϵ}}_{\mathrm{init}}^{\mathrm{2}/\mathrm{3}}{f}^{\mathrm{1}/\mathrm{3}}={\displaystyle \frac{{\mathit{\alpha }}^{\mathrm{2}}}{\mathrm{8.05}{W}^{\mathrm{8}/\mathrm{3}}}}{\mathrm{\Psi }}_{\mathrm{s}}\left(f\right).\end{array}$$

Then, to solve for ϵ_init, we use a least-squares fit (see Appendix C):

$$\begin{array}{}\text{(24)}& {\mathit{ϵ}}_{\mathrm{init}}^{\mathrm{2}/\mathrm{3}}={\displaystyle \frac{{\mathit{\alpha }}^{\mathrm{2}}}{\mathrm{8.05}{W}^{\mathrm{8}/\mathrm{3}}}}\underset{\text{Calculate on board}}{\underbrace{{\displaystyle \frac{{\sum }_f{f}^{\mathrm{1}/\mathrm{3}}{\mathrm{\Psi }}_{\mathrm{s}}}{{\sum }_f{f}^{\mathrm{2}/\mathrm{3}}}}}},\end{array}$$

where the sums are understood to be over the range f_l–f_h. The quantities α, W, and ϵ_init are calculated in post-processing.

4.3 Quality control of the shear spectral fits

Measured shear spectra are often quality controlled either by manual visual inspection or, more objectively, by quantifying the level of mismatch between them and their associated model. Possible mismatch quantities include the mean absolute deviation or the variance of the ratio ${\mathrm{\Phi }}_{u_z}/{\mathrm{\Phi }}_{\mathrm{Na}}$ (e.g., Ruddick et al., 2000; Bluteau et al., 2016). We cannot calculate such quantities with our reduced scheme because we do not know what each spectrum should look like until we calculate its ϵ value in the post-processing stage. (Recall that Φ_Na is a function of ϵ.) By this stage, we have lost information about the spectral shape through the summing operation in Eq. (24).

To retain at least some information about the shape of each voltage spectrum, we will split the 1–5 Hz range and compute two fits rather than one. Doing so allows for a first-order check that the spectrum over the 1–5 Hz range approximately follows the expected shape.

Figure 3 A visual demonstration of how the ϵ fit score (Eq. 26) characterizes better and worse fits. For all three of these hypothetical spectra, ϵ values from the two fits (1–3 and 3–5 Hz) average to $\mathit{ϵ}=\mathrm{1.0}\times {\mathrm{10}}^{-\mathrm{8}}$ W kg⁻¹. Only in panel (a), however, does the measured spectrum agree well with the Nasmyth spectrum for this ϵ value. In practice, the initial fits would be undertaken on voltage spectra; here, we are using physical units for simplicity. In all three examples, the kinematic viscosity is $\mathit{\nu }={\mathrm{10}}^{-\mathrm{6}}$ m² s⁻¹.

Figure 4 Examples of measured shear spectra exhibiting a range of ϵ fit scores (Eq. 26). The best fits are at the top with progressively worse fits (lower scores) moving downward. The examples in each column have (left) $\mathit{ϵ}={\mathrm{10}}^{-\mathrm{9}}$–10⁻⁸, (middle) $\mathit{ϵ}={\mathrm{10}}^{-\mathrm{8}}$–10⁻⁷, and (right) $\mathit{ϵ}={\mathrm{10}}^{-\mathrm{7}}$–10⁻⁶ W kg⁻¹. Each score is only based on spectral coefficients from 1–5 Hz, but lower and higher frequencies are shown for reference.

Mathematically, there is nothing special about our choice f_l–f_h=1–5 Hz. In theory, we can split the 1–5 Hz range into two (1–3 and 3–5 Hz) and obtain a value of ϵ_init for each. These values will differ, but so will the associated values of F_Na. For a measured spectrum that conforms to a Nasmyth spectrum, the two values of ϵ calculated with Eq. (17) will not differ (Fig. 3). We therefore calculate on board the sums in Eq. 24 over both f_l–f_m and f_m–f_h, where the middle frequency f_m=3 Hz. (In our code, f_m is actually 7.5×0.39 Hz = 2.93 Hz for the reason given in Sect. 4.1.) Hence, for each spectrum, we are able to post-process to recover two independent estimates of ϵ, denoted ϵ_l–m and ϵ_m–h. The mean of these two provides a single, final value for ϵ, and their ratio quantifies the match of a measured spectrum to a Nasmyth spectrum over the range f_l–f_h:

$$\begin{array}{}\text{(25)}& {\displaystyle }\mathit{ϵ}=\mathrm{mean}\left({\mathit{ϵ}}_{\text{l–m}},{\mathit{ϵ}}_{\text{m–h}}\right),\text{(26)}& {\displaystyle }\mathit{ϵ}\text{ fit score}={\displaystyle \frac{min\left({\mathit{ϵ}}_{\text{l–m}},{\mathit{ϵ}}_{\text{m–h}}\right)}{max\left({\mathit{ϵ}}_{\text{l–m}},{\mathit{ϵ}}_{\text{m–h}}\right)}}.\end{array}$$

The best possible fit score is 1; the lower the score, the poorer the fit. The example spectra in Fig. 4 show that a high fit score does not necessarily imply small residuals. Rather, fits with high scores are typically those with random residuals, meaning that a given measured spectral coefficient is just as likely to be above the fit as it is to be below it. Fits with low fit scores are typically those with autocorrelated residuals, meaning that the sign and/or magnitude of a residual is correlated with that of its neighbors. In practice, we expect a range of ϵ fit scores: instantaneous and unaveraged spectra differ from the Nasmyth spectrum because they are derived from a limited sampling of a statistical process and because of non-stationarity, anisotropy, and inhomogeneity of the turbulence.

When ϵ is small (≲ 10⁻⁹ W kg⁻¹), the fit score may be consistently low if spectral coefficients in the f_m–f_h range are affected by noise, and consequently ϵ_m–h≫ϵ_l–m. For such cases, we choose to use only the lower-frequency fit. We would rather have a more accurate estimate of ϵ and forgo the fit score than have a biased-high ϵ value with a biased-low fit score. (Either way, the small values of ϵ in question will have minimal effect on any averages, given that turbulence distributions have high kurtosis, so high values dominate means.) Specifically,

$$\begin{array}{}\text{(27)}& \left.\begin{array}{ll}\mathit{ϵ}& ={\mathit{ϵ}}_{\text{l–m}}\\ \mathit{ϵ}\text{ fit score}& =\text{–}\end{array}\right\}\text{if}\mathrm{0.1}W{\left({\displaystyle \frac{{\mathit{ϵ}}_{\text{l–m}}}{{\mathit{\nu }}^{\mathrm{3}}}}\right)}^{\mathrm{1}/\mathrm{4}}<f_{\mathrm{m}},\end{array}$$

where the threshold is equivalent to $k<\mathrm{0.1}/\mathit{\eta }$, with η being the Kolmogorov length scale estimated from ϵ_l–m. For reference, Φ_Na peaks at $k=\mathrm{0.026}/\mathit{\eta }$ and rolls off to 11 % of its maximum by $k=\mathrm{0.1}/\mathit{\eta }$ (see Eq. 15).

5 Test of the reduction scheme for ϵ

To test the accuracy of the shear reduction scheme described in the previous section, we apply it retrospectively to the dataset from the 2019 test cruise (Sect. 2). We compare the results to those obtained with the standard processing scheme. This standard scheme (Appendix D) features a more sophisticated despiking routine than that used for our reduced scheme, which employs a 3-standard-deviations threshold filter (Appendix E).

A profile-by-profile comparison of the two schemes is shown in Fig. 5. The comparison is then extended to all 650 profiles (> 77 000 segments of shear), where we find that ϵ from the reduced scheme (${\mathit{ϵ}}_{\mathrm{init}}/F_{\mathrm{Na}}^{\mathrm{3}/\mathrm{2}}$) is within a factor of 2 of that from the standard scheme 87 % of the time over the full range of measured values, ${\mathrm{10}}^{-\mathrm{10}}<\mathit{ϵ}<{\mathrm{10}}^{-\mathrm{4}}$ W kg⁻¹ (Fig. 6a–b). For comparison, in only 72 % do we obtain a factor-of-2 agreement between the two independent values of ϵ measured on the unit with two working shear probes (not shown). Further, to obtain this 87 % agreement, we clearly need the correction function F_Na: Fig. 6c shows that the uncorrected values ϵ_init only have 1 : 1 agreement with ϵ from the standard scheme if $\mathit{ϵ}\gtrsim {\mathrm{10}}^{-\mathrm{7}}$ W kg⁻¹. For the lowest values of ϵ, the ratio is closer to 1 : 30.

Figure 5 Testing the proposed data reduction scheme for shear measurements against the standard processing approach. One upward and one downward profile from each of the two FCS units were arbitrarily chosen for this comparison.

Figure 6 Statistical test of the proposed data reduction scheme for ϵ based on all 650 profiles (77 000 segments). (a) A comparison that includes the dependence on ϵ. (b) Further summarized data that exclude this dependence. (c) Same as for panel (a) but uncorrected (F_Na = 1).

To demonstrate the ability of the ϵ fit score to characterize spectra, we show two-dimensional histograms of non-dimensionalized spectral coefficients from all 77 000 measured shear spectra separated into three classes based on their ϵ fit score: 0.67–1.00, 0.33–0.67, and 0.00–0.33. Only the lowest-scoring class fails to collapse to the Nasmyth spectrum (Fig. 7c, f).

Figure 7 As the ϵ fit score decreases from top to bottom, there is a corresponding decrease in the level of agreement and tightness of spread between (i) non-dimensionalized, measured shear spectra and (ii) the Nasmyth spectrum. The two-dimensional histograms in the left column include only spectral coefficients with frequencies between f_l and f_h; those in the right column include all frequencies. The total number of spectra in this figure is lower than in Fig. 10 because some low-ϵ spectra do not have scores (Eq. 27) and because one of the shear probes failed (Sect. 2).

6 Reduction of fast-thermistor data

The scheme to reduce fast-thermistor data to enable measurement of χ is much like the scheme to reduce shear data. As in Sect. 4, we first show how we summarize a model spectrum in terms of a power law fit and a correction factor. (In this case, the correction factor partly depends on the values of ϵ calculated in Sect. 4.) Then we derive the implementation in terms of voltages and calculate a spectral-fit metric.

6.1 Summarizing Kraichnan spectra with f¹ fits

Here we take the Kraichnan spectrum Φ_Kr (Kraichnan, 1968) as our model; for its low-wavenumber approximation, we use the viscous–convective subrange, which scales as k⁺¹. In units of K² m⁻² cpm⁻¹, Φ_Kr and its approximation are as follows (e.g., Peterson and Fer, 2014):

$$\begin{array}{}\text{(28)}& {\displaystyle }{\mathrm{\Phi }}_{\mathrm{Kr}}(k,\mathit{ϵ},\mathit{\chi })=\mathrm{4}{\mathit{\pi }}^{\mathrm{2}}k\mathit{\chi }q\sqrt{\mathit{\nu }/\mathit{ϵ}}\exp \left(-\sqrt{\mathrm{6}q}\mathrm{2}\mathit{\pi }k{\mathit{\lambda }}_B\right),\text{(29)}& {\displaystyle }{\mathrm{\Phi }}_{\mathrm{Kr}}(k\ll {\mathit{\lambda }}_B^{-\mathrm{1}},\mathit{ϵ},\mathit{\chi })=\mathrm{4}{\mathit{\pi }}^{\mathrm{2}}k\mathit{\chi }q\sqrt{\mathit{\nu }/\mathit{ϵ}},\end{array}$$

where the Batchelor length scale ${\mathit{\lambda }}_B=(\mathit{\nu }{D}_{\mathrm{T}}^{\mathrm{2}}/\mathit{ϵ}{)}^{\mathrm{1}/\mathrm{4}}$, and q is a constant taken to be 5.26. This expression does not include a $k^{+\mathrm{1}/\mathrm{3}}$ inertial–convective subrange, which we ignore here, as it increases by less than 1 % the integral of the temperature gradient spectrum from k=0 to k=∞ and therefore has negligible effect on our results.

A fit against Eq. (29) can be rearranged to give χ_init, which is related to χ through the correction function F_Kr as

$$\begin{array}{}\text{(30)}& \mathit{\chi }={\mathit{\chi }}_{\mathrm{init}}/F_{\mathrm{Kr}}.\end{array}$$

F_Kr is not raised to a power like F_Na (Eq. 17). For small values of k, Φ_Kr∝χ, whereas ${\mathrm{\Phi }}_{\mathrm{Na}}\propto {\mathit{ϵ}}^{\mathrm{2}/\mathrm{3}}$.

The derivation of F_Kr is equivalent to F_Na. We therefore present only the result:

$$\begin{array}{}\text{(31)}& {\displaystyle \frac{\mathrm{1}}{k_{\mathrm{h}}-k_{\mathrm{l}}}}\underset{k_{\mathrm{l}}}{\overset{k_{\mathrm{h}}}{\int }}{\displaystyle \frac{H_{Tt}^{\mathrm{2}}\left(k\right){\mathrm{\Phi }}_{\mathrm{Kr}}(k,\mathit{ϵ},{\mathit{\chi }}_{\mathrm{init}}/F_{\mathrm{Kr}})}{\mathrm{4}{\mathit{\pi }}^{\mathrm{2}}k({\mathit{\chi }}_{\mathrm{init}}/F_{\mathrm{Kr}})q\sqrt{\mathit{\nu }/\mathit{ϵ}}}}-F_{\mathrm{Kr}}=\mathrm{0}.\end{array}$$

Note that $F_{\mathrm{Kr}}(\mathit{ϵ},{\mathit{\chi }}_{\mathrm{init}},W)$ depends on the underestimate χ_init, whereas it depends on the “true” or “corrected” value of ϵ as calculated in Sect. 4.

6.2 Obtaining χ_init from fast-thermistor voltage spectra

Like we did for ϵ_init in Sect. 4.2, we derive the expression for χ_init in three steps. First, equate the right-hand sides of Eqs. (11) and (29) (excluding the transfer function $H_{Tt}^{\mathrm{2}}\left(k\right)$, which is incorporated into Eq. 31):

$$\begin{array}{}\text{(32)}& {\left({\displaystyle \frac{C_{\mathrm{2}T}+\mathrm{2}{C}_{\mathrm{3}T}\langle V_T\rangle }{C_{Tt}}}\right)}^{\mathrm{2}}{\displaystyle \frac{\mathrm{1}}{W}}{\mathrm{\Psi }}_{Tt}\left(f\right)=\mathrm{4}{\mathit{\pi }}^{\mathrm{2}}k\mathit{\chi }q\sqrt{\mathit{\nu }/\mathit{ϵ}}.\end{array}$$

Then, rearrange while substituting $k=f/W$ to get

$$\begin{array}{}\text{(33)}& {\mathit{\chi }}_{\mathrm{init}}{f}^{\mathrm{1}}={\displaystyle \frac{\mathrm{1}}{\mathrm{4}{\mathit{\pi }}^{\mathrm{2}}q\sqrt{\mathit{\nu }/\mathit{ϵ}}}}{\left({\displaystyle \frac{C_{\mathrm{2}T}+\mathrm{2}{C}_{\mathrm{3}T}\langle V_T\rangle }{C_{Tt}}}\right)}^{\mathrm{2}}{\mathrm{\Psi }}_{Tt}.\end{array}$$

Finally, solve for χ_init using a least-squares fit (Appendix C):

$$\begin{array}{}\text{(34)}& {\mathit{\chi }}_{\mathrm{init}}={\displaystyle \frac{\mathrm{1}}{\mathrm{4}{\mathit{\pi }}^{\mathrm{2}}q\sqrt{\mathit{\nu }/\mathit{ϵ}}}}{\left({\displaystyle \frac{C_{\mathrm{2}T}+\mathrm{2}{C}_{\mathrm{3}T}\langle V_T\rangle }{C_{Tt}}}\right)}^{\mathrm{2}}\underset{\text{Calculate on board}}{\underbrace{{\displaystyle \frac{{\sum }_{f}f{\mathrm{\Psi }}_{Tt}}{{\sum }_f{f}^{\mathrm{2}}}}}}.\end{array}$$

6.3 Quality control of the temperature gradient spectral fits

The approach to quality controlling the fast-thermistor data is the same as that for shear (Sect. 4.3). That is, we fit Ψ_Tt over f_l–f_m and f_m–f_h (1–3 and 3–5 Hz). This ultimately provides two estimates of χ for each spectrum, which are combined as follows:

$$\begin{array}{}\text{(35)}& {\displaystyle }\mathit{\chi }=\mathrm{mean}\left({\mathit{\chi }}_{\text{l–m}},{\mathit{\chi }}_{\text{m–h}}\right),\text{(36)}& {\displaystyle }\mathit{\chi }\text{ fit score}={\displaystyle \frac{min\left({\mathit{\chi }}_{\text{l–m}},{\mathit{\chi }}_{\text{m–h}}\right)}{max\left({\mathit{\chi }}_{\text{l–m}},{\mathit{\chi }}_{\text{m–h}}\right)}}.\end{array}$$

We do not apply a low χ threshold equivalent to Eq. (27).

7 Test of the reduction scheme for χ

Profiles of χ from the reduced scheme compare well to the standard processing, albeit with a small bias in one direction for low values and in the other direction for high values (Fig. 8). Across all values, the two approaches agree within a factor of 2 78 % of the time (Fig. 9). By comparison, 82 % of segments exhibit a factor-of-2 agreement between χ values from the two fast thermistors on the same unit.

Figure 8 Testing the proposed data reduction scheme for fast-thermistor measurements. The profiles used are the same as those chosen in Fig. 5.

Figure 9 Statistical test of the proposed data reduction scheme for χ. Equivalent to Fig. 6, except for the use of χ instead of ϵ. In total, there are 100 000 segments of data.

Compared to shear spectra, non-dimensionalized temperature gradient spectra have lower fit scores. Especially for the lowest fit scores, the measured temperature gradient spectra tend to be too high at lower frequencies and too low at frequencies near f_h (Fig. 10a–c). As frequency increases beyond f_h, the effects of noise and thermal-response corrections (Appendix A) begin to dominate.

Figure 10 Same as for Fig. 7 but for temperature gradient rather than shear.

There are three reasons for the poorer fits to temperature gradient spectra compared to that for shear. First, shapes of temperature gradient spectra are often more variable; the best choice for the non-dimensional spectral model can be debated (e.g., Sanchez et al., 2011). Second, the temperature gradient fits depend on ϵ, so uncertainties in ϵ propagate into the calculation of χ. Third, for our 2019 experiment, the recorded temperature gradient signals were occasionally affected by digitization noise as a consequence of sampling mixed layers. (Shear signals were not affected by digitization noise.)

8 Recommendations

8.1 Setting the scheme's parameters

Our scheme requires a few user-defined parameters: f_l, f_h, N_seg, and N_fft. For this paper, we based these partly on the profiling speed and scientific goals of FCS. For a different profiler, we suggest the following:

• Choose f_h based on a typical profiling speed such that $k_{\mathrm{h}}=f_{\mathrm{h}}/W\approx \mathrm{25}$ cpm for a nominal profiling speed W. For a wide range of ϵ values, 25 cpm is close to, or beyond, the peak of the Nasmyth spectrum (Fig. 1). As an example, if FCS profiled at ∼ 0.5 m s⁻¹, we would consider setting f_h≈12 Hz.

• Ensure that there are no known issues such as vibrations that are likely to adversely affect spectral coefficients within the f_l–f_h range. If, however, there are known issues within the desired frequency range, then an alternative approach (one that we did not test) is to use accelerometer signals to correct spectra that are contaminated by vibrations (Levine and Lueck, 1999; Goodman et al., 2006).

• Define f_l and f_h separately for shear and temperature gradient if appropriate. Although we set them to be equal here, this is not necessary.

• Use more than two fitting bands if desired. We use only two bands (1–3 and 3–5 Hz) so as to minimize the file size to be transmitted, but there is nothing preventing there being three or more bands (e.g., adding a 5–7 Hz band). Indeed, this would enable improved estimates of the fit scores (Eqs. 26 and 36) and more flexibility to discard noise-affected bands as in Eq. (27). If file size is less of an issue such that it is possible to send back fit values for many more than two bands, then the Rainville et al. (2017) scheme outlined in Sect. 1 may be a better choice than ours.

• Choose N_seg and N_fft based on scientific goals and, possibly, any logistical constraints; the data reduction scheme is agnostic to these numbers. For example, at the expense of vertical resolution, we could halve the file size of our transmitted dataset by doubling N_seg from 512 to 1024.

• Reasonable choices for N_fft are $N_{\mathrm{seg}}/\mathrm{2}$ or $N_{\mathrm{seg}}/\mathrm{4}$, which correspond to three or seven half-overlapping subsegments, respectively. There is little to be gained by dividing a segment into even more subsegments so as to produce smoother spectra before fitting. As Ruddick et al. (2000) note, the choice is analogous to fitting a line to 20 points at once or first clumping them into groups of, say, five and then fitting the four averaged points.

8.2 Evaluating the reduced data

One step that cannot be automated is the heuristic evaluation of the reduced turbulence data after they have been converted from voltage quantities to physical ones. For this evaluation, we recommend looking into multiple quantities. First consider the fit scores (Eqs. 26 and 36). We recommend discarding any ϵ or χ values with an associated fit score lower than 0.33. Note, however, that these scores are not a perfect measure of fit. They should be used together with other quality control checks such as comparing

• W and W_min (Eqs. 13 and 14) to check whether the profiling speed is constant over a segment;

• ϵ values from the two shear probes; and

• turbulent features in successive profiles.

The last point is most applicable for a vertical profiler cycling rapidly – for example, twice per hour for FCS. In this case, the profiler is nominally sampling the same vertical fragment of the ocean on a timescale comparable to that over which turbulence evolves. In our experience, many turbulent patches extend over 5–10 profiles.

Recall, also, that all uncertainty in ϵ propagates into the calculation of χ (Sect. 7). If ϵ for a given segment cannot be trusted, neither can χ.

9 Conclusions

We have developed a data reduction scheme applicable to vertical profiling of turbulence variables in which each ∼ 5 s segment is distilled to 12 quantities (Scheme 1). In post-processing, we reconstruct estimates of ϵ and χ, associated quality control metrics, and other quantities such as the temperature and profiling speed. The raw data that go into the 12 quantities are seven different voltages (V_P; V_T and V_Tt for each thermistor; and V_s for each shear probe). Hence, for each 512-element segment, we effectively reduce the data by a factor of $\mathrm{512}\times \mathrm{7}/\mathrm{12}\approx \mathrm{300}$.

Scheme 1 Summary of the data reduction scheme. Each 512-element segment of data is ultimately compressed down to the 12 highlighted quantities that are then transmitted. These are calibrated and/or converted into turbulence quantities in post-processing.

This reduction compresses the output data file size for each dive from megabytes to kilobytes. For example, the total amount of data per dive (two profiles) can be estimated assuming our nominal dive depth and profiling velocity of 120 m and 0.2 m s⁻¹. Each dive creates $\mathrm{12}\times \mathrm{2}\times \mathrm{120}$ m $/$ (0.2 m s${}^{-\mathrm{1}}\times \mathrm{5.12}\mathrm{s})\approx \mathrm{2800}$ quantities. Transmitting each quantity as a 16-bit float or integer equates to approximately 6 kB per dive. This can be reduced by one-third if the spectral-fit metrics are suitably scaled logarithmically and then transmitted as 8-bit integers.

One luxury we lose is the ability to inspect the raw signals. Typically this would help to (i) cultivate faith in the data, (ii) flag which segments to discard, and (iii) inform work-arounds such as filtering out potential narrowband vibrations in shear spectra. Our scheme accounts for this constraint in two ways. First, we fit spectra over relatively low frequencies (1–5 Hz) that are unlikely to be affected by noise or vibration. Second, we reduce the data in a way that uses as little arithmetic as possible. Obviously, we cannot reverse-engineer the raw signals, but by making the onboard calculations simple, we give ourselves the best chance to later fix or identify any unforeseen issues.

Although the onboard reduction eliminates possibilities in how we process turbulence data, it opens up possibilities in how we obtain turbulence data. By visualizing how turbulence evolves over successive dives in near-real time, we can concentrate on regions of interest by adapting the dive schedule to profile more frequently or to different depths. If instead we encounter quiescent periods, we might consider profiling less frequently, thereby conserving battery life. Our ultimate objective is to treat FCS floats as expendable.

Appendix A: Transfer functions for FCS sensors

Voltage signals from shear probes and thermistors are a smoothed representation of the true environmental signal. If the smoothing is a spatial effect, it is described by a transfer function H²(k). If the smoothing is a temporal effect, it is more natural to use H²(f). We can use these interchangeably because f=Wk, and therefore H²(f)=H²(Wk). For FCS, there are three components to the transfer function for each sensor:

$$\begin{array}{}\text{(A1)}& {\displaystyle }{H}_{\mathrm{s}}^{\mathrm{2}}\left(k\right)=H_{\mathrm{SP}}^{\mathrm{2}}\left(k\right)H_{\mathrm{AA}}^{\mathrm{2}}\left(f\right)H_{\mathrm{D}}^{\mathrm{2}}\left(f\right),\text{(A2)}& {\displaystyle }{H}_{Tt}^{\mathrm{2}}\left(k\right)=H_{\mathrm{FT}}^{\mathrm{2}}\left(f\right)H_{\mathrm{AA}}^{\mathrm{2}}\left(f\right)H_{\mathrm{D}}^{\mathrm{2}}\left(f\right),\end{array}$$

where we have used the following shorthand: SP is shear probe, FT is fast thermistor, AA is anti-aliasing, and D is digital. We describe each of these in turn.

Shear probes built and calibrated by the Ocean Mixing Group are very close in dimension to those examined by Ninnis (1984), who measured their wavenumber response and represented it as

$$\begin{array}{}\text{(A3)}& H_{\mathrm{SP}}^{\mathrm{2}}\left(k\right)=\sum _{n=\mathrm{0}}^{\mathrm{4}}{a}_n{\left({\displaystyle \frac{k}{k_{\mathrm{0}}}}\right)}^n,\end{array}$$

where a₀=1.000, $a_{\mathrm{1}}=-\mathrm{0.164}$, $a_{\mathrm{2}}=-\mathrm{4.537}$, a₃=5.503, $a_{\mathrm{4}}=-\mathrm{1.804}$, and k₀=170 cpm.

Temporal averaging of temperature at high frequencies due to the thermal response of the fast thermistor is modeled using a double-pole filter:

$$\begin{array}{}\text{(A4)}& H_{\mathrm{FT}}^{\mathrm{2}}\left(f\right)={\displaystyle \frac{\mathrm{1}}{{\left(\mathrm{1}+(f/f_{\mathrm{c}}{)}^{\mathrm{2}}\right)}^{\mathrm{2}}}},\end{array}$$

where the cut-off frequency f_c=30 Hz. This comes from Nash et al. (1999), who measured the frequency response for two different thermistors on an instrument profiling at 0.3 m s⁻¹ and found cut-off frequencies of 25.1 and 36.7 Hz (see their Fig. A2). The 30 Hz value is the approximate mean of these two values.

Raw shear and thermistor voltage signals are both subject to two filters. First, an analog anti-aliasing filter (two-pole Butterworth) with an f_c=40 Hz cut-off:

$$\begin{array}{}\text{(A5)}& H_{\mathrm{AA}}^{\mathrm{2}}\left(f\right)={\displaystyle \frac{\mathrm{1}}{\mathrm{1}+(f/f_{\mathrm{c}}{)}^{\mathrm{4}}}}.\end{array}$$

After the analog signal is anti-aliased, it is digitized at 400 Hz. Before subsampling to the final 100 Hz output, the signal is digitally filtered. For the 2019 FCS cruise, the signal was convolved with a symmetric 29-element kernel in which the first 15 elements were

$$\begin{array}{}\text{(A6)}& \begin{array}{rl}{g}_i=({\mathrm{2}}^{\mathrm{16}}-\mathrm{1}{)}^{-\mathrm{1}}\times & [\mathrm{52},\mathrm{221},\mathrm{393},\mathrm{427},\mathrm{174},\mathrm{0},\mathrm{0},\mathrm{0},\mathrm{0},\mathrm{0},\\ & \mathrm{1970},\mathrm{5054},\mathrm{8202},\mathrm{10}\mathrm{558},\mathrm{11}\mathrm{433}].\end{array}\end{array}$$

This is a sinc kernel but with negative values set to zero. (We are currently investigating better choices for future implementations). The filter has a half-power (−3 db) point at 25 Hz.

Appendix B: Identifying the start and end of a profile

Early in our processing routine, we partition the raw voltage signals into 512-element segments. In order to discard the segments in which FCS was not profiling, we need robust (yet simple) criteria that demarcate the start and end of a profile. For the start, we search for the first three consecutive segments in which W_min>0.05 m s⁻¹. For the end, we swap the inequality.

A drawback of this approach is the appearance of a quantity in physical units (0.05 m s⁻¹). This is the one instance where we hard code a calibration coefficient in the onboard software rather than applying it in post-processing. Fortunately, the relevant coefficient can be approximated as constant: C_2P = 76.7 psi V⁻¹ (barring a redesign of the circuitry or the use of a different brand or model of pressure sensor). For the two units already built, C_2P=76.81 psi V⁻¹ and 76.53 psi V⁻¹. By comparison, among the four shear probes on the two units, the calibration coefficients vary by 30 %.

At least for the initial implementation of our scheme, we do not include an algorithm to detect the surface to within centimeters. Doing so would let us work backward to put our uppermost depth bin as close to the surface as possible. However, we expect that this could be a fragile part of the scheme. Further, FCS lacks a micro-conductivity sensor, which is likely the sensor best suited for identifying the air–sea interface (e.g., Ward et al., 2014).

Without surface detection, the depths of the uppermost bins will be realized randomly. In the worst cases, we would discard the top ∼ 1 m (5 s at ∼ 0.2 m s⁻¹). To alleviate this, we may use half-overlapping bins near the surface. The exact implementation will be determined later in the development.

Appendix C: Least-squares fitting of power laws

In this paper, we use power law fits to derive turbulence quantities: ${\mathrm{\Psi }}_{\mathrm{s}}=A_{\mathit{ϵ}}{f}^{\mathrm{1}/\mathrm{3}}$ and Ψ_Tt=A_χf¹, where A_ϵ and A_χ are substitutes for the expressions in Eqs. (23) and (33). With only a single parameter for each fit, implementing a least-squares fit is easy.

Assume we are fitting the vector Ψ_i to the function $A{f}_i^n$ where n is either $\mathrm{1}/\mathrm{3}$ or 1. The sum of squared residuals is therefore

$$\begin{array}{}\text{(C1)}& \sum r_i^{\mathrm{2}}=\sum {\left({\mathrm{\Psi }}_i-A{f}_i^n\right)}^{\mathrm{2}}.\end{array}$$

The minimum with respect to A is where the derivative is zero:

$$\begin{array}{}\text{(C2)}& {\displaystyle \frac{\partial }{\partial A}}\sum r_i^{\mathrm{2}}=\sum -\mathrm{2}{f}_i^n\left({\mathrm{\Psi }}_i-A{f}_i^n\right)=\mathrm{0}.\end{array}$$

Hence,

$$\begin{array}{}\text{(C3)}& A={\displaystyle \frac{\sum f_i^n{\mathrm{\Psi }}_i}{\sum f_i^{\mathrm{2}n}}}.\end{array}$$

We had originally intended to find A by following Becherer and Moum (2017), who were fitting $f^{\mathrm{1}/\mathrm{3}}$ spectra. Their simpler method, $A=\sum \left({\mathrm{\Psi }}_i/f_i^n\right)$, is equivalent to a least-squares fit, except the quantity minimized is the sum of the squares of the adjusted residuals, where adjusted means divided by fⁿ. Differences can be ignored when $n=\mathrm{1}/\mathrm{3}$ but not when n=1.

Appendix D: Standard processing of FCS turbulence measurements

The standard processing of FCS turbulence data differs from the reduced scheme in three ways. First, raw data are despiked differently (Appendix E). Second, the 100 Hz raw voltage signals are calibrated into physical quantities right away. Hence, means and spectra are calculated in physical units and not in voltage units. Third, the integration of spectra occurs over a variable wavenumber band, which is found iteratively.

When integrating shear spectra (after correction; see Appendix A) to find ϵ, we follow the approach used for the Chameleon profiler (Moum et al., 1995). A first estimate of ϵ is made by integrating over k=4–10 cpm. This value provides a first estimate of the Kolmogorov wavenumber $k_{\mathrm{s}}=(\mathit{ϵ}/{\mathit{\nu }}^{\mathrm{3}}{)}^{\mathrm{1}/\mathrm{4}}/\mathrm{2}\mathit{\pi }$. (The lower limit for Chameleon is 2 cpm, but we increase this for FCS given its slower profiling speed and hence the possibility of contamination by waves at lower wavenumbers.) The upper integral limit is then set to 0.5k_s (with a minimum of 10 cpm and a maximum of 45 cpm). The Nasmyth spectra (Eq. 15) is integrated over the same wavenumber range. If the measured and Nasmyth integrals are within 1 %, then ϵ is set to be equal to the integral of the Nasmyth spectrum over all k. Otherwise, ϵ and k_s are adjusted iteratively until the two integrals agree.

A similar but non-iterative approach is used for integrating T_z spectra to find χ. The model spectrum is the Kraichnan spectrum (Eq. 28), and again, the lower limit of integration is 4 cpm. The upper limit is the Batchelor wavenumber $k_b=(\mathit{ϵ}/\mathit{\nu }{D}_{\mathrm{T}}^{\mathrm{2}}{)}^{\mathrm{1}/\mathrm{4}}/\mathrm{2}\mathit{\pi }$ (with a maximum defined by kW=15 Hz).

Appendix E: Identifying and removing noise and spikes in the shear signals

To properly despike the raw output of a shear probe requires several steps. Lueck et al. (2018) describe a process in which the signal is high-passed, then rectified, and then low-passed to derive a measure of the local variance. A value is defined as a spike if it is more than 8 times (or a similar threshold) higher than the local variance. Spikes are replaced with an average based on surrounding points. This process is then repeated on the new signal and so on until no spikes are identified.

In our standard processing of FCS data, we use the Lueck et al. (2018) despiking routine. For our data reduction scheme, we use an approach that is easier to implement and quicker to compute, albeit one that is less precise. For each 512-element segment of data, a spike is defined as any data point larger than 3 standard deviations from the mean. These spikes are replaced by the mean of the remaining values in the segment.